c begin file flsh_sub.f c c This file contains iterative flash routines which call the c intermediate-level routines c c contained here are: c subroutine TPRHO (t,p,x,kph,kguess,rho,ierr,herr) c subroutine TPFLSH (t,p,z,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr) c subroutine TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr) c subroutine TDFLSH (t,D,x,p,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,herr) c subroutine TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew, c & p,Dl,Dv,x,y,q,ierr,herr) c subroutine PHFLSH (p,h,x,t,D,Dl,Dv,xl,xv,q,e,s,cv,cp,w,ierr,herr) c subroutine PHFL1 (p,h,x,kph,t,D,ierr,herr) c subroutine PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, c & t,Dl,Dv,x,y,q,ierr,herr) c subroutine PSFLSH (p,s,x,t,D,Dl,Dv,xl,xv,q,e,h,cv,cp,w,ierr,herr) c subroutine PSFL1 (p,s,x,kph,t,D,ierr,herr) c subroutine PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, c & t,Dl,Dv,x,y,q,ierr,herr) c c these routines use the following common blocks from other files c common /Gcnst/ R c common /NCOMP/ nc c common /CHAR/ htab,hnull c c various arrays are dimensioned with parameter statements c parameter (ncmax=5) !max number of components in mixture c parameter (n0=-ncmax,nx=ncmax) c c ====================================================================== c ====================================================================== c subroutine TPRHO (t,p,x,kph,kguess,rho,ierr,herr) c c iterate for density as a function of temperature, pressure, and c composition for a specified phase c c inputs: c t--temperature [K] c p--pressure [kPa] c x--composition [array of mol frac] c kph--phase flag: 1 = liquid c 2 = vapor c N.B.: 0 = stable phase--NOT ALLOWED (use TPFLSH) c kguess--input flag: 1 = first guess for rho provided c 0 = no first guess provided c rho--first guess for molar density [mol/L], only if kguess = 1 c c outputs: c rho--molar density [mol/L] c ierr--error flag: 0 = successful c 200 = CRITP did not converge c 201 = illegal input (kph <= 0) c 202 = liquid-phase iteration did not converge c 203 = vapor-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 04-03-95 MM, original version c 09-11-95 MM, add error string to argument list c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 01-22-95 MM, check that input rho > 0 before using as initial guess c 01-23-95 MM, use Rackett technique for liquid density initial guess c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 03-07-96 MM, change tolr from 1.0d-6 to 1.0d-8 c 10-16-96 MM, change call from DPRHO to DPDD c 12-02-96 MM, add check that vapor density guess does not give p < 0 c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 02-11-97 EWL, reduced tolerance after 10 and after 15 iterations c 04-22-97 MM, add additional iteration for low-p (allow p = 0) c 05-16-97 MM, do not apply step limitation for t > 1.5*tc c 06-03-97 EWL, add second order Newton's solution to increase convergence speed c 06-05-97 EWL, remove checks in liquid solution for large jumps in density c change initial guess in liquid phase to 2*Dc c 07-15-97 MM, renumber and add detail to error messages c 09-22-97 EWL, modify test to select liquid or vapor iteration c 09-23-97 MM, increase itmax from 20 to 30 to aid in near-critical convergence c 10-01-97 MM, add compiler switches to allow access by DLL c 11-13-97 EWL, revert to 1st order Newton's method if 2nd order gives bad result c 11-18-97 EWL, ditto for liquid-phase iteration c 12-31-97 MM, check derivative dP/dD for liquid-phase initial guess c 02-11-98 MM, check for t > Tc before switching liq to vapor iteration c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: TPRHO c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 dimension x(ncmax) common /Gcnst/ R common /NCOMP/ nc common /CHAR/ htab,hnull c common /IT/ it c !temporary--convergence tolerance, put in common block data itmax /30/ tolr = 1.0d-8 c !temporary--end ierr=0 herr=hnull c write (*,1001) t,p,(x(i),i=1,5) c1001 format (1x,' TPRHO--input t,p,x: ',f10.5,f14.6,5f10.6) c c determine which iteration to use: c For liquids (or fluids above the critical pressure) the iteration is c carried out in transformed coordinates of log(V). For vapor (or c fluids at supercritical temperatures but pressures below the critical c value) the iteration is in terms of log(V) and log(p). For very low c pressure vapors a simple iteration in rho and p is used. The iteration c has converged when the pressure calculated from the equation of state c agrees with the input pressure within a convergence tolerance. c if (p.lt.1.0d-14) then c ideal-gas rho=p/(R*t) RETURN end if c call CRITP (x,tc,pc,rhoc,ierr,herr1) if (ierr.ne.0) then ierr=200 write (herr,1200) herr1(1:236),hnull 1200 format ('[TPRHO error 200] ',a236,a1) call ERRMSG (ierr,herr) RETURN end if vc=1.0d0/rhoc vclog=LOG(vc) c EWL modification to aid near-critical convergence; c following approximates the pressure along the critical isochore if (p.lt.pc*(1.0d0+7.0d0*(t/tc-1.0d0)) .and. t.gt.tc) then c use 'vapor phase' iteration lliq=.false. else if (p.gt.pc .or. t.gt.tc) then c use 'liquid phase' iteration lliq=.true. c EWL modification; change from vapor to liquid iteration c else if (t.gt.tc) then c use 'vapor phase' iteration c lliq=.false. else if (kph.le.0) then c illegal input to TPRHO (two-phase state) ierr=201 herr='[TPRHO error 201] illegal input to TPRHO (kph <= 0); ' & //'use TPFLSH instead'//hnull RETURN else c use input phase specification to find solution if (kph.eq.1) then lliq=.true. else lliq=.false. end if end if c c set initial guess for density 10 if (kguess.eq.1 .and. rho.gt.0.0d0) then c use input density as initial guess vlog=LOG(1.0d0/rho) else c initial guess based on region if (p.gt.pc .and. t.gt.tc) then vlog=LOG(0.5d0*vc) else if (lliq) then c vlog=LOG(0.5d0*vc) c initial guess for sub-critical liquid density based on modified c Rackett technique (Reid, Prausnitz, and Poling (1987), Properties c of Gases and Liquids, 4th edition) except that Zra = Zc Rtp=R*tc/pc Zc=vc/Rtp c vlog=LOG(Rtp*Zc**(1.0d0+(1.0d0-t/tc)**(2.0d0/7.0d0))) Dguess=1.0d0/(Rtp*Zc**(1.0d0+(1.0d0-t/tc)**(2.0d0/7.0d0))) iguess=1 80 call DPDD (t,Dguess,x,dpdrho) if (dpdrho.le.0.0d0 .and. iguess.le.8) then c initial guess is in two-phase region, set to higher density Dguess=1.1d0*Dguess iguess=iguess+1 goto 80 end if vlog=LOG(1.0d0/Dguess) if (vlog.gt.LOG(0.5d0*vc)) vlog=LOG(0.5d0*vc) else vlog=LOG(R*t/p) end if if (.not.lliq) vlog=LOG(R*t/p) end if c write (*,1082) kph,kguess,t,p,1.0/EXP(vlog) c1082 format (1x,' TPRHO--kph,kguess,t,p,initial rho: ',2i4,3e14.6) c c enter Newton's method iteration, separate loops for liquid and vapor c rho2=0.0d0 if (lliq) then c liquid phase iteration do 100 it=1,itmax vlog0=vlog c write (*,*) ' TPRHO (liquid) vlog: ',vlog rho=1.0d0/exp(vlog) call PRESS (t,rho,x,p2) call DPDD (t,rho,x,dpd) call DPDD2 (t,rho,x,dpd2) c write (*,1010) it,t,x(1),p,p2,dpd,rho c1010 format (1x,'it,t,x(1),p,p2,dpd,rho (liquid): ',i3,2f8.3,4e16.8) if (dpd.lt.0.0d0) then c unstable portion of two-phase region, make another guess c write (*,1012) it,1.0/exp(vlog) c1012 format (1x,'unstable 2-phase in liquid, it, rho = ',i4,d16.8) vlog=vlog-0.10d0 else dpd2=dpd2+dpd/rho !2nd order Newton's method fvdp=(p-p2)/(dpd2*(p-p2)/2.0d0/dpd+dpd)/rho if (ABS(fvdp).gt.100.0d0) then c revert to 1st order Newton's method if too large step from 2nd order dpdlv=-rho*dpd fvdp=(p2-p)/dpdlv !1st order Newton's method endif c write (*,1014) dpdlv,dpd,fvdp c1014 format (1x,' TPRHO (liquid) dpdlv,dpd,fvdp: ',3d16.8) c if calculation has not converged after 10 or 15 iterations, loosen c tolerance (original tolr was too tight near critical) if (it.eq.10) tolr=tolr*10 if (it.eq.15) tolr=tolr*10 c The rho-rho2 check is only important for very low pressures on the c liquid surface (propane or R124). In this case, the change in density c required to get the correct pressure is less than machine precision. if (ABS(fvdp/p).lt.0.0001*tolr .or. & ABS(rho-rho2).lt.1.0d-11) then c iteration has converged rho=1.0d0/exp(vlog-fvdp) c call PRESS (t,rho,x,p2) !convergence testing only c write (*,1016) it,t,p,p2,rho c1016 format (1x,8x,'% TPRHO converged (liq); it,t,p,p2,rho ', c & i3,f8.3,2e18.10,e18.10) c write (*,*) ' TPRHO final rho-liq: ',rho RETURN else c next guess vlog=vlog0-fvdp c do not allow too great of a change in density between iterations c for liquid states if (ABS(vlog-vlog0).gt.0.1d0 .and. t.lt.1.5d0*tc) then vlog=vlog0+SIGN(0.1d0,vlog-vlog0) end if if (vlog.gt.vclog .and. t.lt.tc) then vlog=0.5d0*(vlog0+vclog) end if c switch to the "vapor" iteration if (vlog.lt.-5.0d0 .and. t.gt.tc) then c write (*,*) ' TPRHO--switching to the vapor iteration' lliq=.false. goto 10 end if end if rho2=rho end if 100 continue c iteration has not converged rho=1.0d0/exp(vlog) ierr=202 write (herr1,1202) t,p,rho,(x(i),i=1,nc) 1202 format ('[TPRHO error 202] liquid iteration has not ', & 'converged for T =',1pe11.3,' K, P =',1pe11.3, & ' kPa, rho (last guess) = ',1pe11.3,' mol/L,', & ' compositions = ',0p5f9.5) herr=herr1(1:193)//hnull call ERRMSG (ierr,herr) c write (*,*) ' TPRHO final rho-liq (not converged): ',rho RETURN c c else if (lowp) then cc low-pressure vapor iteration c disabled--not required c do 200 it=1,itmax c call PRESS (t,rho,x,p2) c call DPDD (t,rho,x,dpd) c delp=p2-p c if (ABS(delp).lt.0.0001*tolr) then cc iteration has converged cc write (*,1024) it,t,p,p2,delp,rho cc1024 format (9x,'% TPRHO converged (low-p); it,t,p,p2,delp,rho ', cc & i3,f8.3,3e18.10,e18.10) c RETURN c else cc next guess c rho=rho-delp/dpd c end if c 200 continue cc iteration has not converged c rho=p/(R*t) c ierr=1 c write (herr,1002) t,p,rho,hnull c 1002 format (' ERROR--TPRHO has not converged (low-p), t,p,rho =', c & f8.3,2e14.6,a1) c RETURN cc else c vapor phase iteration plog=LOG(p) do 300 it=1,itmax vlog0=vlog rho=1.0d0/exp(vlog) call PRESS (t,rho,x,p2) call DPDD (t,rho,x,dpd) call DPDD2 (t,rho,x,dpd2) c write (*,1030) it,t,x(1),p,p2,dpd,rho c1030 format (1x,'it,t,x(1),p,p2,dpd,rho (vapor): ',i3,2f8.3,4e18.10) if (dpd.lt.0.0d0 .or. p2.le.0.0d0) then c unstable portion of two-phase region or p2<0, make another guess c write (*,1032) it,1.0/exp(vlog) c1032 format (1x,'unstable 2-phase in vapor or p2 < 0, it, rho = ', c & i4,d16.8) vlog=vlog+0.10d0 else dpd2=dpd2/p2-(dpd/p2)**2+dpd/rho/p2 fvdpl=(plog-LOG(p2))/(p2*dpd2*(plog-LOG(p2)) & /2.0d0/dpd+dpd/p2)/rho !2nd order Newton's method c if 2nd order Newton's method gives unreasonable result, revert to c first order method (large value of fvdpl corresponds to huge change c in next guess for volume) if (ABS(fvdpl).gt.1.0d0) then dpdlv=-rho*dpd fvdpl=(LOG(p2)-plog)*p2/dpdlv !1st order Newton's method endif if (ABS(fvdpl).gt.1.0d0) fvdpl=0.1d-5 c if calculation has not converged after 10 or 15 iterations, loosen c tolerance (original tolr was too tight near critical) if (it.eq.10) tolr=tolr*10 if (it.eq.15) tolr=tolr*10 if (ABS(fvdpl).lt.0.0001*tolr) then c iteration has converged rho=1.0d0/exp(vlog) call PRESS (t,rho,x,p2) !convergence testing only c write (*,1034) it,t,p,p2,rho c1034 format (1x,8x,'% TPRHO converged (vap); it,t,p,p2,rho ', c & i3,f8.3,2e18.10,e18.10) c write (*,*) ' TPRHO final rho-vap: ',rho RETURN else c next guess vlog=vlog-fvdpl c if (ABS(vlog-vlog0).gt.0.5d0) then c vlog=vlog0+SIGN(0.5d0,vlog-vlog0) c end if c if (vlog.lt.vclog .and. t.lt.tc) then c vlog=0.5d0*(vlog0+vclog) c end if end if end if 300 continue c iteration has not converged rho=1.0d0/exp(vlog) c iteration has not converged rho=1.0d0/exp(vlog) ierr=203 write (herr,1203) t,p,rho,hnull 1203 format ('[TPRHO error 203] vapor iteration has not ', & 'converged for T =',1pe11.3,' K, P =',1pe11.3, & ' kPa, rho (last guess) = ',1pe11.3,' mol/L',a1) call ERRMSG (ierr,herr) c write (*,*) ' TPRHO final rho-vap (not converged): ',rho RETURN c end if c end !subroutine TPRHO c c ====================================================================== c subroutine TPFLSH (t,p,z,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr) c c flash calculation given temperature, pressure, and bulk composition c c This routine accepts both single-phase and two-phase states as the c input; if the phase is known, the subroutine TPRHO is faster. c c inputs: c t--temperature [K] c p--pressure [kPa] c z--overall (bulk) composition [array of mol frac] c c outputs: c D--overall (bulk) molar density [mol/L] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c if only one phase is present, Dl = Dv = D c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c if only one phase is present, x = y = z c q--vapor quality on a MOLAR basis [moles vapor/total moles] c q < 0 indicates subcooled (compressed) liquid c q = 0 indicates saturated liquid c q = 1 indicates saturated vapor c q > 1 indicates superheated vapor c q = 998 superheated vapor, but quality not defined (t > Tc) c q = 999 indicates supercritical state (t > Tc) and (p > Pc) c e--overall (bulk) internal energy [J/mol] c h--overall (bulk) enthalpy [J/mol] c s--overall (bulk) entropy [J/mol-K] c Cv--isochoric (constant V) heat capacity [J/mol-K] c Cp--isobaric (constant p) heat capacity [J/mol-K] c w--speed of sound [m/s] c Cp, w are not defined for 2-phase states c in such cases, a flag = -9.99998d6 is returned c ierr--error flag: 0 = successful c 1 = Tin < Tmin c 4 = Pin < 0 c 5 = T and P out of range c 8 = x out of range (component and/or sum < 0 c or > 1) c 9 = x and T out of range c 12 = x out of range and P < 0 c 13 = x and T and P out of range c 210 = CRITP did not converge c 211 = SATT did not converge at bubble point c 212 = SATT did not converge at dew point c 213 = TPRHO did not converge for liquid state c 214 = TPRHO did not converge for vapor state c 215 = TPRHO did not convg for supercritical state c 216 = TPRHO did not convg for liq in 2-phase it c 217 = TPRHO did not convg for vap in 2-phase it c 218 = liquid frac for 2-phase it did not converge c 219 = composition for 2-phase it did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-11-95 MM, original version c 09-12-95 MM, change argument list (add e, h, s, etc.) c 09-20-95 MM, define output xl, xv for single-component super-critical c 09-25-95 MM, rearrange argument list (outputs in order rho, x, q) c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 12-12-95 MM, restructure logic to accomodate mixtures c first guess for 2-phase mix from ratio of dew & bubble c 12-13-95 MM, set undefined output compositions to zero c initial rho guess for 1-phase vapor based on p/pdew c 12-14-95 MM, bug on call to CRITP: pass z (not x) c 12-28-95 MM, add full 2-phase mixture iteration using fugacities c 12-29-95 MM, move 2-phase iteration to separate routine TPFL2 c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 05-31-96 MM, check input t,p,z against limits c 06-03-96 MM, add 'EOS' to calling list for LIMITX c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 04-09-97 MM, initialize x,y to z c 04-22-97 MM, use h, rather than s, to compute q (h converges at p = 0) c 05-15-97 MM, use V to compute q for comp. liq. (T,h can be double-valued) c add q = 998 case for t > Tc, but p < Pc c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c 02-11-98 MM, check that computed density is within bounds c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: TPFLSH c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1,herr2,herrl dimension z(ncmax),x(ncmax),y(ncmax) dimension xdew(ncmax),ydew(ncmax),ybub(ncmax),xbub(ncmax) common /Gcnst/ R common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc c ierr=0 herr=hnull c c initialize output liquid and vapor compositions to input values c zero output values for undefined components do 100 i=1,nc x(i)=z(i) y(i)=z(i) 100 continue if (nc.lt.ncmax) then do 120 i=nc+1,ncmax x(i)=0.0d0 y(i)=0.0d0 120 continue end if c c check that input conditions are within limits c rhodum=0.0d0 call LIMITX ('EOS',t,rhodum,p,z,tmin,tmax,rhomax,pmax,ierrl,herrl) if (ierrl.lt.0) then c one or inputs are outside limits--if just a warning proceed w/ calc write (herr,1006) ierrl,herrl(1:234),hnull 1006 format ('[TPFLSH warning',i3,'] ',a234,a1) call ERRMSG (ierrl,herr) c if error (as opposed to warning) set output density to zero and return else if (ierrl.gt.0) then write (herr,1007) ierrl,herrl(1:236),hnull 1007 format ('[TPFLSH error',i3,'] ',a236,a1) call ERRMSG (ierrl,herr) D=0.0d0 Dl=0.0d0 Dv=0.0d0 call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt) q=999.d0 !quality undefined ierr=ierrl RETURN end if c call CRITP (z,tc,pc,rhoc,ierr,herr2) if (ierr.ne.0) then ierr=210 write (herr,1008) ierr,herr2(1:235),hnull 1008 format ('[TPFLSH error 210] ',a235,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c if (t.ge.tc) then c supercritical state (x = y = z as set above) c write (*,*) ' TPFLSH--supercritical' kph=2 kguess=0 call TPRHO (t,p,z,kph,kguess,D,ierr,herr2) call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt) Dl=D Dv=D if (p.le.pc) then q=998.d0 !superheated, but cannot define quality else q=999.d0 !quality undefined end if if (ierr.ne.0) then ierr=215 write (herr,1215) herr2(1:183),hnull 1215 format ('[TPFLSH error 215] supercritical density iteration ', & 'did not converge: ',a183,a1) call ERRMSG (ierr,herr) RETURN end if c do 130 i=1,nc c x(i)=z(i) c y(i)=z(i) c 130 continue else c subcritical state--call saturation routine to determine liq or vap call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,ydew,ierr,herr2) if (ierr.ne.0) then ierr=212 write (herr,1212) herr2(1:193),hnull 1212 format ('[TPFLSH error 212] dew point calculation ', & 'did not converge: ',a193,a1) call ERRMSG (ierr,herr) D=0.0d0 Dl=0.0d0 Dv=0.0d0 q=999.d0 !quality undefined RETURN end if c call ENTRO (t,Dldew,xdew,sldew) c call ENTRO (t,Dvdew,z,svdew) call ENTHAL (t,Dldew,xdew,hldew) call ENTHAL (t,Dvdew,z,hvdew) if (p.lt.pdew) then c single-phase vapor (x = y = z as set above) c write (*,*) ' TPFLSH--single-phase vapor' kph=2 kguess=1 D=Dvdew*p/pdew !initial guess for density call TPRHO (t,p,z,kph,kguess,D,ierr,herr2) Dl=D Dv=D if (ierr.ne.0) then ierr=214 write (herr,1214) herr2(1:190),hnull 1214 format ('[TPFLSH error 214] vapor density iteration ', & 'did not converge: ',a190,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt) q=(h-hldew)/(hvdew-hldew) c do 140 i=1,nc c x(i)=z(i) c y(i)=z(i) c 140 continue else if (nc.eq.1) then c special case: pure component single-phase liquid c write (*,*) ' TPFLSH--pure fluid single-phase liquid' kph=1 kguess=1 D=Dldew call TPRHO (t,p,z,kph,kguess,D,ierr,herr2) Dl=D Dv=D if (ierr.ne.0) then ierr=213 write (herr,1213) herr2(1:189),hnull 1213 format ('[TPFLSH error 213] liquid density iteration ', & 'did not converge: ',a189,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt) c q=(h-hlbub)/(hvdew-hlbub) c compute quality based on volumes (possible for T,h to be c double-valued in compressed liquid) q=(1.0d0/D-1.0d0/Dldew)/(1.0d0/Dvdew-1.0d0/Dldew) c write (*,*) ' TPFLSH--comp liq q by volumes: ',q c do 150 i=1,nc c x(i)=z(i) c y(i)=z(i) c 150 continue else c mixture: calculate bubble point to determine if liquid or 2-phase call SATT (t,z,1,pbub,Dlbub,Dvbub,xbub,ybub,ierr,herr2) if (ierr.ne.0) then ierr=211 write (herr,1211) herr2(1:190),hnull 1211 format ('[TPFLSH error 211] bubble point calculation ', & 'did not converge: ',a190,a1) call ERRMSG (ierr,herr) D=0.0d0 Dl=0.0d0 Dv=0.0d0 q=999.d0 !quality undefined RETURN end if c call ENTRO (t,Dlbub,z,slbub) call ENTHAL (t,Dlbub,z,hlbub) if (p.ge.pbub) then c mixture single-phase liquid c write (*,*) ' TPFLSH--mixture single-phase liquid' kph=1 kguess=1 D=Dlbub call TPRHO (t,p,z,kph,kguess,D,ierr,herr2) Dl=D Dv=D if (ierr.ne.0) then ierr=213 write (herr,1213) herr2(1:189),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt) c q=(h-hlbub)/(hvdew-hlbub) c compute quality based on volumes (possible for T,h to be c double-valued in compressed liquid q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub) c write (*,*) ' TPFLSH--comp liq q by volumes: ',q c do 160 i=1,nc c x(i)=z(i) c y(i)=z(i) c 160 continue else c c two-phase mixture c write (*,*) ' TPFLSH--two-phase mixture' c c generate initial guesses for x and y by interpolating input pressure c with dew and bubble point pressures c if (ABS(1.0d0-pbub/pdew).lt.1.0d-6) then c possible azeotrope q=0.5d0 else q=1.0d0-(p-pdew)/(pbub-pdew) end if xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 240 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 240 continue do 244 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum 244 continue Dl=Dlbub !initial guess for liquid density Dv=Dvdew !initial guess for vapor density call TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr) if (ierr.ne.0) then c two-phase iteration did not converge--error message written by TPFL2 q=999.d0 !quality undefined RETURN end if c compute 2-phase properties and load output variables c call TPRHO to ensure consistency of t, p, rho kguess=1 call TPRHO (t,p,x,1,kguess,Dl,ierr1,herr1) call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt) call TPRHO (t,p,y,2,kguess,Dv,ierr1,herr1) call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt) if (ierr1.ne.0) then herr2=' ERROR from TPFLSH (mix, 2-phase): '//herr1 call ERRMSG (ierr1,herr2) RETURN end if cp=xnotd !Cp, w not defined for 2-phase w=xnotd cv=xnotd c bulk properties are weighted average of liquid and vapor phases alpha=1.0d0-q !alpha is liq fraction, D=1.0d0/(alpha/Dl+q/Dv) !q is vapor frac e=alpha*el+q*ev h=alpha*hl+q*hv s=alpha*sl+q*sv end if end if end if end if c c if limits check resulted in warning (as opposed to error) return c that message; do this again in case intermediate iteration did not c converge (thereby overwriting any warning message from LIMITX) c if (ierrl.lt.0) then ierr=ierrl write (herr,2006) ierrl,herrl(1:234),hnull 2006 format ('[TPFLSH warning',i3,'] ',a234,a1) call ERRMSG (ierr,herr) end if c if (D.gt.rhomax) then c check that computed (output) density is within bounds of EOS ierr=2 write (herr,2008) ierr,hnull 2008 format ('[TPFLSH error',i3,']; input T,P correspond to ', & 'a density above the limit of the EOS.',a1) call ERRMSG (ierr,herr) end if c RETURN end !subroutine TPFLSH c c ====================================================================== c subroutine TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr) c c flash calculation given temperature, pressure, and bulk composition c c This routine accepts only two-phase states as input; if the phase is c not known use TPFLSH. Use TPRHO for single-phase states. c c inputs: c t--temperature [K] c p--pressure [kPa] c z--overall (bulk) composition [array of mol frac] c Dl--initial guess for molar density [mol/L] of the liquid phase c Dv--initial guess for molar density [mol/L] of the vapor phase c x--initial guess for composition of liquid phase c [array of mol frac] c y--initial guess for composition of vapor phase c [array of mol frac] c q--initial guess for vapor quality on a MOLAR basis c [moles vapor/total moles] c c outputs: c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c q--vapor quality on a MOLAR basis [moles vapor/total moles] c ierr--error flag: (these are also passed up to TPFLSH for output) c 0 = successful c 216 = TPRHO did not converge for liquid c 217 = TPRHO did not converge for vapor c 218 = inner loop (liquid frac) did not converge c 219 = outer loop (composition) did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 12-29-95 MM, original version, extracted from TPFLSH c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 04-09-97 MM, fix bug: converged value of q not returned c 04-11-97 MM, fix divide by zero on saturation boundary c adjust tolr --> d-7; allow alpha outside 0,1 by d-5 c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: TPFL2 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr2 dimension z(ncmax),x(ncmax),y(ncmax) dimension alpha(3),falpha(2),fliq(ncmax),fvap(ncmax),xk(ncmax) common /NCOMP/ nc common /CHAR/ htab,hnull c !temporary--convergence tolerance, put in common block data tolr /1.0d-7/ data itmax /20/ c !temporary--end c ierr=0 herr=hnull c c begin outer iteration loop for the vapor compostion c alpha(1)=1.0d0-q kguess=1 j=1 do 800 iy=1,itmax c compute densities and fugacities for each phase c write (*,1001) t,p,(x(i),i=1,5) c1001 format (1x,' TPFL2 call TPRHO; liq t,p,x: ',f10.5,f14.6,5f10.6) call TPRHO (t,p,x,1,kguess,Dl,ierr,herr2) if (ierr.ne.0) then ierr=216 write (herr,1216) herr2(1:170),hnull 1216 format ('[TPFLSH error 216] liquid density iteration for ', & '2-phase state did not converge: ',a170,a1) call ERRMSG (ierr,herr) RETURN end if call FGCTY (t,Dl,x,fliq) c write (*,1002) t,p,(y(i),i=1,5) c1002 format (1x,' TPFL2 call TPRHO; vap t,p,x: ',f10.5,f14.6,5f10.6) call TPRHO (t,p,y,2,kguess,Dv,ierr,herr) if (ierr.ne.0) then ierr=217 write (herr,1217) herr2(1:170),hnull 1217 format ('[TPFLSH error 217] vapor density iteration for ', & '2-phase state did not converge: ',a170,a1) call ERRMSG (ierr,herr) RETURN end if call FGCTY (t,Dv,y,fvap) c compute k-factors (assume constant for a given guess of y) do 520 i=1,nc if (z(i).gt.0.0d0) then xk(i)=fliq(i)*y(i)/(x(i)*fvap(i)) else c zero composition for component i xk(i)=0.0d0 end if 520 continue c write (*,1004) iy,(xk(i),i=1,nc) c1004 format (1x,' TPFL2 ity,k-factors: ',i4,5e18.10) c c begin inner iteration loop for the liquid fraction (alpha) c j=1 !iteration flag c lfapos=.false. !flags for reguli-falsi iteration c lfaneg=.false. do 600 ix=1,itmax c write (*,1040) ix,(x(i),i=1,nc) c1040 format (1x,' TPFL2 itx,old liq comps: ',i4,10x,5f16.12) xsum=0.0d0 do 540 i=1,nc if (z(i).gt.0.0d0) then c compute new guess for liquid compostion x(i)=z(i)/(alpha(j)+xk(i)*(1.0d0-alpha(j))) else c zero composition for component i x(i)=0.0d0 end if xsum=xsum+x(i) 540 continue c normalize liquid composition do 560 i=1,nc x(i)=x(i)/xsum 560 continue c write (*,1050) ix,xsum,(x(i),i=1,nc) c1050 format (1x,' TPFL2 itx,xsum,new comps: ',i4,6f16.12) falpha(j)=1.0d0-xsum c write (*,1056) ix,alpha(j),xsum,(x(i),i=1,3) c1056 format (1x,' TPFL2:',i4,35x,2f10.6,3f10.6) if (ABS(falpha(j)).lt.1.0d-2*tolr) then c inner loop has converged c (tolerance on inner loop must be tighter than outer loop) goto 640 else c store variables for possible use in reguli-falsi c reguli-falsi actually seems to slow convergence--comment out c if (falpha(j).lt.0.0d0) then c lfaneg=.true. c aneg=alpha(j) c faneg=falpha(j) c else c lfapos=.true. c apos=alpha(j) c fapos=falpha(j) end if c generate next guess for liquid fraction if (j.eq.1) then c for second iteration, move 5 % of way to alpha = 1.0 c alpha(2)=alpha(1)+0.05d0*(1.0d0-alpha(1)) c above could result in alpha(2) = alpha(1) if alpha(1) = 1 c try moving towards alpha = 0.5 !MM 04-11-97 j=2 alpha(2)=alpha(1)+SIGN(0.05d0,0.5d0-alpha(1)) else c for subsequent iterations, use secant method if (ABS(falpha(2)-falpha(1)).gt.1.0d-2*tolr) then alpha(3)=alpha(2)-falpha(2)*(alpha(2)-alpha(1))/ & (falpha(2)-falpha(1)) else alpha(3)=alpha(1)+SIGN(tolr,0.5d0-alpha(1)) end if c new guess for alpha must lie within bounds of zero and one c but allow small tolerance for slop in TPRHO and thus fugacities c write (*,1060) ix,(alpha(i),i=1,3),(falpha(i),i=1,2) c1060 format (1x,' TPFL2 itx,alpha,falpha: ',i4,5f16.12) if (alpha(3).gt.1.0d0+1.0d2*tolr) then alpha(3)=1.0d0 else if (alpha(3).lt.0.0d0-1.0d2*tolr) then alpha(3)=0.0d0 end if c check that new guess is not outside previous bounds c if (lfaneg .and. lfapos .and. (alpha(3).gt.MAX(fapos,faneg) c & .or. alpha(3).lt.MIN(fapos,faneg))) then c if so, use reguli-falsi c alpha(3)=apos-fapos*(apos-aneg)/(fapos-faneg) c end if c discard oldest iteration alpha(1)=alpha(2) alpha(2)=alpha(3) falpha(1)=falpha(2) end if 600 continue c inner iteration loop has not converged, issue warning and proceed ierr=-218 write (herr,1218) hnull 1218 format ('[TPFLSH warning 218] inner iteration loop for liquid ', & 'fraction in 2-phase state did not converge: ',a1) call ERRMSG (ierr,herr) c c end of inner iteration loop for liquid fraction c 640 continue ysum=0.0d0 c compute next guess for vapor composition, using x(i) from inner loop do 740 i=1,nc y(i)=x(i)*xk(i) ysum=ysum+y(i) 740 continue c normalize vapor compositions do 760 i=1,nc y(i)=y(i)/ysum 760 continue c write (*,1076) iy,t,Dl,Dv,ysum,(y(i),i=1,3) c1076 format (1x,' TPFL2:',i4,f7.2,2e14.6,10x,f10.6,3f10.6) if (ABS(1.0d0-ysum).lt.tolr) then c outer loop has converged (also check for warning from inner loop) if (ierr.ne.0) then ierr=218 write (herr,2218) hnull 2218 format ('[TPFLSH error 218] inner iteration loop for liquid ', & 'fraction in 2-phase state did not converge: ',a1) call ERRMSG (ierr,herr) end if goto 840 end if 800 continue c outer iteration loop has not converged, issue error and return ierr=219 write (herr,1219) hnull 1219 format ('[TPFLSH error 219] outer iteration loop for ', & 'composition in 2-phase state did not converge: ',a1) call ERRMSG (ierr,herr) c c end of outer iteration loop for vapor composition c 840 continue q=1.0d0-alpha(j) c RETURN end !subroutine TPFL2 c c ====================================================================== c subroutine TDFLSH (t,D,z,p,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr) c c flash calculation given temperature, bulk density, & bulk composition c c This routine accepts both single-phase and two-phase states as input, c if the phase is known, the subroutine THERM is (much) faster. c c inputs: c t--temperature [K] c D--overall (bulk) molar density [mol/L] c z--overall (bulk) composition [array of mol frac] c c outputs: c p--pressure [kPa] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c if only one phase is present, Dl = Dv = D c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c if only one phase is present, xl = xv = x c q--vapor quality on a MOLAR basis [moles vapor/total moles] c q < 0 indicates subcooled (compressed) liquid c q = 0 indicates saturated liquid c q = 1 indicates saturated vapor c q > 1 indicates superheated vapor c q = 998 superheated vapor, but quality not defined (t > Tc) c q = 999 indicates supercritical state (t > Tc) and (p > Pc) c e--overall (bulk) internal energy [J/mol] c h--overall (bulk) enthalpy [J/mol] c s--overall (bulk) entropy [J/mol-K] c Cv--isochoric (constant V) heat capacity [J/mol-K] c Cp--isobaric (constant p) heat capacity [J/mol-K] c w--speed of sound [m/s] c Cp, w are not defined for 2-phase states c in such cases, a flag = -9.99998d6 is returned c ierr--error flag: 0 = successful c 1 = Tin < Tmin c 2 = Din > Dmax or Din < 0 c 3 = T and D out of range c 8 = x out of range (component and/or sum < 0 c or > 1) c 9 = x and T out of range c 10 = x and D out of range c 11 = x and T and D out of range c 220 = CRITP did not converge c 221 = SATT did not converge at bubble point c 222 = SATT did not converge at dew point c 223 = SATT (bubble pt) did not converge for 2-ph c 224 = SATT (dew pt) did not converge for 2-phase c 225 = TPFL2 did not converge c 226 = 2-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-20-95 MM, original version c 09-25-95 MM, rearrange argument list (outputs in order rho, x, q) c 10-11-95 MM, RETURN if any error detected c 11-08-95 MM, patch to allow testing of GUI with mixtures c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 12-13-95 MM, set undefined output compositions to zero c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 01-21-97 MM, return error message for mixtures (not yet implemented) c 04-09-97 MM, add mixture calculations; restructure cases c add limits checks c 04-22-97 MM, use h, rather than s, to compute q (h converges at p = 0) c 05-15-97 MM, use V to compute q for comp. liq. (T,h can be double-valued), c add q = 998 case for t > Tc, but p < Pc c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c 10-16-97 MM, fix bug--out of range errors not caught properly c 12-01-97 MM, error code not passed correctly if inputs out of range c 02-11-98 MM, check that computed pressure is within bounds c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: TDFLSH c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1,herr2,herrl dimension x(ncmax),y(ncmax),z(ncmax) dimension xbub(ncmax),xdew(ncmax),ybub(ncmax),ydew(ncmax) common /Gcnst/ R common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc c ierr=0 herr=hnull c c set output liquid and vapor compositions to input values c zero output values for undefined components do 100 i=1,nc x(i)=z(i) y(i)=z(i) 100 continue if (nc.lt.ncmax) then do 120 i=nc+1,ncmax x(i)=0.0d0 y(i)=0.0d0 120 continue end if c c check that input conditions are within limits c pdum=0.0d0 call LIMITX ('EOS',t,D,pdum,z,tmin,tmax,rhomax,pmax,ierrl,herrl) if (ierrl.lt.0) then c one or inputs are outside limits--if just a warning proceed w/ calc write (herr,1006) ierrl,herrl(1:234),hnull 1006 format ('[TDFLSH warning',i3,'] ',a234,a1) call ERRMSG (ierrl,herr) c if error (as opposed to warning) set output density to zero and return else if (ierrl.gt.0) then write (herr,1007) ierrl,herrl(1:236),hnull 1007 format ('[TDFLSH error',i3,'] ',a236,a1) call ERRMSG (ierrl,herr) q=999.d0 !quality undefined ierr=ierrl RETURN end if c call CRITP (z,tc,pc,rhoc,ierr,herr2) if (ierr.ne.0) then ierr=220 write (herr,1008) ierr,herr2(1:235),hnull 1008 format ('[TDFLSH error 220] ',a235,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c if (t.ge.tc .or. D.lt.1.0d-10) then c super-critical state or rho = 0 (x = y = z as set above) call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt) Dl=D Dv=D if (p.le.pc) then q=998.d0 !superheated, but cannot define quality else q=999.d0 !quality undefined end if else c sub-critical state--call saturation routine to determine phase call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=222 write (herr,1222) herr1(1:193),hnull 1222 format ('[TDFLSH error 222] dew point calculation ', & 'did not converge: ',a193,a1) call ERRMSG (ierr,herr) D=0.0d0 Dl=0.0d0 Dv=0.0d0 q=999.d0 !quality undefined RETURN end if if (nc.eq.1) then c if pure, bubble point density provided by call to SATT at dew point Dlbub=Dldew else c if mixture, call saturation routine again at bubble point call SATT (t,z,1,pbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=221 write (herr,1221) herr1(1:190),hnull 1221 format ('[TPFLSH error 221] bubble point calculation ', & 'did not converge: ',a190,a1) call ERRMSG (ierr,herr) D=0.0d0 Dl=0.0d0 Dv=0.0d0 q=999.d0 !quality undefined RETURN end if end if if (D.le.Dvdew .or. D.ge.Dlbub) then c single-phase (liq or vapor) (pure or mixture) (x = y = z as set above) call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt) Dl=D Dv=D call ENTHAL (t,Dvdew,z,hvdew) call ENTHAL (t,Dlbub,z,hlbub) c q=(h-hlbub)/(hvdew-hlbub) c compute quality based on volumes (possible for T,h to be c double-valued in compressed liquid) q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub) else if (nc.eq.1) then c special case: pure-fluid two-phase (x = y = z as set above) p=pdew q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub) Dl=Dlbub Dv=Dvdew else c general case: mixture 2-phase ksat=1 !bubble and dew point data provided to TDFL2 call TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew, & p,Dl,Dv,x,y,q,ierr,herr) if (ierr.ne.0) then c two-phase iteration did not converge--error message written by TDFL2 q=999.d0 !quality undefined RETURN end if c write (*,1225) p,q,(x(i),i=1,3),(y(i),i=1,3) c1225 format (1x,' TDFLSH--TDFL2 return p,q,x,y = ',8f12.7) end if c compute remaining properties for 2-phase states call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt) call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt) if (nc.ge.2) p=ptherm e=q*ev+(1.0d0-q)*el h=q*hv+(1.0d0-q)*hl s=q*sv+(1.0d0-q)*sl w=xnotd !Cp,w not defined for 2-phase states cp=xnotd cv=xnotd end if end if c c if limits check resulted in warning (as opposed to error) return c that message; do this again in case intermediate iteration did not c converge (thereby overwriting any warning message from LIMITX) c if (ierrl.lt.0) then ierr=ierrl write (herr,2006) ierrl,herrl(1:234),hnull 2006 format ('[TDFLSH warning',i3,'] ',a234,a1) call ERRMSG (ierr,herr) end if c if (p.gt.pmax) then c check that computed (output) pressure is within bounds of EOS ierr=4 write (herr,2008) ierr,hnull 2008 format ('[TDFLSH error',i3,']; input T,rho correspond to ', & 'a pressure above the limit of the EOS.',a1) call ERRMSG (ierr,herr) end if c RETURN end !subroutine TDFLSH c c ====================================================================== c subroutine TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew, & p,Dl,Dv,x,y,q,ierr,herr) c c flash calculation given temperature and bulk density and composition c c This routine accepts only two-phase states as input; it is intended c primarily for use by the general temperature-density flash routine c TDFLSH. It may be called independently if the state is known to be c two-phase. But beware--this routine does not check limits, and it c will be significantly faster than TDFLSH only if the bubble and dew c point limits can be provided (ksat = 1 option). c c This routine calls TPFL2 within a secant-method iteration for c pressure to find a solution. Initial guesses are based on liquid c density at the bubble point and vapor density at the dew point. c c inputs: c t--temperature [K] c D--overall (bulk) molar density [mol/L] c z--overall (bulk) composition [array of mol frac] c ksat--flag for bubble and dew point limits c 0 = dew and bubble point limits computed here c 1 = must provide values for following: c pbub--bubble point pressure [kPa] at (t,x=z) c pdew--dew point pressure [kPa] at (t,y=z) c Dlbub--liquid density [mol/L] at bubble point c Dvdew--vapor density [mol/L] at dew point c ybub--vapor composition [array of mol frac] at bubble point c xdew--liquid composition [array of mol frac] at dew point c c outputs: c p--pressure [kPa] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c q--vapor quality on a MOLAR basis [moles vapor/total moles] c ierr--error flag: 0 = successful c 223 = SATT did not converge at bubble point c 224 = SATT did not converge at dew point c 225 = TPFL2 did not converge c 226 = 2-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Physical & Chemical Properties Div, Boulder, CO c 04-10-97 MM, original version, adapted from PSFL1 c 04-28-97 MM, move generation of initial guesses here from TDFLSH; c use dew and bubble points as bounds on the iteration c 04-30-97 MM, if TPFL2 does not converge return intial guesses c 05-01-97 MM, add option to provide bubble, dew limits as inputs c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: TDFL2 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 dimension z(ncmax),x(ncmax),y(ncmax), & xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax), & x1(ncmax),y1(ncmax) dimension pj(3),fpj(2) common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc data itmax /20/ c c calculate bubble and dew points to establish bounds on iteration; c store results for use in possible reguli-falsi if (ksat.ne.1) then call SATT (t,z,1,pbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=223 write (herr,1223) herr1(1:164),hnull 1223 format ('[TDFLSH error 223] bubble point calculation did ', & 'not converge (2-phase initial guess): ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=224 write (herr,1224) herr1(1:167),hnull 1224 format ('[TDFLSH error 224] dew point calculation did ', & 'not converge (2-phase initial guess): ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if end if c base convergence tolerance on pbub/pdew ratio: otherwise, near- c azeotropes present impossible iteration (yet initial guesses are c very good) tolv=1.0d-5 pratio=pbub/pdew-1.0d0 if (ABS(pratio).gt.1.0d-12) then tolv=MAX(tolv/1.0d2,MIN(tolv*1.0d2,tolv/pratio)) c write (*,*) ' TDFL2--adjusted tolv: ',tolv end if pjpos=pbub fpjpos=D*(1.0d0/D-1.0d0/Dlbub) pj(1)=pdew fpj(1)=D*(1.0d0/D-1.0d0/Dvdew) pjneg=pdew fpjneg=fpj(1) c write (*,1004) pjpos,fpjpos,pjneg,fpjneg c1004 format (1x,' TDFL2--p,fp bounds: ',4f12.6) c write (*,1006) -1,0,pjpos,Dlbub,xnotd,z(1),ybub(1),fpjpos c write (*,1006) 0,1,pj(1),xnotd,Dvdew,xdew(1),z(1),fpj(1) c1006 format (' TDFL2: it,j,pj,Dl,Dv,x,y,fpj: ',2i3,5(1x,f12.6),f18.10) c c define third guess for pressure j=2 c initial guesses for quality, liquid & vapor density & composition q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 240 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 240 continue do 244 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum x1(i)=x(i) !save initial guesses for use in case y1(i)=y(i) !TPFL2 does not converge 244 continue c base initial guesses for density and pressure on above guesses for D,x call SATT (t,x,1,pliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=223 write (herr,2223) herr1(1:164),hnull 2223 format ('[TDFLSH error 223] bubble point calculation did ', & 'not converge within 2-phase iteration: ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATT (t,y,2,pvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=224 write (herr,2224) herr1(1:167),hnull 2224 format ('[TDFLSH error 224] dew point calculation did ', & 'not converge within 2-phase iteration: ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes pj(2)=(1.0d0-q)*pliq+q*pvap c write (*,1244) pj(2),q,(x(i),i=1,3),(y(i),i=1,3) c1244 format (1x,' TDFL2 initial guesses--p,q,x,y = ',8f12.7) nfail=0 pj1=pj(2) Dl1=Dl !save initial guesses for use in case Dv1=Dv !TPFL2 does not converge q1=q do 400 it=2,itmax ierr1=0 !reset possible error flag from TPFL2 herr1=hnull c find density and composition for each phase c note that the initial guesses for Dl,Dv,x,y,q are refined in TPFL2 call TPFL2 (t,pj(j),z,Dl,Dv,x,y,q,ierr1,herr1) if (ierr1.ne.0) then nfail=nfail+1 c TPFL2 has not converged c write (*,1016) pj(j),q,herr1 c1016 format (1x,' TDFL2--TPFL2 returns error for p = ',f14.6, c & ' kPa, quality = ',f10.6,': ',a255) if (nfail.ge.2) then c TPFL2 has not converged a second time--and probably never will c return the initial guesses calculated above do 324 i=1,nc x(i)=x1(i) !retrieve initial guesses y(i)=y1(i) 324 continue Dl=Dl1 Dv=Dv1 p=pj1 q=q1 if (q.lt.0.005d0 .or. q.gt.0.995d0) then c initial guesses are excellent very close to saturation ierr=0 herr=hnull c write (*,*) ' TDFL2--returning initial guesses' else ierr=225 herr='[TDFLSH error 225] TPFL2 did not converge in two-'// & 'phase iteration; returning initial guesses based on'// & ' dew and bubble points'//hnull end if RETURN end if c set next guess for pressure by moving 1 % of way towards q = 0.5 q=q+SIGN(0.01d0,0.5d0-q) c write (*,1034) q c1034 format (1x,' TDFL2--next guess for quality = ',33x,f10.6) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 340 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 340 continue do 344 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum 344 continue call SATT (t,x,1,pliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=223 write (herr,2223) herr1(1:164),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATT (t,y,2,pvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=224 write (herr,2224) herr1(1:167),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes call PRESS(t,Dv,y,pj(j)) Vj=(1.0d0-q)/Dl+q/Dv !estimate of new bulk volume fpj(j)=D*(1.0d0/D-Vj) !objective function c write (*,1036) it,j,pj(j),Dl,Dv,x(1),y(1),fpj(j) c1036 format (1x,'TDFL2: reset iteration to: ',2i4,5f12.6,f18.10) goto 400 end if Vj=(1.0d0-q)/Dl+q/Dv !bulk volume for current iteration fpj(j)=D*(1.0d0/D-Vj) !objective function c write (*,1040) it,j,pj(j),Dl,Dv,x(1),y(1),fpj(j) c1040 format (' TDFL2: it,j,pj,Dl,Dv,x,y,fpj: ',2i3,5(1x,f12.6),f18.10) if (ABS(fpj(j)).lt.tolv) then c iteration has converged p=pj(j) ierr=ierr1 herr=herr1 RETURN else if (fpj(j).lt.0.0d0) then c store "negative" guess for use in possible reguli-falsi iteration c lneg=.true. pjneg=pj(j) fpjneg=fpj(j) else c store "positive" guess for use in possible reguli-falsi iteration c lpos=.true. pjpos=pj(j) fpjpos=fpj(j) end if c define next guess for temperature by secant method if (ABS(fpj(2)-fpj(1)).gt.1.0d-12) then !check: divide by zero pj(3)=pj(2)-fpj(2)*(pj(2)-pj(1))/(fpj(2)-fpj(1)) else pj(3)=pj(2)+0.5d0*(pj(2)-pj(1)) end if if (pj(3).gt.max(pjpos,pjneg).or.pj(3).lt.min(pjpos,pjneg)) then c secant method has yielded guess further from solution, use reguli-falsi c write (*,*) 'TDFL2--using reguli-falsi; p by secant: ',pj(3) pj(3)=pjpos-fpjpos*(pjpos-pjneg)/(fpjpos-fpjneg) c write (*,*) ' next pressure by reguli-falsi: ',pj(3) end if pj(1)=pj(2) pj(2)=pj(3) fpj(1)=fpj(2) if (ABS(pj(2)-pj(1)).gt.0.25*pj(1)) then c do not permit too large a step pj(2)=pj(1)+pj(1)*SIGN(0.25d0,pj(2)-pj(1)) end if 400 continue c c iteration has not converged c ierr=226 write (herr,1226) pj(j),fpj(j),hnull 1226 format ('[TDFLSH error 226] 2-phase iteration did not converge, ', & 'P, deltaP =',2(1pe11.3),' kPa.',a1) call ERRMSG (ierr,herr) c RETURN end !subroutine TDFL2 c c ====================================================================== c subroutine PHFLSH (p,h,z,t,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr) c c flash calculation given pressure, bulk enthalpy, and bulk composition c c inputs: c p--pressure [kPa] c h--overall (bulk) enthalpy [J/mol] c z--composition [array of mol frac] c c outputs: c t--temperature [K] c D--overall (bulk) molar density [mol/L] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c if only one phase is present, Dl = Dv = D c x--composition [array of mol frac] for liquid phase c y--composition [array of mol frac] for vapor phase c if only one phase is present, x = y = z c q--vapor quality on a MOLAR basis [moles vapor/total moles] c q < 0 indicates subcooled (compressed) liquid c q = 0 indicates saturated liquid c q = 1 indicates saturated vapor c q > 1 indicates superheated vapor c q = -998 subcooled liquid, but quality not defined (p > Pc) c q = 999 indicates supercritical state (t > Tc) and (p > Pc) c e--overall (bulk) internal energy [J/mol] c s--overall (bulk) entropy [J/mol.K] c Cv--isochoric (constant V) heat capacity [J/mol.K] c Cp--isobaric (constant p) heat capacity [J/mol.K] c w--speed of sound [m/s] c Cp, w are not defined for 2-phase states c in such cases, a flag = -9.99998d6 is returned c ierr--error flag: 0 = successful c 4 = Pin < 0 c 8 = x out of range (< 0 or > 1) c 12 = x out of range and P < 0 c 240 = CRITP did not converge c 241 = SATP did not converge at bubble point c 242 = SATP did not converge at dew point c 243 = SATP (bubble pt) did not converge for 2-ph c 244 = SATP (dew pt) did not converge for 2-phase c 245 = TPFL2 did not converge c 246 = 2-phase iteration did not converge c 247 = TPRHO did not converge for single-phase c 248 = single-phase iteration did not converge c 249 = H out of range c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-20-95 MM, original version c 09-25-95 MM, rearrange argument list (outputs in order rho, x, q) c 10-11-95 MM, correct: t not returned for two-phase pure-component c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 12-13-95 MM, set undefined output compositions to zero c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 04-30-97 MM, restructure logic to parallel TDFLSH c 05-15-97 MM, add q = -998 case for p > Pc, but t < Tc c 05-21-97 MM, call ENTHAL at rho = 1d-10 to define hmax c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c 10-16-97 MM, evaluate hmin at p = 100 kPa, rather than p = pmax c 02-11-98 MM, check that computed density is within bounds c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PHFLSH c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1,herr2 dimension x(ncmax),y(ncmax),z(ncmax), & xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax) common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc c ierr=0 herr=hnull c c set output liquid and vapor compositions to input values c zero output values for undefined components do 100 i=1,nc x(i)=z(i) y(i)=z(i) 100 continue if (nc.lt.ncmax) then do 120 i=nc+1,ncmax x(i)=0.0d0 y(i)=0.0d0 120 continue end if c call CRITP (z,tc,pc,Dc,ierr,herr1) if (ierr.ne.0) then ierr=240 write (herr,1002) ierr,herr1(1:234),hnull 1002 format ('[PHFLSH error 240] ',a234,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c c check that input conditions are within limits c tdum=0.75*tc Ddum=0.0d0 call LIMITX ('EOS',tdum,Ddum,p,z,tmin,tmax,Dmax,pmax,ierr,herr1) c calculate approx enthalpy limits: c lower limit at (tmin,pmax); upper limit at (1.5tmax,rho = 0) c !MM, 10-16-97, because of curvature of isotherms, h at pmax is not c the minimum, use a lower pressure instead p0=100.0d0 !kPa call TPRHO (tmin,p0,z,1,0,D0,ierr2,herr2) call ENTHAL (tmin,D0,z,hmin) tmax=1.5d0*tmax rho0=1.0d-10 call ENTHAL (tmax,rho0,z,hmax) c write (*,1004) tmin,tmax,hmin,hmax c1004 format (1x,' PHFLSH--limits on t,h: ',2f8.2,2f12.2) c if inputs are outside limits set outputs equal to critical c point values and return if (ierr.lt.-1) then !ignore ierr = -1 (t out of range) c one or inputs are outside limits--if just a warning proceed w/ calc write (herr,1006) ierr,herr1(1:234),hnull 1006 format ('[PHFLSH warning',i4,'] ',a234,a1) call ERRMSG (ierr,herr) else if (ierr.gt.1) then c ignore ierr=1 out of LIMITX (t out of range, but t not an input) write (herr,1007) ierr,herr1(1:235),hnull 1007 format ('[PHFLSH error',i4,'] ',a235,a1) else if (h.lt.hmin .or. h.gt.hmax) then ierr=249 write (herr,1008) ierr,h,hmin,hmax,hnull 1008 format ('[PHFLSH error',i4,'] Input enthalpy is ', & 'outside limits; H =',1pe11.3,' J/mol, Hmin,Hmax =', & 2(1pe11.3),' J/mol.',a1) end if if (ierr.gt.1) then call ERRMSG (ierr,herr) t=tc D=Dc Dl=Dc Dv=Dc call THERM (t,D,z,ptherm,e,htherm,s,cv,cp,w,hjt) q=999.d0 !quality undefined RETURN end if end if c if (p.ge.pc .or. p.lt.1.0d-12) then c super-critical state or special case for p = 0 (x = y = z as set above) kph=2 !use vapor iteration t=tc !initial guesses if (p.ge.pc) then D=Dc else D=0.0d0 end if call PHFL1 (p,h,z,kph,t,D,ierr,herr) if (ierr.ne.0) then c single-phase iteration did not converge--error message written by PHFL1 q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,ptherm,e,htherm,s,cv,cp,w,hjt) Dl=D Dv=D if (t.le.tc) then q=-998.d0 !subcooled, but cannot define quality else q=999.d0 !supercritical--quality undefined end if else c sub-critical state--call saturation routine to determine phase call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=242 write (herr,1242) herr1(1:192),hnull 1242 format ('[PHFLSH error 242] dew point calculation ', & 'did not converge: ',a192,a1) RETURN end if if (nc.eq.1) then c if pure, bubble point density provided by call to SATT at dew point Dlbub=Dldew tbub=tdew else c if mixture, call saturation routine again at bubble point call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=241 write (herr,1241) herr1(1:188),hnull 1241 format ('[TPFLSH error 241] bubble point calculation ', & 'did not converge: ',a188,a1) call ERRMSG (ierr,herr) RETURN end if end if c calculate enthalpy at dew and bubble points call ENTHAL (tdew,Dvdew,z,hvdew) call ENTHAL (tbub,Dlbub,z,hlbub) q=(h-hlbub)/(hvdew-hlbub) if (h.le.hlbub .or. h.ge.hvdew) then c single-phase (liq or vapor) (pure or mixture) (x = y = z as set above) c set initial guesses, depending on phase if (h.le.hlbub) then kph=1 !liquid t=tbub !initial guesses D=Dlbub else kph=2 !vapor t=tdew !initial guesses D=Dvdew end if call PHFL1 (p,h,z,kph,t,D,ierr,herr) if (ierr.ne.0) then c single-phase iteration did not converge--error message written by PHFL1 call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt) Dl=D Dv=D else if (nc.eq.1) then c special case: pure-fluid two-phase (x = y = z as set above) t=tbub Dl=Dlbub Dv=Dvdew else c general case: mixture 2-phase ksat=1 !bubble and dew point data provided to PHFL2 call PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, & t,Dl,Dv,x,y,q,ierr,herr) if (ierr.ne.0) then c two-phase iteration did not converge--error message written by PHFL2 call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c write (*,1225) t,q,(x(i),i=1,3),(y(i),i=1,3) c1225 format (1x,' PHFLSH--PHFL2 return t,q,x,y = ',8f12.7) end if c compute remaining properties for 2-phase states call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt) call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt) alpha=1.0d0-q D=1.0d0/(alpha/Dl+q/Dv) e=alpha*el+q*ev s=alpha*sl+q*sv w=xnotd !Cp,w not defined for 2-phase states cp=xnotd cv=xnotd end if end if c if (D.gt.Dmax) then c check that computed (output) density is within bounds of EOS ierr=2 write (herr,2008) ierr,hnull 2008 format ('[PHFLSH error',i3,']; input P,H correspond to ', & 'a density above the limit of the EOS.',a1) call ERRMSG (ierr,herr) end if c RETURN end !subroutine PHFLSH c c ====================================================================== c subroutine PHFL1 (p,h,x,kph,t,D,ierr,herr) c c flash calculation given pressure, enthalpy, and composition c c N.B. This routine accepts only single-phase inputs, it is intended c primarily for use with the more general flash routine PHFLSH. c c inputs: c p--pressure [kPa] c h--enthalpy [J/mol] c x--composition [array of mol frac] c kph--phase flag: 1 = liquid c 2 = vapor c t--initial guess for temperature [K] c D--initial guess for molar density [mol/L] c c outputs: c t--temperature [K] c D--molar density [mol/L] c ierr--error flag: 0 = successful c 247 = TPRHO did not converge c 248 = single-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-22-95 MM, original version c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PHFL1 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 logical lpos,lneg dimension x(ncmax) dimension z(3),fz(2) !z is current guess for temperature common /NCOMP/ nc common /CHAR/ htab,hnull data itmax /12/ !temp--move to common?? data tolh /1.0d-2/ c ierr=0 herr=hnull c j=1 lpos=.false. !initialize flags for reguli-false lneg=.false. zneg=0.0d0 zpos=0.0d0 fzneg=0.0d0 fzpos=0.0d0 c z(1)=t do 200 it=1,itmax t=z(j) call TPRHO (t,p,x,kph,1,D,ierr,herr1) !find density if (ierr.ne.0) then ierr=247 write (herr,1247) herr1(1:172),hnull 1247 format ('[PHFLSH error 247] density iteration for ', & 'single-phase state did not converge: ',a172,a1) call ERRMSG (ierr,herr) RETURN end if call ENTHAL (t,D,x,ht) !find enthalpy for current t fz(j)=h-ht c write (*,1020) it,j,z(j),D,p,h,ht,fz(j) c1020 format (1x,'PHFL1: it,j,t,D,p,h,ht,ft: ',2i4,f10.5,f12.7,4f15.6) if (ABS(fz(j)).lt.tolh) then t=z(j) ierr=0 herr=hnull RETURN else if (fz(j).lt.0.0d0) then c store "negative" guess for use in possible reguli-falsi iteration lneg=.true. zneg=z(j) fzneg=fz(j) else c store "positive" guess for use in possible reguli-falsi iteration lpos=.true. zpos=z(j) fzpos=fz(j) end if if (j.eq.1) then c define second guess for temperature call THERM (z(1),D,x,ptherm,e,ht,s,cv,cp,w,hjt) delz=fz(1)/cp if (ABS(delz).lt.0.0005d0*z(1)) then delz=sign(0.0005d0*z(1),delz) end if z(2)=z(1)+delz j=2 else c general case--define next guess for temperature by secant method if (ABS(fz(2)-fz(1)).gt.1.0d-12) then !check for divide by zero z(3)=z(2)-fz(2)*(z(2)-z(1))/(fz(2)-fz(1)) else z(3)=z(2)+0.5d0*(z(2)-z(1)) end if if (lneg .and. lpos .and. & (z(3).gt.MAX(zpos,zneg) .or. z(3).lt.MIN(zpos,zneg))) then c secant method has yielded guess further from solution, use reguli-falsi c write (*,*) 'PHFL1--using reguli-falsi; t by secant: ',z(3) z(3)=zpos-fzpos*(zpos-zneg)/(fzpos-fzneg) c write (*,*) ' next temperature by reguli-falsi: ',z(3) end if z(1)=z(2) z(2)=z(3) fz(1)=fz(2) end if if (ABS(z(2)-z(1)).gt.0.25*z(1)) then c do not permit too large a step in temperature z(2)=z(1)+z(1)*sign(0.25d0,z(2)-z(1)) end if 200 continue c c iteration has not converged c ierr=248 write (herr,1248) z(j),fz(j),hnull 1248 format ('[PHFLSH error 248] single-phase iteration did not ', & 'converge, T, deltaT =',2(1pe11.3),' K.',a1) call ERRMSG (ierr,herr) c RETURN end !subroutine PHFL1 c c ====================================================================== c subroutine PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, & t,Dl,Dv,x,y,q,ierr,herr) c c flash calculation given pressure and bulk enthalpy and composition c c This routine accepts only two-phase states as input; it is intended c primarily for use by the general pressure-enthalpy flash routine c PHFLSH. It may be called independently if the state is known to be c two-phase. But beware--this routine does not check limits, and it c will be significantly faster than PHFLSH only if the bubble and dew c point limits can be provided (ksat = 1 option). c c This routine calls TPFL2 within a secant-method iteration for c pressure to find a solution. Initial guesses are based on liquid c density at the bubble point and vapor density at the dew point. c c inputs: c p--pressure [kPa] c h--overall (bulk) molar enthalpy [J/mol] c z--overall (bulk) composition [array of mol frac] c ksat--flag for bubble and dew point limits c 0 = dew and bubble point limits computed here c 1 = must provide values for following: c tbub--bubble point temperature [K] at (p,x=z) c tdew--dew point temperature [K] at (p,y=z) c Dlbub--liquid density [mol/L] at bubble point c Dvdew--vapor density [mol/L] at dew point c ybub--vapor composition [array of mol frac] at bubble point c xdew--liquid composition [array of mol frac] at dew point c c outputs: c t--temperature [K] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c q--vapor quality on a MOLAR basis [moles vapor/total moles] c ierr--error flag: 0 = successful c 243 = SATP did not converge at bubble point c 244 = SATP did not converge at dew point c 245 = TPFL2 did not converge c 246 = 2-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Physical & Chemical Properties Div, Boulder, CO c 04-30-97 MM, original version, based on TDFL2 c 05-16-97 MM, skip calls to SATP if ksat = 1 c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PHFL2 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 dimension z(ncmax),x(ncmax),y(ncmax), & xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax), & x1(ncmax),y1(ncmax) dimension tj(3),ftj(2) common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc data itmax /20/ c c bubble and dew points to establish bounds on iteration; c store results for use in possible reguli-falsi if (ksat.ne.1) then call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=243 write (herr,1243) herr1(1:164),hnull 1243 format ('[PHFLSH error 243] bubble point calculation did ', & 'not converge (2-phase initial guess): ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=244 write (herr,1244) herr1(1:167),hnull 1244 format ('[PHFLSH error 244] dew point calculation did ', & 'not converge (2-phase initial guess): ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if end if call ENTHAL (tbub,Dlbub,z,hlbub) call ENTHAL (tdew,Dvdew,z,hvdew) c base convergence tolerance on tdew/tbub ratio: otherwise, near- c azeotropes present impossible iteration (yet initial guesses are c very good) tolh=1.0d-3 tratio=tdew/tbub-1.0d0 if (ABS(tratio).lt.1.0d-8) tratio=1.0d-8 tolh=MAX(tolh/1.0d2,MIN(tolh*1.0d2,0.02d0*tolh/tratio)) c write (*,*) ' PHFL2--adjusted tolh: ',tolh tjpos=tbub ftjpos=h-hlbub tj(1)=tdew ftj(1)=h-hvdew tjneg=tj(1) ftjneg=ftj(1) c write (*,1004) tjpos,ftjpos,tjneg,ftjneg c1004 format (1x,' PHFL2--t,ft bounds: ',2(f12.6,f14.4,3x)) c write (*,1006) -1,0,tjpos,Dlbub,xnotd,z(1),ybub(1),ftjpos c write (*,1006) 0,1,tj(1),xnotd,Dvdew,xdew(1),z(1),ftj(1) c1006 format (' PHFL2: it,j,tj,Dl,Dv,x,y,ftj: ',2i3,5(1x,f12.6),f18.10) c c define third guess for temperature j=2 c initial guesses for quality, liquid & vapor density & composition q=(h-hlbub)/(hvdew-hlbub) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 240 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 240 continue do 244 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum x1(i)=x(i) !save initial guesses for use in case y1(i)=y(i) !TPFL2 does not converge 244 continue c base initial guesses for density and pressure on above guesses for D,x call SATP (p,x,1,tliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=243 write (herr,2243) herr1(1:164),hnull 2243 format ('[PHFLSH error 243] bubble point calculation did ', & 'not converge within 2-phase iteration: ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,y,2,tvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=244 write (herr,2244) herr1(1:167),hnull 2244 format ('[PHFLSH error 244] dew point calculation did ', & 'not converge within 2-phase iteration: ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes tj(2)=(1.0d0-q)*tliq+q*tvap c write (*,1009) tj(2),q,(x(i),i=1,3),(y(i),i=1,3) c1009 format (1x,' PHFL2 initial guesses--t,q,x,y = ',8f12.7) nfail=0 tj1=tj(2) Dl1=Dl !save initial guesses for use in case Dv1=Dv !TPFL2 does not converge q1=q do 400 it=2,itmax ierr1=0 !reset possible error flag from TPFL2 herr1=hnull c find density and composition for each phase c note that the initial guesses for Dl,Dv,x,y,q are refined in TPFL2 call TPFL2 (tj(j),p,z,Dl,Dv,x,y,q,ierr1,herr1) if (ierr1.ne.0) then c TPFL2 has not converged nfail=nfail+1 c write (*,1016) tj(j),q,herr1 c1016 format (1x,' PHFL2--TPFL2 returns error for t = ',f14.6, c & ' K, quality = ',f10.6,': ',a255) if (nfail.ge.2) then c TPFL2 has not converged a second time--and probably never will c return the initial guesses calculated above do 324 i=1,nc x(i)=x1(i) !retrieve initial guesses y(i)=y1(i) 324 continue Dl=Dl1 Dv=Dv1 t=tj1 q=q1 if (q.lt.0.005d0 .or. q.gt.0.995d0) then c initial guesses are excellent very close to saturation ierr=0 herr=hnull c write (*,*) ' PHFL2--returning initial guesses' else ierr=245 herr='[PHFLSH error 245] TPFL2 did not converge in two-'// & 'phase iteration; returning initial guesses based on'// & ' dew and bubble points'//hnull end if RETURN end if c set next guess for temperature by moving 1 % of way towards q = 0.5 q=q+SIGN(0.01d0,0.5d0-q) c write (*,1034) q d1034 format (1x,' PHFL2--next guess for quality = ',31x,f10.6) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 340 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 340 continue do 344 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum 344 continue call SATP (p,x,1,tliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=243 write (herr,2243) herr1(1:164),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,y,2,tvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=244 write (herr,2244) herr1(1:167),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes tj(j)=(1.0d0-q)*tliq+q*tvap call TPFL2 (tj(j),p,z,Dl,Dv,x,y,q,ierr,herr) call ENTHAL (tj(j),Dl,x,hl) call ENTHAL (tj(j),Dv,y,hv) ftj(j)=h-(1.0d0-q)*hl-q*hv c write (*,1036) it,j,tj(j),Dl,Dv,x(1),y(1),ftj(j) c1036 format (1x,'PHFL2: reset iteration to: ',2i4,5f12.6,f18.10) goto 400 end if call ENTHAL (tj(j),Dl,x,hl) call ENTHAL (tj(j),Dv,y,hv) ftj(j)=h-(1.0d0-q)*hl-q*hv !objective function c write (*,1040) it,j,tj(j),Dl,Dv,x(1),y(1),ftj(j) c1040 format (' PHFL2: it,j,tj,Dl,Dv,x,y,ftj: ',2i3,5(1x,f12.6),f18.10) if (ABS(ftj(j)).lt.tolh) then c iteration has converged t=tj(j) ierr=ierr1 herr=herr1 RETURN else if (ftj(j).lt.0.0d0) then c store "negative" guess for use in possible reguli-falsi iteration c lneg=.true. tjneg=tj(j) ftjneg=ftj(j) else c store "positive" guess for use in possible reguli-falsi iteration c lpos=.true. tjpos=tj(j) ftjpos=ftj(j) end if c define next guess for temperature by secant method if (ABS(ftj(2)-ftj(1)).gt.1.0d-12) then !check: divide by zero tj(3)=tj(2)-ftj(2)*(tj(2)-tj(1))/(ftj(2)-ftj(1)) else tj(3)=tj(2)+0.5d0*(tj(2)-tj(1)) end if if (tj(3).gt.max(tjpos,tjneg).or.tj(3).lt.min(tjpos,tjneg)) then c secant method has yielded guess further from solution, use reguli-falsi c write (*,*) 'PHFL2--using reguli-falsi; t by secant: ',tj(3) tj(3)=tjpos-ftjpos*(tjpos-tjneg)/(ftjpos-ftjneg) c write (*,*) ' next temperature by reguli-falsi: ',tj(3) end if tj(1)=tj(2) tj(2)=tj(3) ftj(1)=ftj(2) if (ABS(tj(2)-tj(1)).gt.0.05*tj(1)) then c do not permit too large a step tj(2)=tj(1)+tj(1)*SIGN(0.05d0,tj(2)-tj(1)) end if 400 continue c c iteration has not converged c t=tj(j) ierr=246 write (herr,1246) tj(j),ftj(j),hnull 1246 format ('[PHFLSH error 246] 2-phase iteration did not converge, ', & 'T, deltaT =',2(1pe11.3),' K.',a1) call ERRMSG (ierr,herr) c RETURN end !subroutine PHFL2 c c ====================================================================== c subroutine PSFLSH (p,s,z,t,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr) c c flash calculation given pressure, bulk entropy, and bulk composition c c inputs: c p--pressure [kPa] c s--overall (bulk) entropy [J/mol.K] c z--composition array (mol frac) c c outputs: c t--temperature [K] c D--overall (bulk) molar density [mol/L] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c if only one phase is present, Dl = Dv = D c x--composition [array of mol frac] for liquid phase c y--composition [array of mol frac] for vapor phase c if only one phase is present, xl = xv = x c q--vapor quality on a MOLAR basis [moles vapor/total moles] c q < 0 indicates subcooled (compressed) liquid c q = 0 indicates saturated liquid c q = 1 indicates saturated vapor c q > 1 indicates superheated vapor c q = -998 subcooled liquid, but quality not defined (p > Pc) c q = 999 indicates supercritical state (t > Tc) and (p > Pc) c e--overall (bulk) internal energy [J/mol] c h--overall (bulk) enthalpy [J/mol] c Cv--isochoric (constant V) heat capacity [J/mol.K] c Cp--isobaric (constant p) heat capacity [J/mol.K] c w--speed of sound [m/s] c Cp, w are not defined for 2-phase states c in such cases, a flag = -9.99998d6 is returned c ierr--error flag: 0 = successful c 4 = Pin < 0 c 8 = x out of range (< 0 or > 1) c 12 = x out of range and P < 0 c 260 = CRITP did not converge c 261 = SATP did not converge at bubble point c 262 = SATP did not converge at dew point c 263 = SATP (bubble pt) did not converge for 2-ph c 264 = SATP (dew pt) did not converge for 2-phase c 265 = TPFL2 did not converge c 266 = 2-phase iteration did not converge c 267 = TPRHO did not converge for single-phase c 268 = single-phase iteration did not converge c 269 = S out of range c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-22-95 MM, original version c 09-25-95 MM, rearrange argument list (outputs in order rho, x, q) c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 12-13-95 MM, set undefined output compositions to zero c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 05-01-97 MM, restructure logic to parallel TDFLSH, PHFLSH c 05-15-97 MM, add q = -998 case for p > Pc, but t < Tc c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c 10-16-97 MM, evaluate smin at p = 100 kPa, rather than p = pmax c 02-11-98 MM, check that computed density is within bounds c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PSFLSH c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1,herr2 dimension x(ncmax),y(ncmax),z(ncmax), & xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax) common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc c ierr=0 herr=hnull c c set output liquid and vapor compositions to input values c zero output values for undefined components do 100 i=1,nc x(i)=z(i) y(i)=z(i) 100 continue if (nc.lt.ncmax) then do 120 i=nc+1,ncmax x(i)=0.0d0 y(i)=0.0d0 120 continue end if c call CRITP (z,tc,pc,Dc,ierr,herr1) if (ierr.ne.0) then ierr=260 write (herr,1002) ierr,herr1(1:234),hnull 1002 format ('[PSFLSH error 260] ',a234,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c c check that input conditions are within limits c tdum=0.75*tc Ddum=0.0d0 call LIMITX ('EOS',tdum,Ddum,p,z,tmin,tmax,Dmax,pmax,ierr,herr1) c calculate approx entropy limits: c lower limit at (tmin,pmax); upper limit at (1.5tmax,rho = 0) c !MM, 10-16-97, because of curvature of isotherms, h at pmax is not c the minimum, use a lower pressure instead p0=100.0d0 !kPa call TPRHO (tmin,p0,z,1,0,D0,ierr2,herr2) call ENTRO (tmin,D0,z,smin) tmax=1.5d0*tmax rho0=1.0d-6 call ENTRO (tmax,rho0,z,smax) c write (*,1004) tmin,tmax,smin,smax c1004 format (1x,' PSFLSH--limits on t,s: ',2f8.2,2f12.4) c if inputs are outside limits set outputs equal to critical c point values and return if (ierr.lt.-1) then !ignore ierr = -1 (t out of range) c one or inputs are outside limits--if just a warning proceed w/ calc write (herr,1006) ierr,herr1(1:234),hnull 1006 format ('[PSFLSH warning',i4,'] ',a234,a1) call ERRMSG (ierr,herr) else if (ierr.gt.1) then c ignore ierr=1 out of LIMITX (t out of range, but t not an input) write (herr,1007) ierr,herr1(1:235),hnull 1007 format ('[PSFLSH error',i4,'] ',a235,a1) else if (s.lt.smin .or. s.gt.smax) then ierr=269 write (herr,1008) ierr,s,smin,smax,hnull 1008 format ('[PHFLSH error',i4,'] Input entropy is outside ', & 'limits; S =',1pe11.3,' J/mol-K, Smin,Smax =', & 2(1pe11.3),' J/mol.',a1) end if if (ierr.gt.1) then call ERRMSG (ierr,herr) t=tc D=Dc Dl=Dc Dv=Dc call THERM (t,D,z,ptherm,e,h,stherm,cv,cp,w,hjt) q=999.d0 !quality undefined RETURN end if end if c if (p.ge.pc .or. p.lt.1.0d-12) then c super-critical state or special case for p = 0 (x = y = z as set above) kph=2 !use vapor iteration t=tc !initial guesses if (p.ge.pc) then D=Dc else D=0.0d0 end if call PSFL1 (p,s,z,kph,t,D,ierr,herr) if (ierr.ne.0) then c single-phase iteration did not converge--error message written by PSFL1 q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,ptherm,e,h,stherm,cv,cp,w,hjt) Dl=D Dv=D if (t.le.tc) then q=-998.d0 !subcooled, but cannot define quality else q=999.d0 !supercritical--quality undefined end if else c sub-critical state--call saturation routine to determine phase call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=262 write (herr,1262) herr1(1:192),hnull 1262 format ('[PSFLSH error 262] dew point calculation ', & 'did not converge: ',a192,a1) RETURN end if if (nc.eq.1) then c if pure, bubble point density provided by call to SATT at dew point Dlbub=Dldew tbub=tdew else c if mixture, call saturation routine again at bubble point call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=241 write (herr,1261) herr1(1:188),hnull 1261 format ('[TSFLSH error 261] bubble point calculation ', & 'did not converge: ',a188,a1) call ERRMSG (ierr,herr) RETURN end if end if c calculate entropy at dew and bubble points call ENTRO (tdew,Dvdew,z,svdew) call ENTRO (tbub,Dlbub,z,slbub) q=(s-slbub)/(svdew-slbub) if (s.le.slbub .or. s.ge.svdew) then c single-phase (liq or vapor) (pure or mixture) (x = y = z as set above) c set initial guesses, depending on phase if (s.le.slbub) then kph=1 !liquid t=tbub !initial guesses D=Dlbub else kph=2 !vapor t=tdew !initial guesses D=Dvdew end if call PSFL1 (p,s,z,kph,t,D,ierr,herr) if (ierr.ne.0) then c single-phase iteration did not converge--error message written by PSFL1 call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt) Dl=D Dv=D else if (nc.eq.1) then c special case: pure-fluid two-phase (x = y = z as set above) t=tbub Dl=Dlbub Dv=Dvdew else c general case: mixture 2-phase ksat=1 !bubble and dew point data provided to PHFL2 call PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, & t,Dl,Dv,x,y,q,ierr,herr) if (ierr.ne.0) then c two-phase iteration did not converge--error message written by PSFL2 call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if c write (*,1225) t,q,(x(i),i=1,3),(y(i),i=1,3) c1225 format (1x,' PSFLSH--PSFL2 return t,q,x,y = ',8f12.7) end if c compute remaining properties for 2-phase states call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt) call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt) alpha=1.0d0-q D=1.0d0/(alpha/Dl+q/Dv) e=alpha*el+q*ev h=alpha*hl+q*hv w=xnotd !Cp,w not defined for 2-phase states cp=xnotd cv=xnotd end if end if c if (D.gt.Dmax) then c check that computed (output) density is within bounds of EOS ierr=2 write (herr,2008) ierr,hnull 2008 format ('[PSFLSH error',i3,']; input P,S correspond to ', & 'a density above the limit of the EOS.',a1) call ERRMSG (ierr,herr) end if c RETURN end !subroutine PSFLSH c c ====================================================================== c subroutine PSFL1 (p,s,x,kph,t,D,ierr,herr) c c flash calculation given pressure, entropy, and composition c c N.B. This routine accepts only single-phase inputs, it is intended c primarily for use with the more general flash routine PSFLSH. c c inputs: c p--pressure [kPa] c s--entropy [J/mol.K] c x--composition [array of mol frac] c kph--phase flag: 1 = liquid c 2 = vapor c t--initial guess for temperature [K] c D--initial guess for molar density [mol/L] c c outputs: c t--temperature [K] c D--molar density [mol/L] c ierr--error flag: 0 = successful c 267 = TPRHO did not converge c 268 = single-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado c 09-22-95 MM, original version c 10-11-95 MM, RETURN if any error detected c 11-29-95 MM, variable lower limit on coefficient/constant arrays c to accomodate ECS reference fluid c 02-27-96 MM, parameter n0=-ncmax to accomodate ECS-thermo model c 01-07-97 MM, error message bug (do not concatenate herr with itself) c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PSFL1 implicit double precision (a-h,o-z) implicit integer (i-k,m,n) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 logical lpos,lneg dimension x(ncmax) dimension z(3),fz(2) !z is current guess for temperature common /NCOMP/ nc common /CHAR/ htab,hnull data itmax /12/ !temp--move to common?? data tols /1.0d-6/ c ierr=0 herr=hnull c j=1 lpos=.false. !initialize flags for reguli-false lneg=.false. zneg=0.0d0 zpos=0.0d0 fzneg=0.0d0 fzpos=0.0d0 c z(1)=t do 200 it=1,itmax t=z(j) call TPRHO (t,p,x,kph,1,D,ierr,herr1) !find density if (ierr.ne.0) then ierr=267 write (herr,1267) herr1(1:172),hnull 1267 format ('[PSFLSH error 267] density iteration for ', & 'single-phase state did not converge: ',a172,a1) call ERRMSG (ierr,herr) RETURN end if call ENTRO (t,D,x,sj) !find enthalpy for current t fz(j)=s-sj c write (*,1020) it,j,z(j),D,p,s,sj,fz(j) c1020 format (1x,'PSFL1: it,j,t,D,p,s,sj,ft: ',2i4,f10.5,2f14.5,3f15.8) if (ABS(fz(j)).lt.tols) then t=z(j) ierr=0 herr=hnull RETURN else if (fz(j).lt.0.0d0) then c store "negative" guess for use in possible reguli-falsi iteration lneg=.true. zneg=z(j) fzneg=fz(j) else c store "positive" guess for use in possible reguli-falsi iteration lpos=.true. zpos=z(j) fzpos=fz(j) end if if (j.eq.1) then c define second guess for temperature call THERM (z(1),D,x,ptherm,e,h,sj,cv,cp,w,hjt) delz=fz(1)/cv if (ABS(delz).lt.0.0005d0*z(1)) then delz=sign(0.0005d0*z(1),delz) end if z(2)=z(1)+delz j=2 else c general case--define next guess for temperature by secant method if (ABS(fz(2)-fz(1)).gt.1.0d-12) then !check for divide by zero z(3)=z(2)-fz(2)*(z(2)-z(1))/(fz(2)-fz(1)) else z(3)=z(2)+0.5d0*(z(2)-z(1)) end if if (lneg .and. lpos .and. & (z(3).gt.max(zpos,zneg) .or. z(3).lt.min(zpos,zneg))) then c secant method has yielded guess further from solution, use reguli-falsi c write (*,*) 'PSFL1--using reguli-falsi; t by secant: ',z(3) z(3)=zpos-fzpos*(zpos-zneg)/(fzpos-fzneg) c write (*,*) ' next temperature by reguli-falsi: ',z(3) end if z(1)=z(2) z(2)=z(3) fz(1)=fz(2) end if if (ABS(z(2)-z(1)).gt.0.25*z(1)) then c do not permit too large a step in temperature z(2)=z(1)+z(1)*sign(0.25d0,z(2)-z(1)) end if 200 continue c c iteration has not converged c ierr=268 write (herr,1268) z(j),fz(j),hnull 1268 format ('[PSFLSH error 268] single-phase iteration did not ', & 'converge, T, deltaT =',2(1pe11.3),' K.',a1) call ERRMSG (ierr,herr) c RETURN end !subroutine PSFL1 c c ====================================================================== c subroutine PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew, & t,Dl,Dv,x,y,q,ierr,herr) c c flash calculation given pressure and bulk entropy and composition c c This routine accepts only two-phase states as input; it is intended c primarily for use by the general pressure-entropy flash routine c PSFLSH. It may be called independently if the state is known to be c two-phase. But beware--this routine does not check limits, and it c will be significantly faster than PSFLSH only if the bubble and dew c point limits can be provided (ksat = 1 option). c c This routine calls TPFL2 within a secant-method iteration for c pressure to find a solution. Initial guesses are based on liquid c density at the bubble point and vapor density at the dew point. c c inputs: c p--pressure [kPa] c s--overall (bulk) molar entropy [J/mol-K] c z--overall (bulk) composition [array of mol frac] c ksat--flag for bubble and dew point limits c 0 = dew and bubble point limits computed here c 1 = must provide values for following: c tbub--bubble point temperature [K] at (p,x=z) c tdew--dew point temperature [K] at (p,y=z) c Dlbub--liquid density [mol/L] at bubble point c Dvdew--vapor density [mol/L] at dew point c ybub--vapor composition [array of mol frac] at bubble point c xdew--liquid composition [array of mol frac] at dew point c c outputs: c t--temperature [K] c Dl--molar density [mol/L] of the liquid phase c Dv--molar density [mol/L] of the vapor phase c x--composition of liquid phase [array of mol frac] c y--composition of vapor phase [array of mol frac] c q--vapor quality on a MOLAR basis [moles vapor/total moles] c ierr--error flag: 0 = successful c 263 = SATP did not converge at bubble point c 264 = SATP did not converge at dew point c 265 = TPFL2 did not converge c 266 = 2-phase iteration did not converge c herr--error string (character*255 variable if ierr<>0) c c written by M. McLinden, NIST Physical & Chemical Properties Div, Boulder, CO c 05-01-97 MM, original version, same logic as PHFL2 c 05-16-97 MM, skip calls to SATP if ksat = 1 c 07-15-97 MM, renumber and add detail to error messages c get flags for 'not defined' from common /FLAGS/ c 10-01-97 MM, add compiler switches to allow access by DLL c c compiler switches to allow access by DLL; for use with Digital Visual c Fortran; these should be treated as comments by all other compilers c !MS$ATTRIBUTES DLLEXPORT :: PSFL2 c implicit double precision (a-h,o-z) implicit integer (i-k,m,n) implicit logical (l) parameter (ncmax=5) !max number of components in mixture character*1 htab,hnull character*255 herr,herr1 dimension z(ncmax),x(ncmax),y(ncmax), & xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax), & x1(ncmax),y1(ncmax) dimension tj(3),ftj(2) common /NCOMP/ nc common /CHAR/ htab,hnull c flags to GUI indicating 'not applicable', '2-phase', etc. common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc data itmax /20/ c c bubble and dew points to establish bounds on iteration; c store results for use in possible reguli-falsi if (ksat.ne.1) then call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1) if (ierr.ne.0) then ierr=263 write (herr,1263) herr1(1:164),hnull 1263 format ('[PSFLSH error 263] bubble point calculation did ', & 'not converge (2-phase initial guess): ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr1) if (ierr.ne.0) then ierr=264 write (herr,1264) herr1(1:167),hnull 1264 format ('[PSFLSH error 264] dew point calculation did ', & 'not converge (2-phase initial guess): ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if end if call ENTRO (tbub,Dlbub,z,slbub) call ENTRO (tdew,Dvdew,z,svdew) c base convergence tolerance on tdew/tbub ratio: otherwise, near- c azeotropes present impossible iteration (yet initial guesses are c very good) tols=4.0d-5 tratio=tdew/tbub-1.0d0 if (ABS(tratio).lt.1.0d-8) tratio=1.0d-8 tols=MAX(tols/1.0d2,MIN(tols*1.0d2,0.02d0*tols/tratio)) c write (*,*) ' PSFL2--adjusted tols: ',tols tjpos=tbub ftjpos=s-slbub tj(1)=tdew ftj(1)=s-svdew tjneg=tj(1) ftjneg=ftj(1) c write (*,1004) tjpos,ftjpos,tjneg,ftjneg c1004 format (1x,' PSFL2--t,ft bounds: ',2(f12.6,f14.4,3x)) c write (*,1006) -1,0,tjpos,Dlbub,xnotd,z(1),ybub(1),ftjpos c write (*,1006) 0,1,tj(1),xnotd,Dvdew,xdew(1),z(1),ftj(1) c1006 format (' PSFL2: it,j,tj,Dl,Dv,x,y,ftj: ',2i3,5(1x,f12.6),f18.10) c c define third guess for temperature j=2 c initial guesses for quality, liquid & vapor density & composition q=(s-slbub)/(svdew-slbub) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 240 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 240 continue do 244 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum x1(i)=x(i) !save initial guesses for use in case y1(i)=y(i) !TPFL2 does not converge 244 continue c base initial guesses for density and pressure on above guesses for D,x call SATP (p,x,1,tliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=263 write (herr,2263) herr1(1:164),hnull 2263 format ('[PSFLSH error 263] bubble point calculation did ', & 'not converge within 2-phase iteration: ',a164,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,y,2,tvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=264 write (herr,2264) herr1(1:167),hnull 2264 format ('[PSFLSH error 264] dew point calculation did ', & 'not converge within 2-phase iteration: ',a167,a1) call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes tj(2)=(1.0d0-q)*tliq+q*tvap c write (*,1244) tj(2),q,(x(i),i=1,3),(y(i),i=1,3) c1244 format (1x,' PSFL2 initial guesses--t,q,x,y = ',8f12.7) nfail=0 tj1=tj(2) Dl1=Dl !save initial guesses for use in case Dv1=Dv !TPFL2 does not converge q1=q do 400 it=2,itmax ierr1=0 !reset possible error flag from TPFL2 herr1=hnull c find density and composition for each phase c note that the initial guesses for Dl,Dv,x,y,q are refined in TPFL2 call TPFL2 (tj(j),p,z,Dl,Dv,x,y,q,ierr1,herr1) if (ierr1.ne.0) then c TPFL2 has not converged nfail=nfail+1 c write (*,1016) tj(j),q,herr1 c1016 format (1x,' PSFL2--TPFL2 returns error for t = ',f14.6, c & ' K, quality = ',f10.6,': ',a255) if (nfail.ge.2) then c TPFL2 has not converged a second time--and probably never will c return the initial guesses calculated above do 324 i=1,nc x(i)=x1(i) !retrieve initial guesses y(i)=y1(i) 324 continue Dl=Dl1 Dv=Dv1 t=tj1 q=q1 if (q.lt.0.005d0 .or. q.gt.0.995d0) then c initial guesses are excellent very close to saturation ierr=0 herr=hnull c write (*,*) ' PSFL2--returning initial guesses' else ierr=265 herr='[PSFLSH error 265] TPFL2 did not converge in two-'// & 'phase iteration; returning initial guesses based on'// & ' dew and bubble points'//hnull end if RETURN end if c set next guess for temperature by moving 1 % of way towards q = 0.5 q=q+SIGN(0.01d0,0.5d0-q) c write (*,1034) q c1034 format (1x,' PSFL2--next guess for quality = ',31x,f10.6) xsum=0.0d0 ysum=0.0d0 !sums for normalization of compositions do 340 i=1,nc x(i)=(1.0d0-q)*z(i)+q*xdew(i) y(i)=q*z(i)+(1.0d0-q)*ybub(i) xsum=xsum+x(i) ysum=ysum+y(i) 340 continue do 344 i=1,nc x(i)=x(i)/xsum y(i)=y(i)/ysum 344 continue call SATP (p,x,1,tliq,Dlx,Dvx,xbub,ybub,ierr1,herr1) if (ierr1.ne.0) then ierr=263 write (herr,2263) herr1(1:164),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if call SATP (p,y,2,tvap,Dly,Dvy,xdew,ydew,ierr1,herr1) if (ierr1.ne.0) then ierr=264 write (herr,2264) herr1(1:167),hnull call ERRMSG (ierr,herr) q=999.d0 !quality undefined RETURN end if Dl=1.0d0/((1.0d0-q)/Dlx+q/Dly) !initial guesses for use in TPFL2 Dv=1.0d0/((1.0d0-q)/Dvx+q/Dvy) !based on volumes tj(j)=(1.0d0-q)*tliq+q*tvap call TPFL2 (tj(j),p,z,Dl,Dv,x,y,q,ierr,herr) call ENTRO (tj(j),Dl,x,sl) call ENTRO (tj(j),Dv,y,sv) ftj(j)=s-(1.0d0-q)*sl-q*sv c write (*,1036) it,j,tj(j),Dl,Dv,x(1),y(1),ftj(j) c1036 format (1x,'PSFL2: reset iteration to: ',2i4,5f12.6,f18.10) goto 400 end if call ENTRO (tj(j),Dl,x,sl) call ENTRO (tj(j),Dv,y,sv) ftj(j)=s-(1.0d0-q)*sl-q*sv !objective function c write (*,1040) it,j,tj(j),Dl,Dv,x(1),y(1),ftj(j) c1040 format (' PSFL2: it,j,tj,Dl,Dv,x,y,ftj: ',2i3,5(1x,f12.6),f18.10) if (ABS(ftj(j)).lt.tols) then c iteration has converged t=tj(j) ierr=ierr1 herr=herr1 RETURN else if (ftj(j).lt.0.0d0) then c store "negative" guess for use in possible reguli-falsi iteration c lneg=.true. tjneg=tj(j) ftjneg=ftj(j) else c store "positive" guess for use in possible reguli-falsi iteration c lpos=.true. tjpos=tj(j) ftjpos=ftj(j) end if c define next guess for temperature by secant method if (ABS(ftj(2)-ftj(1)).gt.1.0d-12) then !check: divide by zero tj(3)=tj(2)-ftj(2)*(tj(2)-tj(1))/(ftj(2)-ftj(1)) else tj(3)=tj(2)+0.5d0*(tj(2)-tj(1)) end if if (tj(3).gt.max(tjpos,tjneg).or.tj(3).lt.min(tjpos,tjneg)) then c secant method has yielded guess further from solution, use reguli-falsi c write (*,*) 'PSFL2--using reguli-falsi; t by secant: ',tj(3) tj(3)=tjpos-ftjpos*(tjpos-tjneg)/(ftjpos-ftjneg) c write (*,*) ' next temperature by reguli-falsi: ',tj(3) end if tj(1)=tj(2) tj(2)=tj(3) ftj(1)=ftj(2) if (ABS(tj(2)-tj(1)).gt.0.05*tj(1)) then c do not permit too large a step tj(2)=tj(1)+tj(1)*SIGN(0.05d0,tj(2)-tj(1)) end if 400 continue c c iteration has not converged c t=tj(j) ierr=246 write (herr,1266) tj(j),ftj(j),hnull 1266 format ('[PHFLSH error 266] 2-phase iteration did not converge, ', & 'T, deltaT =',2(1pe11.3),' K.',a1) call ERRMSG (ierr,herr) c RETURN end !subroutine PSFL2 c c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 c c ====================================================================== c end file flsh_sub.f c ======================================================================