SUBROUTINE TYPE71 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE A Toolkit for Primary HVAC System Energy C* Calculation C*********************************************************************** C* PROGRAM: TYPE71 (ABSCHIID) C* C* LANGAGE: FORTRAN 77 C* C* PURPOSE: Parameter identification based on C* absorption chiller performances in C* steady-state regime. C*********************************************************************** C* INPUT VARIABLES C* N Number of available working points (-) C* xin(1) (-) C* Choice If Choice (-) C* =1: the absorber and the condenser are set in series C* =2: the absorber and the condenser are not set in C* series C* xin(2) (-) C* Sum If Sum C* =0: the heats rejected out of the condenser (-) C* and absorber are known C* =1: the sum of the heats rejected out of both C* condenser and absorber is known C* xin(3) (-) C* AUevi Lower bound of the evaporator heat transfer (W/K) C* coefficient C* xin(4) (kJ/hr/øC) C* AUevf Upper bound of the evaporator heat transfer (W/K) C* coefficient C* xin(5) (kJ/hr/øC) C* dAUev Increment (W/K) C* xin(6) (kJ/hr/øC) C* AUcdi Lower bound of the condenser heat transfer (W/K) C* coefficient C* xin(7) (kJ/hr/øC) C* AUcdf Upper bound of the condenser heat transfer (W/K) C* coefficient C* xin(8) (kJ/hr/øC) C* dAUcd Increment (W/K) C* xin(9) (kJ/hr/øC) C* AUabsi Lower bound of the absorber heat transfer (W/K) C* coefficient C* xin(10) (kJ/hr/øC) C* AUabsf Upper bound of the absorber heat transfer (W/K) C* coefficient C* xin(11) (kJ/hr/øC) C* dAUabs Increment (W/K) C* xin(12) (kJ/hr/øC) C* AUgeni Lower bound of the generator heat transfer (W/K) C* coefficient C* xin(13) (kJ/hr/øC) C* AUgenf Upper bound of the generator heat transfer (W/K) C* coefficient C* xin(14) (kJ/hr/øC) C* dAUgen Increment (W/K) C* xin(15) (kJ/hr/øC) C* AUexchi Lower bound of the heat-exchanger heat transfer (W/K) C* coefficient C* xin(16) (kJ/hr/øC) C* AUexchf Upper bound of the heat-exchanger heat transfer (W/K) C* coefficient C* xin(17) (kJ/hr/øC) C* dAUexch Increment (W/K) C* xin(18) (kJ/hr/øC) C* Vi Lower bound of the volume flow rate (m**3/s) C* of the pump of the absorption compressor C* xin(19) (m**3/hr) C* Vf Upper bound of the volume flow rate (m**3/s) C* of the pump of the absorption compressor C* xin(20) (m**3/hr) C* dV Increment (m**3/s) C* xin(21) (m**3/hr) C* C* C* For the working point considered : C* Mfrwev Evaporator water mass flow rate (kg/s) C* xin(22) (kg/hr) C* Mfrwabs Absorber water mass flow rate (kg/s) C* xin(23) (kg/hr) C* Mfrwcd Condenser water mass flow rate (kg/s) C* ( = 0 if Choice = 1 ) C* xin(24) (kg/hr) C* Twsuev Evaporator supply water temperature (K) C* xin(25) (øC) C* Twsuabs Absorber supply water temperature (K) C* xin(26) (øC) C* Twsucd Condenser supply water temperature (K) C* ( = 0 if Choice = 1 ) C* xin(27) (øC) C* psteam Steam pressure (Pa) C* xin(28) (atm) C* Pev Cooling capacity (W) C* xin(29) (kJ/hr) C* Pcd Heat rejected in the condenser ( = 0 if Sum = 1 ) (W) C* xin(30) (kJ/hr) C* Pabs Heat rejected out of the absorber ( = 0 if Sum=1) (W) C* xin(31) (kJ/hr) C* Pcdpabs Sum of the heats rejected out of both condenser (W) C* and absorber ( = 0 if Sum = 0 ) C* xin(32) (kJ/hr) C* C* OUTPUT VARIABLES C* AUev Evaporator heat transfer coefficient (W/K) C* out(1) (kJ/hr/øC) C* AUcd Condenser heat transfer coefficient (W/K) C* out(2) (kJ/hr/øC) C* AUabs Absorber heat transfer coefficient (W/K) C* out(3) (kJ/hr/øC) C* AUgen Generator heat transfer coefficient (W/K) C* out(4) (kJ/hr/øC) C* AUexch Heat-exchanger heat transfer coefficient (W/K) C* out(5) (kJ/hr/øC) C* VLow Volume flow rate of the pump of the (m**3/s) C* absorption compressor C* out(6) (m**3/hr) C* Fmin Minimal value of function F (-) C* out(7) (-) C* ErrDetec This variable is equal to 1 if the routine does (-) C* not converge, is equal to 2 if cristallization C* of the high-concentrated solution occurs at C* the heat-exchanger exhaust and is equal to 3 if C* the determinant of the matrix to be inversed C* is equal to 0. C* out(8) (-) C* PevCal Cooling capacity calculated by means of the (W) C* identified parameters for the working point C* considered C* out(9) (kJ/hr) C* Pev Actual cooling capacity (W) C* for the working point considered C* out(10) (kJ/hr) C* PcdCal Heat rejected in the condenser calculated by (W) C* means of the identified parameters for the C* working point considered ( = 0 if Sum = 1) C* out(11) (kJ/hr) C* Pcd Actual heat rejected in the condenser (W) C* for the working point considered ( = 0 if Sum = 1) C* out(12) (kJ/hr) C* PabsCal Heat rejected out of the absorber calculated (W) C* by means of the identified parameters for the C* working point considered ( = 0 if Sum = 1) C* out(13) (kJ/hr) C* Pabs Actual heat rejected out of the absorber (W) C* for the working point considered ( = 0 if Sum = 1) C* out(14) (kJ/hr) C* PcdpabsCal Sum of the heats rejected out of both condenser (W) C* and absorber calculated by means of the identified C* parameters for the working point considered C* ( = 0 if Sum = 0 ) C* out(15) (kJ/hr) C* Pcdpabs Sum of the actual heats rejected out of both (W) C* condenser and absorber ( = 0 if Sum = 0 ) C* out(16) (kJ/hr) C* C*********************************************************************** C MAJOR RESTRICTIONS: The surrounding heat exchanges are C neglected. The refrigerant leaves the C evaporator and the condenser as saturated C vapor and saturated liquid respectively. C The steam is assumed to be C saturated. Perfect fluid properties C are used.The absorbent properties are C valid in the following range of mass C concentration of LiBr: 45 % to 70 % . C The compression of the solution is assumed C to be isenthalpic. C The maximum number of working points is C equal to 100. C C DEVELOPERS: Jean Lebrun C Bruno Boulanger C Jean-Pascal Bourdouxhe C Marc Grodent C University of LiŠge,Belgium C C DATE: March 1, 1995 C C SUBROUTINES CALLED: ABSFLSIM C LINKCK C C Modified for TRNSYS: January 29, 1995 C Jean-Pascal Bourdouxhe C*********************************************************************** C INTERNAL VARIABLES: C F Function to minimize C Err Array containing the value of the error detection C variable corresponding to each working point C C Fmin is a storage variable and Id is a storage array C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 INFO,INFO73 DOUBLE PRECISION XIN,OUT,XIN73,OUT73 REAL*8 MfrSteam,Mfrwev(100),Mfrwcd(100),Mfrwabs(100),Id(406) REAL*4 TIME,PAR73 DIMENSION Twsuev(100),Twsuabs(100),Twsucd(100),psteam(100), & Pev(100),Pevcal(100),Pcd(100),PcdCal(100),Pabs(100), & PabsCal(100),Pcdpabs(100),PcdpabsCal(100),Err(100) DIMENSION XIN(32),OUT(16),INFO(15),PAR73(6),XIN73(11), & OUT73(13),INFO73(15) COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, & PRWARN INFO(6)=16 INFO73(6)=13 C2*** Initialization ErrDetec=0. C1*** Number of working points to be treated and boundaries of the C1*** search method ITIME=INT(TIME) N=IDINT(xin(1)) N2=2*N Choice=xin(2) Sum=xin(3) AUevi=xin(4)/3.6 AUevf=xin(5)/3.6 dAUev=xin(6)/3.6 AUcdi=xin(7)/3.6 AUcdf=xin(8)/3.6 dAUcd=xin(9)/3.6 AUabsi=xin(10)/3.6 AUabsf=xin(11)/3.6 dAUabs=xin(12)/3.6 AUgeni=xin(13)/3.6 AUgenf=xin(14)/3.6 dAUgen=xin(15)/3.6 AUexchi=xin(16)/3.6 AUexchf=xin(17)/3.6 dAUexch=xin(18)/3.6 Vi=xin(19)/3600. Vf=xin(20)/3600. dV=xin(21)/3600. IF (ITIME.LE.N) THEN C1*** FIRST PART : reception of the remaining inputs associated with C1*** the working point considered. These inputs are stored in arrays. C2*** REMAINING INPUTS 12 (converted in SI units) C2*** ******************* Mfrwev(ITIME)=xin(22)/3600. Mfrwabs(ITIME)=xin(23)/3600. Mfrwcd(ITIME)=xin(24)/3600. Twsuev(ITIME)=xin(25)+273.15 Twsuabs(ITIME)=xin(26)+273.15 Twsucd(ITIME)=xin(27)+273.15 psteam(ITIME)=xin(28)*101325. Pev(ITIME)=xin(29)/3.6 IF (Sum.EQ.0.) THEN Pcd(ITIME)=xin(30)/3.6 Pabs(ITIME)=xin(31)/3.6 Pcdpabs(ITIME)=xin(32)/3.6 ELSE Pcdpabs(ITIME)=xin(32)/3.6 Pcd(ITIME)=Pcdpabs(ITIME)/2. Pabs(ITIME)=Pcdpabs(ITIME)/2. ENDIF ENDIF IF (ITIME.EQ.N) THEN C1*** SECOND PART : identification and printing of the identified C1*** parameters and the error status indicator Fmin=999D25 Iter=0 ErrDetec=0 Nev=IDINT((AUevf-AUevi)/dAUev+1) Ncd=IDINT((AUcdf-AUcdi)/dAUcd+1) Nabs=IDINT((AUabsf-AUabsi)/dAUabs+1) Ngen=IDINT((AUgenf-AUgeni)/dAUgen+1) Nexch=IDINT((AUexchf-AUexchi)/dAUexch+1) Npump=IDINT((Vf-Vi)/dV+0.5)+1 C1*** For each value of the evaporator heat transfer coefficient, C1*** we calculate : AUev=AUevi DO 20 Iev=1,Nev C1*** For each value of the condenser heat transfer coefficient, C1*** we calculate : AUcd=AUcdi DO 30 Icd=1,Ncd C1*** For each value of the absorber heat transfer coefficient, C1*** we calculate : AUabs=AUabsi DO 40 Iabs=1,Nabs C1*** For each value of the generator heat transfer coefficient, C1*** we calculate : AUgen=AUgeni DO 50 Igen=1,Ngen C1*** For each value of the heat-exchanger heat transfer coefficient, C1*** we calculate : AUexch=AUexchi DO 60 Iexch=1,Nexch C1*** For each value of the volume flow rate of the pump of the C1*** compressor, we calculate : VLow=Vi DO 70 Ipump=1,Npump C1*** For each working point, call the subroutine ABSFLSIM in order C1*** to calculate the cooling capacity associated with these C1*** parameters of the chiller F=0. DO 80 I=1,N xin73(1)=Mfrwev(I)*3600. xin73(2)=Choice xin73(3)=Mfrwabs(I)*3600. xin73(4)=Mfrwcd(I)*3600. xin73(5)=Twsuev(I)-273.15 xin73(6)=Twsuabs(I)-273.15 xin73(7)=Twsucd(I)-273.15 xin73(8)=psteam(I)/101325. xin73(9)=Pev(I)*3.6 xin73(10)=Pabs(I)*3.6 xin73(11)=Pcd(I)*3.6 par73(1)=SNGL(AUev*3.6) par73(2)=SNGL(AUcd*3.6) par73(3)=SNGL(AUabs*3.6) par73(4)=SNGL(AUgen*3.6) par73(5)=SNGL(AUexch*3.6) par73(6)=SNGL(VLow*3600.) CALL TYPE73 (TIME,XIN73,OUT73,T,DTDT,PAR73,INFO73,ICNTRL,*75) CALL LINKCK('TYPE71','TYPE73 ',1,99) 75 CONTINUE PevCal(I)=out73(1)/3.6 PabsCal(I)=out73(2)/3.6 PcdCal(I)=out73(3)/3.6 Pgen=out73(4)/3.6 COP=out73(5) Twexev=out73(6)+273.15 Twexabs=out73(7)+273.15 Twexcd=out73(8)+273.15 MfrSteam=out73(9)/3600. XUp=out73(10) h4=out73(11) XLow=out73(12) Err(I)=out73(13) IF ((Err(I).EQ.1).OR.(Err(I).EQ.3)) THEN ErrDetec=Err(I) GOTO 999 ENDIF IF (Sum.EQ.0) THEN F=F+((PevCal(I)-Pev(I))/Pev(I))**2 & +((PcdCal(I)-Pcd(I))/Pcd(I))**2 & +((PabsCal(I)-Pabs(I))/Pabs(I))**2 ELSE F=F+((PevCal(I)-Pev(I))/Pev(I))**2 & +((PcdCal(I)+PabsCal(I)-Pcdpabs(I))/Pcdpabs(I))**2 ENDIF 80 CONTINUE C1*** IF F is lower than the smallest value found so far THEN store C1*** the values of the variables associated with the parameters C1*** considered IF (F.LT.Fmin) THEN C2*** Store the value of the function to minimize Fmin=F C2*** Store the chiller parameters Id(1)=AUev Id(2)=AUcd Id(3)=AUabs Id(4)=AUgen Id(5)=AUexch Id(6)=VLow C2*** For each working point, store the cooling capacity C2*** calculated by means of the parameters considered DO 90 I=1,N Id(6+I)=PevCal(I) Id(6+N+I)=PcdCal(I) Id(6+2*N+I)=PabsCal(I) Id(6+3*N+I)=Err(I) 90 CONTINUE ENDIF C1*** Calculate a new value of the volume flow rate of the C1*** pump of the absorption compressor VLow=VLow+dV 70 CONTINUE C1*** Calculate a new value of the heat-exchanger heat transfer C1*** coefficient AUexch=AUexch+dAUexch 60 CONTINUE C1*** Calculate a new value of the generator heat transfer C1*** coefficient AUgen=AUgen+dAUgen 50 CONTINUE C1*** Calculate a new value of the absorber heat transfer C1*** coefficient AUabs=AUabs+dAUabs 40 CONTINUE C1*** Calculate a new value of the condenser heat transfer C1*** coefficient AUcd=AUcd+dAUcd 30 CONTINUE C1*** Calculate a new value of the evaporator heat transfer C1*** coefficient AUev=AUev+dAUev 20 CONTINUE DO 100 I=1,N IF (Id(6+3*N+I).EQ.2) THEN ErrDetec=2 ENDIF 100 CONTINUE C1*** "Printing" of the identified parameters out(1)=Id(1)*3.6 out(2)=Id(2)*3.6 out(3)=Id(3)*3.6 out(4)=Id(4)*3.6 out(5)=Id(5)*3.6 out(6)=Id(6)*3600. out(7)=Fmin ENDIF IF ((ITIME.GT.N).AND.(ITIME.LE.N2)) THEN C1*** THIRD PART : "printing" of the actual and calculated cooling C1*** capacity, heat rejected in the condenser, heat rejected out C1*** of the absorber (the last two ones are identical if Sum=1) IF (Sum.EQ.0.) THEN out(9)=Id(6+ITIME-N)*3.6 out(10)=Pev(ITIME-N)*3.6 out(11)=Id(6+ITIME)*3.6 out(12)=Pcd(ITIME-N)*3.6 out(13)=Id(6+ITIME+N)*3.6 out(14)=Pabs(ITIME-N)*3.6 out(15)=0. out(16)=0. ELSE out(9)=Id(6+ITIME-N)*3.6 out(10)=Pev(ITIME-N)*3.6 out(11)=0. out(12)=0. out(13)=0. out(14)=0. PcdpabsCal(ITIME-N)=Id(6+ITIME)+Id(6+ITIME+N) out(15)=PcdpabsCal(ITIME-N)*3.6 out(16)=Pcdpabs(ITIME-N)*3.6 ENDIF ENDIF C1*** FOURTH PART: For TIME= 2*N+1 to last simulation time, C1*** everything has been done; so, just wait for the last C1*** simulation time. C2*** Remaining output 1 (the error status indicator) C2*** **************** 999 out(8)=ErrDetec RETURN 1 END SUBROUTINE LINKCK(ENAME1,ENAME2,ILINK,LNKTYP) C*************************************************************************** C THIS SUBROUTINE WAS WRITTEN FOR TRNSYS 14.0 LINK CHECKING - THIS ROUTINE C IS CALLED BY OTHER SUBROUTINES WHEN AN UNLINKED SUBROUTINE HAS BEEN C FOUND. LINKCK IS NEEDED IN ORDER TO AVOID PUTTING COMMON BLOCKS LUNITS C AND CONFIG IN THE TRNSYS TYPES - JWT -- 3/93 C*************************************************************************** COMMON /LUNITS/ LUR,LUW,IFORM,LUK COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, 1 PRWARN COMMON /SIM/TIME0,TFINAL,DELT,IWARN CHARACTER*1 TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK, 1 PRWARN CHARACTER*6 ENAME1,ENAME2 INTEGER ILINK,LNKTYP C ILINK = 1 --> GENERATE AN ERROR MESSAGE AND STOP TRNSYS C ILINK = 2 --> GENERATE A WARNING BUT DON'T STOP TRNSYS C ILINK = 3 --> TRNSYS HAS FOUND AN UNLINKED TYPE - GENERATE AN ERROR AND C STOP THE PROGRAM C ILINK = 4 --> WARN THE USER THAT A ROUTINE REQUIRES AN EXTERNAL FUNCTION C ENAME1 --> CALLING PROGRAM THAT NEEDED THE UNLINKED FILE C ENAME2 --> FILE THAT WAS NOT FOUND BY ENAME1 SUBROUTINE C LNKTYP --> TYPE NUMBER THAT IS UNLINKED IF((LNKCHK.EQ.'Y').OR.(LNKCHK.EQ.'y')) THEN IF(ILINK.EQ.1) THEN WRITE(LUW,20) 104,ENAME1,ENAME2 WRITE(LUW,15) CALL MYSTOP(104) ELSE IF(ILINK.EQ.2) THEN WRITE(LUW,20) 104,ENAME1,ENAME2 IWARN=IWARN+1 ELSE IF(ILINK.EQ.3) THEN WRITE(LUW,25) 105,LNKTYP,LNKTYP WRITE(LUW,15) CALL MYSTOP(105) ELSE IF(ILINK.EQ.4) THEN WRITE(LUW,35) LNKTYP,ENAME1,ENAME2 IWARN=IWARN+1 ELSE IF(ILINK.EQ.5) THEN WRITE(LUW,40) 105,LNKTYP,LNKTYP WRITE(LUW,15) CALL MYSTOP(105) ELSE WRITE(LUW,30) ENAME1 IWARN=IWARN+1 ENDIF ENDIF 15 FORMAT(//2X,47H*** SIMULATION TERMINATED WITH ERROR STATUS ***/) 20 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X,A6, 1' REQUIRES THE FILE "',A6,'" WHICH WAS CALLED BUT NOT LINKED.',/1X 1,'PLEASE LINK IN THE REQUIRED FILE AND RERUN THE SIMULATION.') 25 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X, 1'TYPE ',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.', 1/1X,'LINK TYPE ',I3,' BEFORE RUNNING THIS SIMULATION.') 30 FORMAT(/1X,'*****WARNING*****',/1X,'THE LINKCK SUBROUTINE WAS CALL 1ED WITH AN INVALID OPERAND.',/1X,'THE PROGRAM WHICH CALLED LINKCK 1WITH THE IMPROPER OPERAND WAS ',A6,'.',/1X,'PLEASE MAKE SURE THAT 1THE CALLING PROGRAM IS FIXED OR UNLINKED SUBROUTINES MAY ',/1X,'GO 1 UNNOTICED.') 35 FORMAT(/1X,'*****WARNING*****',/1X,'UNIT ',I2,' ',A6,' REQUIRES TH 1E SUBROUTINE ',A6,/1X,'MAKE SURE THAT THIS SUBROUTINE IS LINKED IN 1 TO AVOID PROBLEMS. IT MAY ALREADY BE LINKED IN.',/) 40 FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X, 1'TYPE',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.', 1/1X,'A DUMMY TYPE SUBROUTINE WAS CALLED IN ITS PLACE. PLEASE LINK' 1,/1X,'TYPE',I3,' BEFORE RUNNING THIS SIMULATION OR TURN OFF THE CH 1ECK'/1X,'FOR UNLINKED SUBROUTINES OPTION IN THE CONFIGURATION FILE 1.') RETURN END