SUBROUTINE TYPE230 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) ! ---------------------------------------------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------------------------------------------- ! This subroutine computes the bulk temperature a surface-water body. The ! water body can be used as a heat source/sink in a closed-loop water- ! source heat pump system. ! ! The bulk temperature is determined using a lumped capacitance approach. ! ! This model assumes that one spool is one flow circuit. ! ! ---------------------------------------------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------------------------------------------- ! Developer: Andrew Chiasson and Dr. Jeffrey Spittler ! School of Mechanical and Aerospace Engineering, ! Oklahoma State University ! ! Date: March-April 1998 ! Edited: April 1999 - added convection coefficient subroutine and changed ! convection correlations to those for a horizontal flat ! given by Incropera & DeWitt (1996) ! - added solar absorption calculations after Duffie and ! Beckman (1991) ! ! ! ---------------------------------------------------------------------------------------------------------------------- ! *** MODEL PARAMETERS *** ! ---------------------------------------------------------------------------------------------------------------------- ! ! PAR(1) = MODE (FOR FUTURE DEVELOPMENTS) ! PAR(2) = TSTART (INITIAL TEMPERATURE OF POND) (C) ! PAR(3) = PONDORIENTATION (LENGTH AZIMUTH IN TERMS OF ! 0 - 90 DEGREES FROM NORTH) (degrees) ! PAR(4) = LENGTH (POND LENGTH) (m) ! PAR(5) = WIDTH (POND WIDTH) (m) ! PAR(6) = DEPTH (POND DEPTH) (m) ! PAR(7) = EMMIS (EMMISSIVITY COEFFICIENT) (--) ! PAR(8) = EXTINCT (EXTINCTION COEFFIFIENT OF WATER) (m^-1) ! PAR(9) = SPOOLS (NUMBER OF COIL SPOOLS) (--) ! PAR(10) = SPOOLLENGTH (LENGTH OF EACH SPOOL) (m) ! PAR(11) = PIPDIA (OUTER PIPE DIAMETER) (m) ! PAR(12) = KPIPE (THERMAL CONDUCTIVITY OF PIPE MATERIAL) (W/(m-C)) ! PAR(13) = WALLT (WALL THICKNESS OF PIPE) (m) ! PAR(14) = FLUID (FLUID TYPE = 1 FOR WATER, ELSE IS GS-4) (--) ! PAR(15) = CONC (% WT. OF GS-4 IN SOLUTION IF USED) (%) ! PAR(16) = KSOIL (THERMAL CONDUCTIVITY OF UNDERLYING SOIL) (W/(m-C)) ! PAR(17) = FF (FOULING FACTOR FOR PIPE SURFACES) ((m^2-C)/W) ! PAR(18) = TGW (CONSTANT TEMPERATURE OF GROUNDWATER) (C) ! PAR(19) = KHSOIL (HYDRAULIC CONDUCTIVITY OF SOIL) (m/s) ! PAR(20) = GRADIENT (HYDRAULIC GRADIENT) (--) ! PAR(21) = DGW (DEPTH TO WATER TABLE) (m) ! ---------------------------------------------------------------------------------------------------------------------- ! *** MODEL INPUTS *** ! ---------------------------------------------------------------------------------------------------------------------- ! XIN(1) = TAMB (AMBIENT AIR TEMPERATURE) (C) ! XIN(2) = WAIR (HUMIDITY RATIO OF AIR) (kg WATER/kg DRY AIR) ! XIN(3) = TSKY (SKY TEMPERATURE) (C) ! XIN(4) = WINDSPEED (m/s) ! XIN(5) = WIND_DIR (WIND DIRECTION) (degrees from north 0-90) ! XIN(6) = SOLARRAD (SOLAR RADIATION ON HORIZONTAL) (KJ/h-m^2)-->(W/m^2) ! XIN(7) = THETA (SOLAR ANGLE OF INCIDENCE) (deg) -->(radians) ! XIN(8) = TIN (INLET FLUID TEMPERATURE) (C) ! XIN(9) = MDOT (MASS FLOW RATE OF FLUID) (kg/hr) ! ---------------------------------------------------------------------------------------------------------------------- ! *** MODEL OUTPUTS *** ! ---------------------------------------------------------------------------------------------------------------------- ! OUT(1) = TPOND (BULK TEMPERATURE OF THE POND) (C) ! OUT(2) = TOUT (OUTLET FLUID TEMPERATURE) (C) ! OUT(3) = MDOT (MASS FLOW RATE OF FLUID) (kg/hr) ! OUT(4) = QTRANSFER (HEAT ABSORBED/REJECTED TO POND) (kJ/hr) ! ---------------------------------------------------------------------------------------------------------------------- ! *** INTERNAL VARIABLES *** ! ---------------------------------------------------------------------------------------------------------------------- ! __AIR = PERTAINING TO AIR ABOVE POND SURFACE ! __botfilm = PERTAINING TO WATER FILM AT POND BOTTOM ! __FILM = PERTAINING TO AIR FILM AT POND SURFACE ! __FL = PERTAINING TO HEAT TRANSFER FLUID ! __GW = PERTAINING TO GROUNDWATER ! __PIPEFILM = PERTAINING TO WATER FILM AT PIPE OUTER SURFACE ! __POND = PERTAINING TO POND WATER ! __WIND = PERTAINING TO WIND ABOVE POND SURFACE ! ALPHA____ = THERMAL DIFFUSIVITY OF ____ (m^2/s) ! AREA = POND AREA (m^2) ! BETA____ = VOLUMETRIC THERMAL EXPANSION COEFFICIENT FOR ____ (K^-1) ! CP____ = HEAT CAPACITY OF ____ (J/(kg-C)) ! DCOEFF = BINARY DIFFUSION COEFFICIENT - WATER TO AIR (m^2/s) ! DELT = TRNSYS VARIABLE INDICATING THE USER-DEFINED TIME STEP (hr) ! GWFLUX = VOLUMETRIC GROUNDWATER FLOW RATE (m^3/s) ! GWMASS = MASS FLOW RATE OF GROUNDWATER (kg/s) ! Hair = FREE CONVECTION HEAT TRANSFER COEFF. AT POND SURFACE (W/m^2-C) ! HC = POND SURFACE CONVECTION HEAT TRANSFER COEFF. (W/m^2-C) ! HD = CONVECTION MASS TRANSFER COEFF. BASED ON J-FACTOR ANALOGY (kg/(m^2-s)) ! Hin = PIPE INTERNAL CONVECTION HEAT TRANSFER COEFFICIENT (W/m^2-C) ! Hout = PIPE EXTERNAL CONVECTION HEAT TRANSFER COEFFICIENT (W/m^2-C) ! HR = LINEARIZED RADIATION COEFFICIENT (W/m^2-C) ! Hwind = FORCED CONVECTION HEAT TRANSFER COEFF. AT POND SURFACE (W/m^2-C) ! K____ = THERMAL CONDUCTIVITY OF ____ (W/(m-C)) ! KHSOIL = HYDRAULIC CONDUCTIVITY OF SOIL (m/s) ! LE = LEWIS NUMBER (--) ! LENGTHCHAR = CHARACTERISTIC POND LENGTH FOR FORCED CONVECTION (m) ! LENGTHCRIT = CRITICAL LENGTH FOR LAMINAR/TURBULENT TRANSITION (m) ! MASSFLUX = MASS FLUX DUE TO EVAPORATION (kg/(m^2-s)) ! NAIR = INDEX OF REFRACTION FOR AIR (--) ! NU____ = NUSSELT NUMBER FOR ____ (--) ! NWATER = INDEX OF REFRACTION FOR WATER (--) ! OLD_T = POND TEMPERATURE AT PREVIOUS ITERATION (C) ! PERIMETER = POND PERIMETER (m) ! PR____ = PRANDTL NUMBER FOR ____ (--) ! QCIRCUIT = CONVECTION HEAT TRANSFER FROM HEAT EXCHANGE FLUID ! PER FLOW CIRCUIT (ONE SPOOL) (W) ! QCONVS = CONVECTION HEAT TRANSFER FROM POND SURFACE (W) ! QEMIS = THERMAL RADIATION HEAT TRANSFER FROM POND SURFACE (W) ! QEVAP = HEAT OF EVAPORATION (MASS TRANSFER) FROM POND SURFACE (W) ! QFLUID = TOTAL HEAT TRANSFER OF ALL SPOOLS (W) ! QGROUND = HEAT TRANSFER TO SURROUNDING GROUND FROM POND (W) ! QGROUNDWATER = HEAT TRANSFER FROM GROUNDWATER INFLUX (W) ! QRAD = SOLAR RADIATION HEAT GAIN BY POND (W) ! RA____ = RAYLEIGH NUMBER FOR ____ (--) ! RE____ = REYNOLD'S NUMBER FOR ____ (--) ! REFLECTANCE= REFLECTANCE OF THE POND WATER SURFACE (--) ! REFRACT = REFRACTION ANGLE - AIR TO WATER (radians) ! RHO____ = DENSITY OF ____ (kg/m^3) ! ri = INNER PIPE RADIUS (m) ! Rin = PIPE INTERNAL RESISTANCE (m^2-C/W) ! ro = OUTER PIPE RADIUS (m) ! Rout = PIPE EXTERNAL RESISTANCE (m^2-C/W) ! RPARA = PARALLEL COMPONENT OF UNPOLARIZED RADIATION (--) ! RPERP = PERPENDICULAR COMPONENT OF UNPOLARIZED RADIATION (--) ! RPIPE = PIPE WALL RESISTANCE (m^2-C/W) ! RTOTAL = TOTAL PIPE RESISTANCE (m^2-C/W) ! TAOa = ABSORPTANCE OF SOLAR RADIATION BY THE POND WATER (--) ! TAOr = TRANSMITTANCE OF UNPOLARIZED RADIATION (--) ! TF = FINAL POND TEMPERATURE AT PREVIOUS TIME STEP (C) ! TFILM = AVERAGE OF POND AND AIR TEMPERATURES (C) ! TFLUID = AVERAGE HEAT EXCHANGE FLUID TEMPERATURE (C) ! TI = INITIAL POND TEMPERATURE AT NEW TIME STEP (C) ! TIME = TRNSYS VARIABLE INDICATING CURRENT TIME OF SIMULATION (hr) ! TPIPE = FILM TEMPERATURE AT PIPE SURFACE (C) ! UAPIPE = OVERALL HEAT TRANSFER COEFFICIENT FOR PIPE (W/C) ! WSURF = HUMIDITY RATIO OF SATURATED AIR AT POND SURFACE (--) ! X = EMPIRICAL EXPONENT FOR Pr IN DITTUS-BOELTER EQUATION (--) ! '1' in Q terms denotes dependence on TPOND for use in DIFFEQ subroutine ! '2' in Q terms denotes independence of TPOND for use in DIFFEQ subroutine ! ---------------------------------------------------------------------------------------------------------------------- ! *** SUBROUTINES CALLED *** ! ---------------------------------------------------------------------------------------------------------------------- ! ! (1) WATER_PROPERTIES ! (2) GS4_PROPERTIES ! (3) AIR_PROPERTIES ! ---------------------------------------------------------------------------------------------------------------------- ! *** FUNCTIONS CALLED *** ! ---------------------------------------------------------------------------------------------------------------------- ! ! (1) HFG ! ! ---------------------------------------------------------------------------------------------------------------------- ! Use of the Storage array ! ---------------------------------------------------------------------------------------------------------------------- ! ! stored(1) = T_pond during previous time step - only updated at convergence call, at the end of time step ! stored(2) = T_pond during previous iteration - updated at each iterative call ! ! ---------------------------------------------------------------------------------------------------------------------- ! ! !************************************************************************ !dec$attributes dllexport :: type230 use TrnsysConstants use TrnsysFunctions implicit none ! --- TRNSYS declarations ---------------------------------------------------------------------------------------------- integer, parameter :: NI=9, NP =21, ND=0, NO=10, NS = 2 real(8), intent(in) :: time, xin, par, t real(8), intent(inout) :: dtdt real(8), intent(out) :: out integer, intent(inout) :: info, iCntrl dimension xin(NI), par(NP), OUT(NO), info(15) character*3 OCHECK(NO) !AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE OUTPUTS character*3 YCHECK(NI) !AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE INPUTS ! --- Local variables -------------------------------------------------------------------------------------------------- ! --- Parameters: integer MODE real(8) TSTART real(8) PONDORIENTATION real(8) LENGTH real(8) WIDTH real(8) DEPTH real(8) EMMIS real(8) EXTINCT real(8) SPOOLS real(8) SPOOLLENGTH real(8) PIPDIA real(8) KPIPE real(8) WALLT real(8) FLUID real(8) CONC real(8) KSOIL real(8) FF real(8) TGW real(8) KHSOIL real(8) GRADIENT real(8) DGW ! --- Outputs: real(8) TPOND real(8) TOUT real(8) QTRANSFER real(8) Rin real(8) Rpipe real(8) Rout real(8) QEVAP2 real(8) qground real(8) Refl ! --- Inputs real(8) TAMB real(8) WAIR real(8) TSKY real(8) WINDSPEED real(8) WIND_DIR real(8) SOLARRAD real(8) THETA real(8) TIN real(8) MDOT ! --- Intermediate variables real(8) stored(NS) real(8) PSYDAT(9),STATUS real(8) pi, GRAVITY,NAIR,NWATER character*1 TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHECK,PRWARN real(8) OLD_T, TF,TI, time0 real(8) AREA, PERIMETER,RO,RI,PONDVOL,MDOTT,QCIRCUIT1,QCIRCUIT2 real(8) QFLUID1,QFLUID2,RHOPOND,CPPOND,KPOND,VISCPOND,BETAPOND real(8) RHOGW,CPGW,KGW,VISCGW,BETAGW real(8) TFLUID,RHOFL,CPFL,KFL,VISCFL,BETAFL,NUFL,PRFL real(8) TFILM,RHOFILM,CPFILM,KFILM,VISCFILM,BETAFILM,ALPHAFILM real(8) HIN,TPIPE,X,HOUT,RTOTAL,UAPIPE,REFRACT,RPERP,RPARA,TAOR,TAOA real(8) RHOPIPEFILM,CPPIPEFILM,KPIPEFILM,VISCPIPEFILM,BETAPIPE,ALPHAPIPEFILM,PRPIPEFILM real(8) RAPIPE,RAPIPEFILM,NUPIPEFILM,RAFILM real(8) RAAIR,HAIR,NUAIR,WIND_DIR_R,PONDORIENTATION_R,LENGTHCHAR,REWINDY,PRWIND,LENGTHCRIT,NUWIND real(8) HWIND,HC,QCONVS1,QCONVS2,UBG,QGROUND1,QGROUND2,GWFLUX,GWMASS,QGROUNDWATER1,QGROUNDWATER2 real(8) WSURF,DCOEFF,LE,HD,MASSFLUX,HFG,AA,BB,TBAR,QCIRCUIT,TOUTNEW real(8) REFLECTANCE,QRAD2,HR,QEMIS1,QEMIS2 integer TCHECK data GRAVITY/9.81/ data NAIR/1.00/ data NWATER/1.33/ pi= 4.0*datan(1.d0) !----------------------------------------------------------------------------------------------------------------------- ! GET THE VALUES OF THE PARAMETERS FOR THIS COMPONENT MODE = INT(PAR(1) + 0.1) TSTART = PAR(2) PONDORIENTATION = PAR(3) LENGTH = PAR(4) WIDTH = PAR(5) DEPTH = PAR(6) EMMIS = PAR(7) EXTINCT = PAR(8) SPOOLS = PAR(9) SPOOLLENGTH = PAR(10) PIPDIA = PAR(11) KPIPE = PAR(12) WALLT = PAR(13) FLUID = PAR(14) CONC = PAR(15) KSOIL = PAR(16) FF = PAR(17) TGW = PAR(18) KHSOIL = PAR(19) GRADIENT = PAR(20) DGW = PAR(21) !----------------------------------------------------------------------------------------------------------------------- ! --- Initial call to detect the TRNSYS version for which this Type is written ----------------------------------------- if (info(7) .eq. -2) then info(12) = 16 ! This component is a TRNSYS 16 Type return 1 endif !----------------------------------------------------------------------------------------------------------------------- ! --- Very last call in simulation ------------------------------------------------------------------------------------- if (info(8).eq.-1) then return 1 ! endif !----------------------------------------------------------------------------------------------------------------------- ! --- Extra call after convergence at a time step - The type should not modify its outputs here ------------------------ if (info(13) .gt. 0) then !Update pond temperature with converged value call getStorageVars(stored,NS,info) stored(1) = stored(2) call setStorageVars(stored,NS,info) return 1 endif !----------------------------------------------------------------------------------------------------------------------- ! --- Initialization call (not a simulation call) ---------------------------------------------------------------------- if (info(7) .eq. -1) then info(6) = NO ! Set outputs number info(9) = 1 ! This type's outputs depend upon the passage of time info(10) = 0 call typeck(1,info,NI,NP,ND) ! Check the number of inputs/parameters/derivatives call setStorageSize(NS,info) !Set number of storage spots return 1 endif !----------------------------------------------------------------------------------------------------------------------- ! --- First time step - No iterations ---------------------------------------------------------------------------------- if ( time .lt. (getSimulationStartTime()+getSimulationTimeStep()/2.0) ) then !Check values of parameters if (MODE /= 1) call typeck(-4,info,0,1,0) ! if (TSTAR < 20) call typeck(-4,info,0,2,0) ! if (PONDORIENTATION < 0) call typeck(-4,info,0,3,0) if (LENGTH <= 0.D0) call typeck(-4,info,0,4,0) if (WIDTH <= 0.D0) call typeck(-4,info,0,5,0) if (DEPTH <= 0.D0) call typeck(-4,info,0,6,0) if (EMMIS <= 0.D0) call typeck(-4,info,0,7,0) if (EXTINCT <= 0.D0) call typeck(-4,info,0,8,0) if (SPOOLS <= 0.D0) call typeck(-4,info,0,9,0) if (SPOOLLENGTH <= 0.D0) call typeck(-4,info,0,10,0) if (PIPDIA <= 0.D0) call typeck(-4,info,0,11,0) if (KPIPE < 0.D0) call typeck(-4,info,0,12,0) if (WALLT <= 0.D0) call typeck(-4,info,0,13,0) ! if (FLUID <= 0.D0) call typeck(-4,info,0,14,0) if (CONC < 0.D0) call typeck(-4,info,0,15,0) if (KSOIL < 0.D0) call typeck(-4,info,0,16,0) if (FF < 0.D0) call typeck(-4,info,0,17,0) ! if (TGW <= 0.D0) call typeck(-4,info,0,18,0) if (KHSOIL < 0.D0) call typeck(-4,info,0,19,0) ! if (GRADIENT <= 0.D0) call typeck(-4,info,0,20,0) if (DGW <= 0.D0) call typeck(-4,info,0,21,0) STORED(1) = TSTART STORED(2) = TSTART call setStorageVars(stored,NS,info) !Initialize outputs TF=TSTART TOUT=TIN OUT(1)=XIN(1) !TPOND OUT(2)=XIN(8) !TOUT OUT(3)=XIN(9) !MDOT OUT(4)=0.D0 !QTRANSFER out(5)=Rin out(6)=Rpipe out(7)=Rout out(8)=QEVAP2 out(9)=qground out(10)=Refl return 1 endif !----------------------------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------------------------- ! GET THE VALUES OF THE INPUTS FOR TRNSYS TIME STEP TAMB = XIN(1) WAIR = XIN(2) TSKY = XIN(3) WINDSPEED = XIN(4) WIND_DIR = XIN(5) SOLARRAD = XIN(6)*3.6 !Convert from W/m^2 to kJ/h-m^2 THETA = XIN(7)*pi/180.0 !Convert from degrees to radians TIN = XIN(8) MDOT = XIN(9) ! Get stored values ---------------------------------------------------------------------------------------------------- call getStorageVars(stored,NS,info) OLD_T=stored(1) MDOTT=MDOT/3600./SPOOLS !convert mass flow rate to kg/sec/circuit QCIRCUIT1=0.0 QCIRCUIT2=0.0 QFLUID1=0.0 QFLUID2=0.0 !compute some constants AREA = LENGTH*WIDTH PERIMETER = 2.*LENGTH+2.*WIDTH ro = PIPDIA/2. ri = ro-WALLT pondvol = AREA*DEPTH ! SET TEMPERATURE OF POND EQUAL TO AMBIENT AIR TEMPERATURE ! ON FIRST PASS OF SIMULATION IF (INFO(7).EQ.0) THEN TI=TF !set initial temp for diffeq on 1st iteration TPOND=TF !set pond initial temp to previous final temp END IF ! START OF TEMPERATURE ITERATION LOOP DO TCHECK=1,1000 ! GET TEMPERATURE-DEPENDENT PROPERTIES OF POND WATER CALL WATER_PROPERTIES (TPOND,RHOPOND,CPPOND,KPOND,VISCPOND,BETAPOND) ! COMPUTE TEMPERATURE-DEPENDENT PROPERTIES OF GROUNDWATER !for groundwater heat transport calculations CALL WATER_PROPERTIES (TGW,RHOGW,CPGW,KGW,VISCGW,BETAGW) IF (MDOT.GT.0) THEN ! COMPUTE TEMPERATURE-DEPENDENT PROPERTIES OF HEAT EXCHANGE FLUID !*FOR WATER* TFLUID=(TOUT+TIN)/2. IF (FLUID.EQ.1) THEN CALL WATER_PROPERTIES(TFLUID,RHOFL,CPFL,KFL,VISCFL,BETAFL) ELSE !*FOR GS-4 ANTIFREEZE SOLUTION* CALL GS4_PROPERTIES(TFLUID,CONC,RHOFL,CPFL,KFL,VISCFL) END IF ! *** COMPUTE CONVECTION HEAT FLUX FROM HEAT EXCHANGE FLUID *** !Compute internal pipe resistance REfl=4*MDOTT/(pi*(2*ri)*VISCFL) !Check flow conditions IF (REfl.LE.2000.) THEN !flow is laminar NUfl=4.36 ELSE !flow is turbulent PRfl=VISCfl*CPfl/Kfl ! Check for heating or cooling case IF (TIN.GT.TPOND) THEN X=0.3 ELSE X=0.4 END IF !Compute Nu from Dittus-Boelter Relation NUfl=0.023*REfl**0.8*PRfl**X ! !Compute Nu from Petukhov Relation ! fff=(0.79*log(REfl)-1.34)**(-2.) ! NUfl=REfl*PRfl*fff/8. ! . /( 1.07+12.7* (fff/8.)**0.5* (PRfl**(2./3.)-1.) ) END IF Hin=NUfl*Kfl/(2.*ri) RIN=1./hin !Compute pipe wall resistance RPIPE=ri/KPIPE*LOG(ro/ri) !Compute external pipe resistance TPIPE=(TFLUID+TPOND)/2. !approximation of film temperature CALL WATER_PROPERTIES (TPIPE,RHOPIPEFILM,CPPIPEFILM,KPIPEFILM,VISCPIPEFILM,BETAPIPE) ALPHAPIPEFILM=KPIPEFILM/(RHOPIPEFILM*CPPIPEFILM) PRPIPEFILM=VISCPIPEFILM*CPPIPEFILM/KPIPEFILM RAPIPE=GRAVITY*BETAPIPE*(TFLUID-TPOND)*PIPDIA**3./(VISCPIPEFILM/RHOPIPEFILM*ALPHAPIPEFILM) RAPIPEFILM=ABS(RAPIPE) NUPIPEFILM=(0.6+0.387*RAPIPEFILM**(1./6.)/(1+(0.559/PRPIPEFILM)**(9./16.))**(8./27.))**2. hout=NUPIPEFILM*KPIPEFILM/PIPDIA ROUT=(ri/ro)*(1./hout) RTOTAL=RIN+RPIPE+ROUT+FF UAPIPE=2.*pi*ri*SPOOLLENGTH/RTOTAL !in terms of inside pipe area QCIRCUIT1=UAPIPE*(-1.) QCIRCUIT2=UAPIPE*TFLUID QFLUID1=QCIRCUIT1*SPOOLS !total heat transfer QFLUID2=QCIRCUIT2*SPOOLS !total heat transfer END IF ! *** COMPUTE SOLAR RADIATION HEAT GAIN *** REFRACT=ASIN(SIN(THETA)*NAIR/NWATER)+0.001 RPERP=(SIN(REFRACT-THETA))**2./(SIN(REFRACT+THETA))**2. RPARA=(TAN(REFRACT-THETA))**2./(TAN(REFRACT+THETA))**2. TAOr=0.5*((1.-RPERP)/(1.+RPERP)+(1.-RPARA)/(1.+RPARA)) TAOa=EXP(-1.*EXTINCT*DEPTH/COS(REFRACT)) reflectance=TAOa-TAOa*TAOr QRAD2=SOLARRAD*(1.-reflectance)*AREA ! *** COMPUTE THERMAL RADIATION HEAT TRANSFER *** HR=4.*EMMIS*5.67E-8*(((TPOND+273.15)+(TSKY+273.15))/2.)**3 QEMIS1=AREA*HR*(-1.) QEMIS2=AREA*HR*TSKY ! *** COMPUTE CONVECTION HEAT TRANSFER *** !Compute temperature_dependent properties of air TFILM=(TPOND+TAMB)/2. !the air film temperature CALL AIR_PROPERTIES(TFILM,RHOFILM,CPFILM,KFILM,VISCFILM,BETAFILM) !*FREE CONVECTION CONTRIBUTION ALPHAFILM=KFILM/(RHOFILM*CPFILM) RAair=GRAVITY*BETAFILM*(TPOND-TAMB)*(AREA/PERIMETER)**3./(VISCFILM/RHOFILM*ALPHAFILM) RAfilm=ABS(RAair) IF (RAfilm.LT.1.0E+7) THEN !check validity range NUair=0.54*RAfilm**(0.25)+0.0001 ELSE NUair=0.15*RAfilm**(1./3.)+0.0001 END IF Hair=NUair*KFILM/(AREA/PERIMETER) !*FORCED CONVECTION CONTRIBUTION !first determine characteristic length IF (WIND_DIR.GE.180.) THEN !adjust all wind directions to 0 to 180 WIND_DIR=WIND_DIR-180. IF (WIND_DIR.GE.90.) THEN !adjust all wind directions to 0 to 90 WIND_DIR=180.-WIND_DIR END IF END IF WIND_DIR_R=WIND_DIR*pi/180. !convert to radians PONDORIENTATION_R=PONDORIENTATION*pi/180. !convert to radians LENGTHCHAR=LENGTH*MAX( SIN(PONDORIENTATION_R), COS(WIND_DIR_R) )+WIDTH*MIN( SIN(WIND_DIR_R),COS(PONDORIENTATION_R) ) REwindy=RHOFILM*WINDSPEED*LENGTHCHAR/VISCFILM PRwind=VISCFILM*CPFILM/KFILM !find the critical length assuming RE=5x10^5 is critical LENGTHCRIT=VISCFILM*(5.E+5)/(RHOFILM*WINDSPEED+0.001) IF (LENGTHCRIT.GE.LENGTHCHAR) THEN !all flow is laminar NUwind=0.664*REwindy**0.5*PRwind**(1./3.) ELSE !flow is mixed NUwind=0.037*REwindy**0.8*PRwind**(1./3.) END IF Hwind=NUwind*KFILM/LENGTHCHAR !! HC=Hwind+Hair HC=MAX(Hwind,Hair) !-->> HC=5.7+3.8*WINDSPEED !old convection model from Kishore and Joshi (1984) QCONVS1=AREA*HC*(-1.) QCONVS2=AREA*HC*TAMB ! *** COMPUTE GROUND CONDUCTION HEAT TRANSFER *** !ground conduction model from Hull (1984) UBg=0.999*KSOIL/(DGW-DEPTH)+1.37*KSOIL*PERIMETER/AREA QGROUND1=UBg*AREA*(-1.) QGROUND2=UBg*AREA*TGW QGROUND=QGROUND1*TPOND+QGROUND2 ! *** COMPUTE HEAT TRANSFER FROM GROUNDWATER INFLUX *** IF (DGW.LT.DEPTH) THEN GWFLUX=KHSOIL*GRADIENT*(PERIMETER*(DEPTH-DGW)+AREA) !Darcy's Law GWMASS=GWFLUX*RHOGW !mass flow QGROUNDWATER1=GWMASS*CPGW*(-1.) QGROUNDWATER2=GWMASS*CPGW*TGW END IF ! *** COMPUTE EVAPORATION HEAT/MASS TRANSFER *** !Get humidity ratio of air at pond surface from TRNSYS psych routine PSYDAT(1)=1.0 !pressure in atmospheres PSYDAT(2)=TPOND !Tdb for TRNSYS psychrometric routine PSYDAT(3)=TPOND !Twb for TRNSYS psychrometric routine CALL PSYCHROMETRICS(TIME,INFO,1,1,0,PSYDAT,1,STATUS,*200) 200 CONTINUE WSURF=PSYDAT(6) !humidity ratio (kg WATER/kg DRY AIR) DCOEFF=1.87E-10*(TFILM+273.15)**2.072 LE=ALPHAFILM/DCOEFF ! HD=HC/(1005.*LE**0.66666667) !1005=cp of dry air in J/(kgC) HD=HC/(CPFILM*LE**(2./3.) ) MASSFLUX=HD*(WAIR-WSURF) QEVAP2=AREA*MASSFLUX*HFG(TPOND) ! *** COMPUTE BULK TEMPERATURE OF POND *** !using TRNSYS analytical differential equation solver AA=(QFLUID1+QEMIS1+QCONVS1+QGROUND1+QGROUNDWATER1) *3600./(PONDVOL*RHOPOND*CPPOND) BB=(QRAD2+QFLUID2+QEMIS2+QCONVS2+QGROUND2+QGROUNDWATER2 +QEVAP2)*3600./(PONDVOL*RHOPOND*CPPOND) CALL DIFFERENTIAL_EQN(TIME,AA,BB,TI,TF,TBAR) TPOND=TBAR ! *** COMPUTE RETURN FLUID TEMPERATURE *** IF (MDOT.GT.0) THEN !COMPUTE OUTLET FLUID TEMPERATURE QCIRCUIT=QCIRCUIT2+QCIRCUIT1*TPOND TOUTnew=TFLUID-QCIRCUIT/(2.*MDOTT*CPFL) TOUT=TOUTnew + 0.75*(TOUT-TOUTnew) !relax TOUT END IF !Check for convergence IF (ABS(OLD_T-TPOND).GT.0.00001) GOTO 75 GOTO 90 75 OLD_T=TPOND ENDDO !-------------------------------------------- ! END OF TEMPERATURE ITERATION LOOP ! COMPUTE TOTAL AMOUNT OF HEAT EXCHANGED BY FLUID IN kJ/hr 90 IF (MDOT.GT.0) THEN QTRANSFER=MDOT*(TIN-TOUT)*CPFL/1000. ELSE QTRANSFER=0.0 END IF ! SET THE OUTPUTS OUT(1)=TPOND OUT(2)=TOUT OUT(3)=MDOT OUT(4)=QTRANSFER out(5)=Rin out(6)=Rpipe out(7)=Rout out(8)=QEVAP2 out(9)=qground out(10)=Refl !Update pond temperature during this iteration stored(2) = TF call setStorageVars(stored,NS,info) RETURN 1 END ! ************************************************************** SUBROUTINE WATER_PROPERTIES (Tw,RHOw,CPw,Kw,VISCw,BETAw) ! ************************************************************** ! PURPOSE: To compute thermal properties of water ! INPUTS: ! (1) Twater = WATER TEMPERATURE (C) ! OUTPUTS: ! (1) RHOw = DENSITY (kg/m^3) ! (2) CPw = HEAT CAPACITY (J/(kg-C)) ! (3) Kw = THERMAL CONDUCTIVITY (W/(m-C)) ! (4) VISCw = DYNAMIC VISCOSITY (Pa-s) ! (5) BETAw = VOLUMETRIC THERMAL EXPANSION (K^-1) ! COEFFICIENT ! References: ! ! (1) CRC Handbook of Chemistry and Physics, 61st Ed. (1980-1981). ! (2) Irvine, T.F.Jr., and Liley, P.E., "Steam and Gas Tables ! with Computer Equations," Academic Press, Inc., 1984. ! !*********************************************************************** IMPLICIT REAL(8) (A-Z) ! ** Density ** DATA AR0,AR1,AR2,AR6/999.83952,16.945176,-.0079870401,.01687985/ DATA AR3,AR4,AR5/-46.170461E-06,105.56302E-09,-280.54253E-12/ RHOw=(AR0+TW*(AR1+TW*(AR2+TW*(AR3+TW*(AR4+TW*AR5)))))/(1.+AR6*TW) ! ** Dynamic Viscosity ** 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)) VISCw=0.001*WMU ! ** Thermal Conductivity ** DATA AK0,AK1,AK2/0.560101,0.00211703,-1.05172E-05/ DATA AK3,AK4/1.497323E-08,-1.48553E-11/ Kw=AK0+TW*(AK1+TW*(AK2+TW*(AK3+TW*AK4))) ! ** Specific Heat Capacity ** DATA ACP0,ACP1,ACP2/4.21534,-0.00287819,7.4729E-05/ DATA ACP3,ACP4/-7.79624E-07,3.220424E-09/ WCP=ACP0+TW*(ACP1+TW*(ACP2+TW*(ACP3+TW*ACP4))) CPw=WCP*1000. ! ** Volumetric Thermal Expansion Coeff. ** !from curve-fit to data in Table A.6 of Incropera and DeWitt, 3rd ed. BETAw=1.060228571E-5*TW+3.9560E-5 RETURN END ! ************************************************************** SUBROUTINE GS4_PROPERTIES (Tgs,CONC,RHOgs,CPgs,Kgs,VISCgs) ! ************************************************************** ! PURPOSE: To compute thermal properties of GS4 antifreeze ! INPUTS: ! (1) Tgs = GS4 TEMPERATURE (C) ! (2) CONC = % WT. OF GS-4 IN SOLUTION (%) ! OUTPUTS: ! (1) RHOgs = DENSITY (kg/m^3) ! (2) CPgs = HEAT CAPACITY (J/(kg-C)) ! (3) Kgs = THERMAL CONDUCTIVITY (W/(m-C)) ! (4) VISCgs = DYNAMIC VISCOSITY (Pa-s) ! ************************************************************** IMPLICIT REAL (A-Z) !from Wadivkar (1997) Master's Thesis RHOgs=(-0.000226*Tgs+0.00606*0.506*CONC+0.99)*1000. CPgs=(-6.77E-3*0.506*CONC+1.4E-4*Tgs+1.01)*4186.8 Kgs=(0.000717*Tgs+0.00085*0.506*CONC+1.8/0.506/CONC+0.107)*1.73 VISCgs=EXP(-0.135*(Tgs+273.) + 0.000262*(Tgs+273.)**2. + 0.172*0.506*CONC - 0.000564*(Tgs+273.)*0.506*CONC + 17.5)*0.001 RETURN END ! ************************************************************** SUBROUTINE AIR_PROPERTIES (Tair,RHOa,CPa,Ka,VISCa,BETAa) ! ************************************************************** ! PURPOSE: To compute thermal properties of air ! INPUTS: ! (1) Tair = AIR TEMPERATURE (C) ! OUTPUTS: ! (1) RHOa = DENSITY (kg/m^3) ! (2) CPa = HEAT CAPACITY (J/(kg-C)) ! (3) Ka = THERMAL CONDUCTIVITY (W/(m-C)) ! (4) VISCa = DYNAMIC VISCOSITY (Pa-s) ! (5) BETAa = VOLUMETRIC THERMAL EXPANSION (K^-1) ! COEFFICIENT ! Ref: Irvine, T.F.Jr., and Liley, P.E., "Steam and Gas Tables ! with Computer Equations," Academic Press, Inc., 1984. !*********************************************************************** IMPLICIT REAL(8) (A-Z) ! ** Density ** !from curve-fit to data in Table B-4b of McQuiston and Parker, 4th ed. RHOa=-4.9714E-3*Tair+1.33342 ! ** Dynamic Viscosity ** 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=Tair+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 ! ** Thermal Conductivity ** 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/ T=Tair+TCONV Ka=C0+T*(C1+T*(C2+T*(C3+T*(C4+T*C5)))) ! ** Specific Heat Capacity ** DATA D0,D1,D2/1.03409,-0.284887E-3,0.7816818E-6/ DATA D3,D4/-0.4970786E-9,0.1077024E-12/ T=Tair+TCONV CPa=D0+T*(D1+T*(D2+T*(D3+T*D4))) CPa=CPa*1000. ! ** Volumetric Thermal Expansion Coeff. ** BETAa=1./(Tair+273.15) RETURN END !*********************************************************************** FUNCTION HFG(TC) ! ! PURPOSE: To computr the latent heat of vaporization of water (J/kg) ! given sat. T (C) ! ! References: ! 1) Irvine, T.F.Jr., and Liley, P.E., "Steam and Gas Tables ! with Computer Equations," Academic Press, Inc., 1984. !*********************************************************************** IMPLICIT REAL (A-Z) 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*1000. RETURN END