C*********************************************************************** SUBROUTINE TYPE15(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C ---------------------------------------------------------------------- C C TYPE 15 : ROOM MODEL C C*********************************************************************** DOUBLE PRECISION XIN,OUT REAL PAR,S INTEGER INFO DIMENSION XIN(8),OUT(10),PAR(6),INFO(15) C REAL PAR,SAVED C INTEGER IOSTAT C DIMENSION XIN(8),OUT(6),PAR(6),SAVED(9) C DIMENSION IOSTAT(8) C* COMMON /PROPER/ RHOA,RHOW,CPA,CPW C* COMMON /CHRONO/ TIME,TSTEP,TTIME,TMIN,ITIME C The RHOA,CPA properties @300 K, 1 atm C The units are: C RHOA - kg/m^3 C CPA - kJ/kg-K COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) DATA RHOA,CPA /1.1614,1.007/ IF (INFO(7) .EQ. 0) THEN T1I=OUT(7) TM1I=OUT(8) TM2I=OUT(9) T2BI=OUT(10) ENDIF IF (INFO(7) .EQ. -1) THEN INFO(10)=9 CALL TYPECK(1,INFO,8,6,0) T1I=XIN(3) TM1I=XIN(4) TM2I=XIN(5) T2BI=XIN(6) ENDIF ISAV=INFO(10) WA= XIN(1) TI= XIN(2) T1= XIN(3) TM1= XIN(4) TM2= XIN(5) T2B= XIN(6) Q1= XIN(7) Q2= XIN(8) VOL= PAR(1) WCAP=PAR(2) RCAP=PAR(3) H1= PAR(4) H2= PAR(5) F= PAR(6) C IF (SAVED(9).GT.TIME)SAVED(9)=0. IF (S(ISAV+8) .GT. TIME) S(ISAV+8)=0. F1=1.-F C DO 100 I=1,6 C 100 IOSTAT(I)=0 ACAP=VOL*RHOA*CPA C DT1DT=(Q2/F+H1*(TM1-T1)+H2*(TM2-T1))/ACAP T1AA=-(H1+H2)/ACAP T1BB=((Q2/F)+(H1*TM1)+(H2*TM2))/ACAP CALL DIFFEQ(TIME,T1AA,T1BB,T1I,T1F,T1BAR) TS=(H1*TM1+H2*TM2)/(H1+H2) IF (WA.GT.1.E-6) GOTO 1000 TO=TS TBAR=TS TAU2=F1*(H1+H2)/ACAP GOTO 2000 1000 TAU1=VOL*RHOA/WA TAU2=F1*TAU1 TAU1=F*TAU1 DT=TS-T1 IF (ABS(DT).LT.1.E-10) DT=0. A=F1*(H1+H2)/(WA*CPA) E=0. IF (A.LT.20.) E=EXP(-A) T2=TS-DT*E C TO=DELAY(T2,TAU2,SAVED(1),SAVED(7),SAVED(8),SAVED(9)) TO=DELAY(T2,TAU2,S(ISAV),S(ISAV+6),S(ISAV+7),S(ISAV+8),INFO) C DT1DT=DT1DT+(TI-T1)/TAU1 T1AA=(-1/TAU1)-((H1+H2)/ACAP) T1BB=(((Q2/F)+(H1*TM1)+(H2*TM2))/ACAP)+(TI/TAU1) CALL DIFFEQ(TIME,T1AA,T1BB,T1I,T1F,T1BAR) E1A=1. IF (E.LT.1.0) E1A=(1.0-E)/A TBAR=TS-DT*E1A 2000 TR=F*T1+F1*T2B C DTM1DT=(H1*(TR-TM1)+Q1)/WCAP C DTM2DT=H2*(TR-TM2)/RCAP C DT2BDT=(TBAR-T2B)/TAU2 TM1AA=(-H1/WCAP) TM1BB=(Q1+(H1*TR))/WCAP CALL DIFFEQ(TIME,TM1AA,TM1BB,TM1I,TM1F,TM1BAR) TM2AA=-H2/RCAP TM2BB=(H2*TR)/RCAP CALL DIFFEQ(TIME,TM2AA,TM2BB,TM2I,TM2F,TM2BAR) T2BAA=-1/TAU2 T2BBB=TBAR/TAU2 CALL DIFFEQ(TIME,T2BAA,T2BBB,T2BI,T2BF,T2BBAR) C IF (ABS(DT1DT).LT.1.E-6) DT1DT =0.0 ! added 2/24/94 C IF (ABS(DTM1DT).LT.1.E-6) DTM1DT=0.0 C IF (ABS(DTM2DT).LT.1.E-6) DTM2DT=0.0 C IF (ABS(DT2BDT).LT.1.E-6) DT2BDT=0.0 C OUT(1)= DT1DT C OUT(2)= DTM1DT C OUT(3)= DTM2DT C OUT(4)= DT2BDT OUT(1)=T1BAR OUT(2)=TM1BAR OUT(3)=TM2BAR OUT(4)=T2BBAR OUT(5)= TR OUT(6)= TO OUT(7)=T1F OUT(8)=TM1F OUT(9)=TM2F OUT(10)=T2BF 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