C*********************************************************************** ! This component has been assigned Type Number 211. If that number conflicts with ! another user Type number, you will need to change it and recompile the appropriate ! dll. SUBROUTINE TYPE211(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) !DEC$ATTRIBUTES DLLEXPORT :: TYPE211 C --------------------------------------------------------------------- C C TYPE 211 : STEAM TO AIR HEATING COIL C C Parameters: C C 1) TSC: degrees of subcooling of condensate (C) C 2) AA: minimum air flow area (m2) C 3) AIS: inside heat transfer area (m2) C 4) AOP: primary (tube) outside area (m2) C 5) AOS: secondary (fin) heat transfer area (m2) C 6) DI: inside tube diameter (m) C 7) DO: outside tube diameter (m) C 8) AKT: tube thermal conductivity [kJ/(kg K)] C 9) DF: fin diameter (m) C 10) FCON: fin thermal conductivity [kJ/(kg K)] C 11) FTH: fin thickness (m) C 12) FPM: number of fins per meter C 13) CM: thermal capacitance of fins and tubes (kJ/K) C C*********************************************************************** DOUBLE PRECISION XIN,OUT C REAL PAR,SAVED REAL PAR,S REAL KA,KW,MUA,MUW,NUSSLT,NTU INTEGER INFO DIMENSION XIN(5),OUT(5),PAR(13),INFO(15) C INTEGER IOSTAT C DIMENSION XIN(5),OUT(4),PAR(13),SAVED(5) C DIMENSION IOSTAT(5) C* COMMON /SOSCOM/ RTOLX,ATOLX,XTOL C* COMMON /CHRONO/ TIME,TSTEP,TTIME,TMIN,ITIME C* COMMON /PROPER/ RHOA,RHOW,CPA,CPW C The RHOA,RHOW,CPA,CPW properties @300 K, 1 atm C The units are: C RHOW - kg/m^3 C CPA,CPW - kJ/kg-K COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) DATA RHOW,CPA,CPW /999,1.007,4.179/ DATA GRAV/9.80665/ ! Set the version information for TRNSYS IF (INFO(7).EQ.-2) THEN INFO(12) = 15 RETURN 1 ENDIF IF (INFO(7) .EQ. 0) THEN TAOI=OUT(5) ENDIF IF (INFO(7) .EQ. -1) THEN INFO(10)=5 CALL TYPECK(1,INFO,5,13,0) TAOI=XIN(4) ENDIF ISAV=INFO(10) PSO= XIN(1) TSI= XIN(2) TAI= XIN(3) TAO= XIN(4) WA= XIN(5) TSC= PAR(1) AA= PAR(2) AIS= PAR(3) AOP= PAR(4) AOS= PAR(5) DI= PAR(6) DO= PAR(7) AKT= PAR(8) DF= PAR(9) FCON= PAR(10) FTH= PAR(11) FPM= PAR(12) CM= PAR(13) TAIR=0.5*(TAI+TAO) TSSAT=TSATS(PSO) ROH=DO/DF RT=(DO-DI)/(2.*AKT*AIS) C Determine outside heat transfer coefficient for air, using correla- C tion from H.Hausen, "Heat Transfer in Counterflow, Parallel Flow and C Cross Flow," McGraw-Hill, 1983. KA=AKA(TAIR) MUA=VISCA(TAIR) CALL CPCVA(TAIR,CPAIR,X,X,X) RE=WA*DO/(MUA*AA) PR=MUA*CPAIR/KA C FFO is the ratio of finned tube surface area to the surface area of C a similar but unfinned tube. FFO=(0.5*(DF*DF-DO*DO)+FTH*(DF-DO)+DO*(1./FPM-FTH))*FPM/DO NUSSLT=0.3*RE**0.625*FFO**(-0.375)*PR**(1./3.) HTCOS=NUSSLT*KA/DO C Determine dry fin efficiency using coefficients saved on first call C* IF (ITIME.LE.1) CALL SUFED(ROH,SAVED) IF (INFO(7).EQ.-1) CALL SUFED(ROH,S,ISAV) FAI=0.5*(DF-DO)*SQRT(2.*HTCOS/(FTH*FCON)) C* FINEFF=SAVED(1)+FAI*(SAVED(2)+FAI*(SAVED(3)+FAI*(SAVED(4) C* & +FAI*SAVED(5)))) FINEFF=S(ISAV)+FAI*(S(ISAV+1)+FAI*(S(ISAV+2)+FAI*(S(ISAV+3) & +FAI*S(ISAV+4)))) RO=1./(HTCOS*(AOS*FINEFF+AOP)) C If inlet steam temp. is below saturation, ignore input pressure C and assume inlet steam is sat. vapor at inlet temp. IF (TSI.LE.TSSAT) TSSAT=TSI DH=HFG(TSSAT) HSL=HSATW(TSSAT) MUW=WMU(TSSAT) KW=WK(TSSAT) C** RHOW=WRHO(TSSAT) C Iterate for inside surface temperature and inside film coefficient RR=0.1 50 TS=TSSAT-RR*(TSSAT-TAIR) HTCIS=0.612*(KW**3*RHOW*RHOW*DH*GRAV/(MUW*DI*(TSSAT-TS)))**0.25 RI=1./(HTCIS*AIS) RSUM=RI+RT+RO RRR=RI/RSUM C* IF (ABS(RRR-RR)/RRR .GT. RTOLX) THEN IF (ABS(RRR-RR)/RRR .GT. 0.01) THEN RR=RRR GOTO 50 ENDIF RR=(RO+RT)/RI NTU=1./(RSUM*WA*CPA) C Find mass flow rate of steam and condensate Q=HTCIS*AIS*(TSSAT-TS) WS=Q/(HS(PSO,TSI)-(HSL-CPW*TSC)) TAOSS=TAI+Q/(WA*CPAIR) TWO=TSSAT-TSC C Air temperature dynamics, based on Myers et al., Trans. ASME, J. Heat C Transfer, May 1970, pp. 269-275. Z=0.5*NTU*(1.+RR)/RR A=2.*NTU*(1.+RR)*EXP(-Z)*SINH(Z)/((SINH(Z)+COSH(Z))*EXP(-Z) & -EXP(-NTU)) TAU=CM/(A*WA*CPAIR) C DTAODT=(TAOSS-TAO)/TAU TAOAA=-1/TAU TAOBB=TAOSS/TAU CALL DIFFEQ(TIME,TAOAA,TAOBB,TAOI,TAOF,TAOBAR) Q=WA*CPAIR*(TAO-TAI) C IOSTAT(1)=0 C IF ((ABS(TAOSS-TAO).LT.RTOLX*TAOSS+ATOLX) .AND. C & (IOSTAT(1) .LT. -1 .OR. IOSTAT(1) .EQ. 2) .AND. C & (IOSTAT(2) .LT. -1 .OR. IOSTAT(2) .EQ. 2) .AND. C & (IOSTAT(3) .LT. -1 .OR. IOSTAT(3) .EQ. 2) .AND. C & (IOSTAT(5) .LT. -1 .OR. IOSTAT(5) .EQ. 2)) C & IOSTAT(1)=1 C IOSTAT(2)=1 C IOSTAT(3)=1 C IOSTAT(4)=1 C IF (ABS(DTAODT).LT.1.E-6) DTAODT=0.0 ! added 2/25/94 C OUT(1)=DTAODT OUT(1)=TAOBAR OUT(2)=TWO OUT(3)=WS OUT(4)=Q OUT(5)=TAOF RETURN 1 END C*********************************************************************** C* SUBROUTINE SUFED(ROH,COF) SUBROUTINE SUFED(ROH,COF,ISAV) C*********************************************************************** INTEGER ISAV DIMENSION X(60),Y(60),COF(5) COMMON /STORE/ NStore,IAV,S(5000) C*** ROH = (OUTSIDE TUBE DIAM.)/(EFFECTIVE FIN DIAM.) DO 5 K=1,5 COF(K)=S(ISAV+K-1) 5 CONTINUE FAI=0.02 DO 10 I=1,60 FAI=FAI+0.035 R1=FAI/(1.-ROH) R2=R1*ROH RO=2.*ROH/(FAI*(1.+ROH)) CALL BESI(R1,1,R1I1,IE1) CALL BESK(R2,1,R2K1,IE2) CALL BESI(R2,1,R2I1,IE3) CALL BESK(R1,1,R1K1,IE4) CALL BESI(R2,0,R2I0,IE5) CALL BESK(R2,0,R2K0,IE6) FED=RO*(R1I1*R2K1-R2I1*R1K1)/(R2I0*R1K1+R1I1*R2K0) X(I)=FAI Y(I)=FED 10 CONTINUE CALL POLFIT(X,Y,60,1.E-05,4,COF) DO 15 K=1,5 S(ISAV+K-1)=COF(K) 15 CONTINUE RETURN END C*********************************************************************** SUBROUTINE BESI(X,N,BI,IER) C ---------------------------------------------------------------------- C C*** SUBROUTINE BESI TO CALCULATE THE MODIFIED BESSEL FUNCTION C*** OF ORDER 0 TO N C*** C*** CALL BESI(X,N,BI,IER) C*** C*** X ARGUMENT OF BESSEL FUNCTION C*** N ORDER OF BESSEL FUNCTION, GREATER THAN OR EQUAL TO ZERO C*** BI RESULTANT VALUE OF I BESSEL FUNCTION C*** IER RESULTANT ERROR CODE: C*** IER = 0 NO ERROR C*** IER = 1 N .LT. 0 C*** IER = 2 X .LT. 0 C*** IER = 3 BI .LT. 10**(-30), BI IS SET TO 0 C*** IER = 4 X .GT. N & X .GT. 90, BI IS SET TO 10**38 C*** C*********************************************************************** DATA PI,TOL/3.141592653,1.E-06/ IER=0 BI=1. IF (X .EQ. 0. .AND. N .EQ. 0) RETURN IF (N .LT. 0) THEN IER=1 RETURN ELSE IF (X .LT. 0.) THEN IER = 2 RETURN ELSE IF (X .GT. 12. .AND. X .GT. N) THEN IF (X .GT. 90.) THEN IER = 4 BI=10.**30 RETURN ENDIF FN=4*N*N XX=0.125/X TERM=1. BI=1. DO 10 K=1,30 IF (ABS(TERM).LE.ABS(TOL*BI)) THEN BI=BI*EXP(X)/SQRT(2.*PI*X) RETURN ENDIF FK=(2*K-1)**2 TERM=TERM*XX*(FK-FN)/REAL(K) BI=BI+TERM 10 CONTINUE ENDIF XX=X/2. TERM=1. IF (N.GT.0) THEN DO 20 I=1,N FI=I IF (ABS(TERM).LT.1.E-30*FI/XX) THEN IER=3 BI=0. RETURN ENDIF TERM=TERM*XX/FI 20 CONTINUE ENDIF BI=TERM XX=XX*XX DO 30 K=1,1000 IF (ABS(TERM).LE.ABS(BI*TOL)) RETURN FK=K*(N+K) TERM=TERM*(XX/FK) BI=BI+TERM 30 CONTINUE RETURN END C*********************************************************************** SUBROUTINE BESK(X,N,BK,IER) C ---------------------------------------------------------------------- C C*** SUBROUTINE BESK TO COMPUTE THE K BESSEL FUNCTION FOR A GIVEN C*** ARGUMENT AND ORDER C*** C*** CALL BESK(X,N,BK,IER) C*** C*** X THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED C*** N THE ORDER OF THE K BESSEL FUNCTION DESIRED C*** BK THE RESULTANT K BESSEL FUNCTION C*** IER RESULTANT ERROR CODE: C*** IER=0 NO ERROR C*** IER=1 N IS NEGATIVE C*** IER=2 X IS ZERO OR NEGATIVE C*** IER=3 X .GT. 85, BK .LT. 10**-38; BK SET TO 0. C*** IER=4 BK .GT. 10**38; BK SET TO 10**38 C*** C*** NOTE: N MUST BE GREATER THAN OR EQUAL TO ZERO C*** C*** METHOD: C*** COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING C*** SERIES APPROXIMATIONS AND THEN COMPUTES N TH ORDER FUNCTION C*** USING RECURRENCE RELATION. C*** RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE C*** AS DESCRIBED BY A.J.M. HITCHCOCK, 'POLYNOMIAL APPROXIMATIONS C*** TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED C*** FUNCTIONS,' M.T.A.C., V.11, 1957, PP. 86-88, AND G.N. WATSON, C*** 'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS,' CAMBRIDGE C*** UNIVERSITY PRESS, 1958, P.62 C*** C*********************************************************************** DIMENSION T(12) DATA GJMAX/ 1.E38/ BK=0. G0=0. GJ=0. IF (N.LT.0) THEN IER=1 RETURN ELSE IF (X.LE.0.) THEN IER=2 RETURN ELSE IF (X.GT.85.) THEN IER=3 BK=0. RETURN ENDIF IER=0 C Use polynomial approximation if x > 1. IF (X.GT.1.) THEN A=EXP(-X) B=1./X C=SQRT(B) T(1)=B DO 10 L=2,12 T(L)=T(L-1)*B 10 CONTINUE IF (N.NE.1) THEN C Compute K0 using polynomial approximation G0=A*(1.2533141-.1566642*T(1)+.08811128*T(2)-.09139095 & *T(3)+.1344596*T(4)-.2299850*T(5)+.3792410*T(6)-.5247277 & *T(7)+.5575368*T(8)-.4262633*T(9)+.2184518*T(10) & -.06680977*T(11)+.009189383*T(12))*C IF (N .EQ. 0) THEN BK=G0 RETURN ENDIF ENDIF C Compute K1 using polynomial approximation G1=A*(1.2533141+.4699927*T(1)-.1468583*T(2)+.1280427*T(3) & -.1736432*T(4)+.2847618*T(5)-.4594342*T(6)+.6283381*T(7) & -.6632295*T(8)+.5050239*T(9)-.2581304*T(10)+.07880001*T(11) & -.01082418*T(12))*C IF (N .EQ. 1) THEN BK=G1 RETURN ENDIF ELSE C Use series expansion if x <= 1. B=X/2. A=.5772157+ALOG(B) C=B*B IF (N.NE.1) THEN C Compute K0 using series expansion G0=-A X2J=1. FACT=1. HJ=0. DO 20 J=1,6 RJ=1./FLOAT(J) X2J=X2J*C FACT=FACT*RJ*RJ HJ=HJ+RJ G0=G0+X2J*FACT*(HJ-A) 20 CONTINUE IF (N.EQ.0) THEN BK=G0 RETURN ENDIF ENDIF C Compute K1 using series expansion X2J=B FACT=1. HJ=1. G1=1./X+X2J*(.5+A-HJ) DO 30 J=2,8 X2J=X2J*C RJ=1./FLOAT(J) FACT=FACT*RJ*RJ HJ=HJ+RJ G1=G1+X2J*FACT*(.5+(A-HJ)*FLOAT(J)) 30 CONTINUE IF (N.EQ.1) THEN BK=G1 RETURN ENDIF ENDIF C From K0 and K1 compute KN using recurrence relation DO 40 J=2,N GJ=2.*(FLOAT(J)-1.)*G1/X+G0 IF (GJ-GJMAX) 33,33,32 32 IER=4 GJ=GJMAX GOTO 34 33 G0=G1 G1=GJ 40 CONTINUE 34 BK=GJ RETURN END C*********************************************************************** SUBROUTINE POLFIT(X,Y,N,TOL,LAST,COF) C --------------------------------------------------------------------- C C*** SUBROUTINE POLFIT FITS POLYNOMIAL OF ORDER FROM 1 TO LAST TO THE C*** ORDERED PAIRS OF DATA POINTS X,Y C*** C*********************************************************************** DIMENSION X(60),Y(60), A(10,10),SUMX(10),SUMY(10),COF(10) SUMX(1)=N SUMX(2)=0. SUMX(3)=0. SUMY(1)=0. SUMY(2)=0. L=LAST+1 DO 80 I=1,L 80 COF(I)=0. DO 90 I=1,N SUMX(2)=SUMX(2)+X(I) SUMX(3)=SUMX(3)+X(I)*X(I) SUMY(1)=SUMY(1)+Y(I) SUMY(2)=SUMY(2)+X(I)*Y(I) 90 CONTINUE NORD=1 91 L=NORD+1 KK=L+1 DO 101 I=1,L DO 100 J=1,L IK=J-1+I A(I,J)=SUMX(IK) 100 CONTINUE A(I,KK)=SUMY(I) 101 CONTINUE DO 140 I=1,L A(KK,I)=-1. KKK=I+1 DO 110 J=KKK,KK A(KK,J)=0. 110 CONTINUE C=1./A(1,I) DO 120 II=2,KK DO 120 J=KKK,KK A(II,J)=A(II,J)-A(1,J)*A(II,I)*C 120 CONTINUE DO 140 II=1,L DO 140 J=KKK,KK A(II,J)=A(II+1,J) 140 CONTINUE S2=0. DO 160 J=1,N S1=A(1,KK) DO 150 I=1,NORD S1=S1+A(I+1,KK)*X(J)**I 150 CONTINUE S2=S2+(S1-Y(J))*(S1-Y(J)) 160 CONTINUE B=N-L IF (S2 .GT. 0.0001) S2=SQRT(S2/B) DO 164 I=1,L J=I-1 COF(I)=A(I,KK) 164 CONTINUE IF (NORD-LAST) 170,173,173 170 IF (S2-TOL) 173,173,171 171 NORD=NORD+1 J=2*NORD SUMX(J)=0. SUMX(J+1)=0. SUMY(NORD+1)=0. DO 172 I=1,N SUMX(J)=SUMX(J)+X(I)**(J-1) SUMX(J+1)=SUMX(J+1)+X(I)**J SUMY(NORD+1)=SUMY(NORD+1)+Y(I)*X(I)**NORD 172 CONTINUE GOTO 91 173 CONTINUE RETURN END C*********************************************************************** SUBROUTINE CPCVA(TCA,CPA,CVA,GAMMA,SONIC) C ---------------------------------------------------------------------- C C This subroutine takes Celsius air temperature, TCA, and computes: C CPA: Specific heat of air at constant pressure [KJ/(kg K)] C CVA: Specific heat of air at constant volume [KJ/(kg K)] C GAMMA: The ratio CPA/CVA [dimensionless] C SONIC: Speed of sound in air (M/S) C C*********************************************************************** DATA A0,A1,A2/1.03409,-0.284887E-3,0.7816818E-6/ DATA A3,A4,TCONV,R/-0.4970786E-9,0.1077024E-12,273.15,0.287040/ T=TCA+TCONV CPA=A0+T*(A1+T*(A2+T*(A3+T*A4))) CVA=CPA-R GAMMA=CPA/CVA SONIC=SQRT(GAMMA*R*T) RETURN END C*********************************************************************** FUNCTION AKA(TC) C ---------------------------------------------------------------------- C C Thermal conductivity of air [KW/(M K)], given Celsius temperature C C*********************************************************************** DATA C0,C1,C2/-2.276501E-3,1.2598485E-4,-1.4815235E-7/ DATA C3,C4,C5/1.73550646E-10,-1.066657E-13,2.47663035E-17/ DATA TCONV/273.15/ T=TC+TCONV AKA=0.001*(C0+T*(C1+T*(C2+T*(C3+T*(C4+T*C5))))) RETURN END C*********************************************************************** FUNCTION HFG(TC) C ---------------------------------------------------------------------- C C Latent heat of vaporization of water (KJ/KG) given sat. T (C) C C*********************************************************************** DATA E1,E2,E3,E4,E5/-3.87446,2.94553,-8.06395,11.5633,-6.02884/ DATA B,C,D,HFGTP,TCR/0.779221,4.62668,-1.07931,2500.9,647.3/ DATA TCNV/273.15/ IF (TC.LT.0.) TC=0. TR=(TCR-TC-TCNV)/TCR IF (TR.LT.0.) THEN HFG=0. RETURN ENDIF Y=B*TR**(1./3.)+C*TR**(5./6.)+D*TR**0.875 Y=Y+TR*(E1+TR*(E2+TR*(E3+TR*(E4+TR*E5)))) HFG=Y*HFGTP RETURN END C*********************************************************************** FUNCTION HS(PKPA,TC) C ---------------------------------------------------------------------- C C Enthalpy of superheated steam (KJ/KG) given P (KPA) and T (C) C C*********************************************************************** DATA B11,B12,B13,B21/2041.21,-40.4002,-0.48095,1.610693/ DATA B22,B23,B31/5.472051E-2,7.517537E-4,3.383117E-4/ DATA B32,B33,B41,B42/-1.975736E-5,-2.87409E-7,1707.82,-16.99419/ DATA B43,B44,B45,EM/6.2746295E-2,-1.0284259E-4,6.4561298E-8,45./ DATA C1,C2,C3,PCNV/42.6776,-3892.70,-9.48654,0.001/ DATA C4,C5,C6,TCNV/-387.592,-12587.5,-15.2578,273.15/ P=PKPA*PCNV T=TC+TCNV TS=C1+C2/(ALOG(P)+C3) IF (P .GE. 12.33) TS=C4+C5/(ALOG(P)+C6) A0=B11+P*(B12+P*B13) A1=B21+P*(B22+P*B23) A2=B31+P*(B32+P*B33) A3=B41+TS*(B42+TS*(B43+TS*(B44+TS*B45))) HS=A0+T*(A1+T*A2)-A3*EXP((TS-T)/EM) RETURN END C*********************************************************************** FUNCTION HSATW(TC) C ---------------------------------------------------------------------- C C Sat. enthalpy of liquid water (KJ/KG) given sat. T (C) C C*********************************************************************** DATA E11,E21,E31,HFCR/624.698837,-2343.85369,-9508.12101,2099.3/ DATA E41,E51,E61,TCNV/71628.7928,-163535.221,166531.093,273.15/ DATA E71,A2,E12/-64785.4585,0.8839230108,-2.67172935/ DATA E22,E32,E42/6.22640035,-13.1789573,-1.91322436/ DATA E52,E62,E72,A3/68.7937653,-124.819906,72.1435404,1.0/ DATA B3,C3,D3/-0.441057805,-5.52255517,6.43994847/ DATA E13,E23,TCR/-1.64578795,-1.30574143,647.3/ TK=TC+TCNV TR=(TCR-TK)/TCR IF (TK.LT.300)THEN Y=TR*(E11+TR*(E21+TR*(E31+TR*(E41+TR*(E51+TR*(E61+TR*E71)))))) ELSE IF (TK.LT.600)THEN Y=TR*(E12+TR*(E22+TR*(E32+TR*(E42+TR*(E52+TR*(E62+TR*E72)))))) Y=Y+A2 ELSE Y=A3+B3*TR**(1./3.)+C3*TR**(5./6.)+D3*TR**0.875+TR*(E13+TR*E23) ENDIF HSATW=Y*HFCR RETURN END C*********************************************************************** FUNCTION TSATS(PKPA) C ---------------------------------------------------------------------- C C Saturation temp. of steam (C) as a function of pressure (KPA) C C*********************************************************************** DATA A1,B1,C1,TCONV/42.6776,-3892.70,-9.48654,-273.15/ DATA A2,B2,C2,PCONV/-387.592,-12587.5,-15.2578,0.001/ P=PKPA*PCONV IF (P .LT. 12.33) THEN TSATS=TCONV+A1+B1/(ALOG(P)+C1) ELSE TSATS=TCONV+A2+B2/(ALOG(P)+C2) ENDIF RETURN END C*********************************************************************** FUNCTION VISCA(TC) C ---------------------------------------------------------------------- C C Dynamic viscosity [(N S)/(M*M)] of air, from celsius temperature C C*********************************************************************** DATA A0,A1,A2/-0.98601,9.080125E-2,-1.17635575E-4/ DATA A3,A4/1.2349703E-7,-5.7971299E-11/ DATA B0,B1,B2/4.8856745,5.43232E-2,-2.4261775E-5/ DATA B3,B4,TCONV/7.9306E-9,-1.10398E-12,273.15/ T=TC+TCONV VISCA=A0+T*(A1+T*(A2+T*(A3+T*A4))) IF (T .GE. 600.) VISCA=B0+T*(B1+T*(B2+T*(B3+T*B4))) VISCA=VISCA*1.E-6 RETURN END C*********************************************************************** FUNCTION WK(TW) C ---------------------------------------------------------------------- C C Thermal conductivity equation from linear least-squares fit to data C in CRC Handbook (op.cit.), page E-11; temps. from 270 K to 620 K. C Temperature in Celsius, WK in [kW/(m K)]. Values at one atmosphere C for T from 0 to 100 C, at saturation for T above 100. C C*********************************************************************** DATA AK0,AK1,AK2/0.560101,0.00211703,-1.05172E-05/ DATA AK3,AK4/1.497323E-08,-1.48553E-11/ WK=0.001*(AK0+TW*(AK1+TW*(AK2+TW*(AK3+TW*AK4)))) RETURN END C*********************************************************************** FUNCTION WMU(TW) C ---------------------------------------------------------------------- C C Viscosity equations for water at 1 atm., from CRC Handbook (op.cit.), C page F-51. WMU in kg/meter-sec; for centipoise, multiply by 1000. C For temps > 100 C, fit to data from Karlekar & Desmond (saturated). C C*********************************************************************** DATA AM0,AM1,AM2,AM3,AM4/-3.30233,1301.,998.333,8.1855,.00585/ DATA AM5,AM6,AM7,AM8/1.002,-1.3272,-0.001053,105./ DATA A10,A11,A12,A13/.68714,-.0059231,2.1249E-05,-2.69575E-08/ WMU=AM5*10.**((TW-20.)*(AM6+(TW-20.)*AM7)/(TW+AM8)) IF (TW.LT.20.) & WMU=10.**(AM0+AM1/(AM2+(TW-20.)*(AM3+AM4*(TW-20.))))*100. IF (TW.GT.100.)WMU=A10+TW*(A11+TW*(A12+TW*A13)) WMU=0.001*WMU RETURN END