C*********************************************************************** SUBROUTINE TYPE11(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C ---------------------------------------------------------------------- C C TYPE 11 : HOT WATER COIL MODEL (CROSSFLOW) C C*********************************************************************** DOUBLE PRECISION XIN,OUT REAL PAR,S C REAL PAR,SAVED REAL KA,KW,NTU,LAG INTEGER INFO DIMENSION XIN(9),OUT(5),PAR(8),INFO(15) C LOGICAL FLO1,FLO2,TEM1,TEM2,FRZDE C INTEGER IOSTAT C DIMENSION XIN(9),OUT(5),PAR(8),SAVED(11) C DIMENSION IOSTAT(9) C* COMMON /CHRONO/ TIME,TSTEP,TTIME,TMIN,ITIME C* COMMON /PROPER/ RHOA,RHOW,CPA,CPW C* COMMON /SOSCOM/ RTOLX,ATOLX,XTOL C* COMMON /XINIT/ INIT,NSAVED 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 INIT /0/ IF (INFO(7) .EQ. -1) THEN INFO(10)=11 CALL TYPECK(1,INFO,9,8,0) ENDIF ISAV=INFO(10) PAO= XIN(1) WA= XIN(2) TAI= XIN(3) TAO= XIN(4) PWO= XIN(5) WW= XIN(6) TWI= XIN(7) TWO= XIN(8) TD= XIN(9) HAA=PAR(1)*WA**PAR(2) HAW=PAR(3)*WW**PAR(4) KA= PAR(5) KW= PAR(6) VOL=PAR(7) CM= PAR(8) IF (WA.LT.1.E-6) WA=0. ! added 2/24/94 IF (WW.LT.1.E-6) WW=0. C* IF (ITIME.GT.1) GOTO 50 IF (INFO(7) .GT. -1) GOTO 50 C IF (INIT.EQ.0 .OR. SAVED(9).GT.TIME) SAVED(9)=0. C IF (INIT.EQ.0)SAVED(11)=TWI C 50 IF (TIME.GT.SAVED(9)) SAVED(10)=SAVED(11) IF (INIT .EQ. 0 .OR. S(ISAV+8) .GT. TIME) S(ISAV+8)=0. IF (INIT .EQ. 0) S(ISAV+10)=TWI 50 IF (TIME .GT. S(ISAV+8)) S(ISAV+9)=S(ISAV+10) C FLO1=(IOSTAT(2).LT.-1).OR.(IOSTAT(2).EQ.2) C TEM1=(IOSTAT(3).LT.-1).OR.(IOSTAT(3).EQ.2) C FLO2=(IOSTAT(6).LT.-1).OR.(IOSTAT(6).EQ.2) C TEM2=(IOSTAT(7).LT.-1).OR.(IOSTAT(7).EQ.2) C FRZDE=FLO1.AND.FLO2.AND.TEM1.AND.TEM2 SS=1. IF (WA.NE.0.) SS=SIGN(1.,WA) PAI=PAO+SS*KA*WA**2 SS=1. IF (WW.NE.0.) SS=SIGN(1.,WW) PWI=PWO+SS*KW*WW**2 C* RATE=1./TSTEP RATE=1./DELT IF (WA.GT.0. .AND. WW.GT.0.) GOTO 200 DTAODT=0. DTDDT=0. DTWODT=0. GOTO 3000 200 CA=CPA*WA CW=CPW*WW CMIN=AMIN1(CA,CW) CMAX=AMAX1(CA,CW) UA=HAA*HAW/(HAA+HAW) C CALCULATE EFFECTIVENESS: UNMIXED BOTH SIDES, APPROX. EQUATION. NTU=UA/CMIN ETA=NTU**0.22 R=CMIN/CMAX E1=0. E2=0. A=R*NTU/ETA IF (A.GT.20.) GOTO 2000 E1=EXP(-A) 2000 A=ETA*(1.-E1)/R IF (A.GT.20.) GOTO 2001 E2=EXP(-A) 2001 EFEC=1.-E2 C--the next 3 lines are effectiveness for fully mixed fluids C** IF (NTU.LE.20.) E1=EXP(-NTU) C** IF (R*NTU.LE.20.) E2=EXP(-R*NTU) C** EFEC=NTU/(NTU/(1.-E1)+R*NTU/(1.-E2)-1.) TAOSS=TAI+EFEC*CMIN*(TWI-TAI)/CA TWOSS=TWI-CA*(TAOSS-TAI)/CW C DTWIDT=(TWI-SAVED(10))*RATE DTWIDT=(TWI-S(ISAV+9))*RATE TAUX=RHOW*VOL/WW B1=HAW/CW B2=TAUX*HAA/CM B3=TAUX*HAW/CM B4=HAA/CA A=B3+2.*B2/(2.+B4) PHI=B1-B1*B3/A E=EXP(-PHI) DPHI=1.+B1*B3/(A*A) c* Z=(A+B1)/(A*PHI)-DPHI*E/(1.-E) E1=E/(1.-E) AB1=A+B1 APHI=A*PHI A2=AB1/APHI-DPHI*E1 A3=A2*A2-(AB1*AB1-APHI)/(APHI*APHI) & +(DPHI*(AB1*(2.+PHI)+PHI*(2.-PHI))-2.*PHI)*E1/(2.*APHI) DTDDT=(TAOSS-TAO)/TAUX DTAODT=(TD-A2*TAO)/(A3*TAUX) TAUW=TAUX*EXP(B1*B3/(2.*A))/A LAG=DTWIDT*EXP(-B1)*TAUW+(TWOSS-TAI) C DTWODT=((TAI-TWO)+DELAY(LAG,TAUX,SAVED(1),SAVED(7), C & SAVED(8),SAVED(9)))/TAUW DTWODT=((TAI-TWO)+DELAY(LAG,TAUX,S(ISAV),S(ISAV+6), & S(ISAV+7),S(ISAV+8),INFO))/TAUW C 3000 SAVED(11)=TWI 3000 S(ISAV+10)=TWI C IF (ABS(DTAODT).LT.ATOLX) DTAODT=0.0 ! added 2/24/94 C IF (ABS(DTWODT).LT.ATOLX) DTWODT=0.0 C IF (ABS(DTDDT).LT.ATOLX) DTDDT=0.0 IF (ABS(DTAODT) .LT. 0.01) DAODT=0.0 IF (ABS(DTWODT) .LT. 0.01) DTWODT=0.0 IF (ABS(DTDDT) .LT. 0.01) DTDDT=0.0 OUT(1)=DTAODT OUT(2)=DTWODT OUT(3)=DTDDT OUT(4)=PAI OUT(5)=PWI C IOSTAT(4)=1 C IOSTAT(5)=1 C IOSTAT(1)=0 C IOSTAT(2)=0 C IOSTAT(3)=0 C EITHER ALL THREE D.E.'S FREEZE, OR NONE OF THEM FREEZE. C IF ( (ABS(DTAODT) .LE. ATOLX) .AND. C & (ABS(DTDDT) .LE. ATOLX) .AND. C & (ABS(DTWODT) .LE. ATOLX) .AND. FRZDE) THEN C IOSTAT(1)=1 C IOSTAT(2)=1 C IOSTAT(3)=1 C ENDIF RETURN 1 END C*********************************************************************** C FUNCTION DELAY(TINEW,TAU,A,TIN,DZOLD,TIME0) FUNCTION DELAY(TINEW,TAU,A,TIN,DZOLD,TIMEZERO,INFO) C --------------------------------------------------------------------- C C Models transport delays in ducting components. C The temperature distribution in the duct is modeled by a fifth C order time dependent polynomial. The coefficients of the C polynomial are calculated at each time step. C C Modified: March 5, 1987 C.P. C C*********************************************************************** C PARAMETER (M=5,N=M+1) C DIMENSION A(N),AM(M,N) DIMENSION A(6),AM(5,6),INFO(15) C COMMON /CHRONO/ TIME,TSTEP,TTIME,TMIN,ITIME C COMMON /SOSCOM/ RTOLX,ATOLX,XTOL C COMMON /XINIT/ INIT,NSAVED COMMON /SIM/ TIME0,TFINAL,DELT,IWARN DATA INIT,M,N /0,5,6/ C Initialize coefficients if this is the first call. C If the transport delay is zero, return the input value. C IF (ITIME.EQ.1 .AND. INIT.EQ.0) THEN IF (INFO(7) .EQ. -1 .AND. INIT .EQ. 0) THEN TIN=TINEW A(1)=TIN C TIME0=TIME TIMEZERO=TIME DELAY=TIN IF (TAU.GT.0.) THEN C DZOLD=TSTEP/TAU DZOLD=DELT/TAU ELSE DZOLD=2.0 ENDIF DO 10 I=2,N A(I)=0. 10 CONTINUE RETURN ENDIF C Set the time step equal to minimum time step if TSTEP < TMIN C IF (TSTEP.LT.TMIN) THEN C TSTEP=TMIN C ENDIF IF (TAU.GT.0.) THEN IF (TAU.GT.1.0E+10) THEN TAU=1.0E+10 ENDIF C DZ=TSTEP/TAU DZ=DELT/TAU ELSE DZ=2.0 ENDIF C If first call of the timestep, evaluate new coefficients. C IF (TIME.GT.TIME0) THEN IF (TIME .GT. TIMEZERO) THEN C Set up matrix and constant vector to solve for new set of C coefficients. IF (DZOLD .LE. 0.5) THEN DO 30 I=1,M AM(I,N)=0. X=1. XX=1. IF (I.EQ.1) THEN IF (DZOLD.GE.0.2) THEN Z=0.5*DZOLD ELSE Z=DZOLD ENDIF ELSEIF (I.EQ.2 .AND. DZOLD.GE.0.2) THEN Z=DZOLD ELSE I2=I Z=0.2*REAL(I2) ENDIF DO 20 J=1,M XX=XX*Z AM(I,J)=XX IF (Z.GT.DZOLD) THEN X=X*(Z-DZOLD) AM(I,N)=AM(I,N)+A(J+1)*X ENDIF 20 CONTINUE IF (Z.LE.DZOLD) THEN AM(I,N)=(A(1)-TIN)*Z/DZOLD ELSE AM(I,N)=AM(I,N)+A(1)-TIN ENDIF 30 CONTINUE C Solve for coefficients by Gaussian elimination. DO 40 I=2,M DO 40 J=I,M R=AM(J,I-1)/AM(I-1,I-1) DO 40 K=I,N AM(J,K)=AM(J,K)-R*AM(I-1,K) 40 CONTINUE DO 50 I=2,M K=N+1-I R=AM(K,N)/AM(K,K) DO 50 J=I,M L=N-J AM(L,N)=AM(L,N)-R*AM(L,K) 50 CONTINUE A(1)=TIN DO 60 I=1,M A(I+1)=AM(I,N)/AM(I,I) 60 CONTINUE C Last DZ between 0.5 and 1: reset to linear temperature C distribution. ELSEIF (DZOLD .GT. 0.5 .AND. DZOLD .LT. 1.0) THEN DO 70 I=2,N A(I)=0. 70 CONTINUE A(2)=A(1)-TIN A(1)=TIN C Last DZ beyond inlet: reset to uniform temperature. ELSE DO 80 I=2,N A(I)=0. 80 CONTINUE A(1)=TIN ENDIF ENDIF C Calculate output. C TIME0=TIME TIMEZERO=TIME TIN=TINEW DZOLD=DZ IF (DZ.GE.1.0) THEN DELAY=TINEW ELSE DELAY=A(1) X=1. DO 90 I=2,N X=X*(1.-DZ) DELAY=DELAY+A(I)*X 90 CONTINUE ENDIF RETURN END