[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [TRNSYS-users] Strange results from Type52
Hello, David.
Thanks for your reply. I have made a slight modification of the model
code in order to trace the problem. Then I find that when the control
signal drops below 0.4, the water flow inside the tube lies in the
transition region. The resulting heat transfer coefficient on the tube
inside surface decreases by more than 80% from the "just turbulent"
situation to the "just laminar" condition. As a linear interpolation is
used in the model to calculate the heat transfer coefficient within the
transition region, the reduction in the coil capacity becomes much
larger than the reduction in the water flow. Hence, the outlet water
temperature starts to decrease when the control signal drops below 0.4.
The simplest way is to rectify is to choose a smaller tube size. It
helps increase the Reynolds number inside the tubes. Another problem is
the adoption of a constant Nusselt number when the flow becomes laminar
inside the tubes in the model code. I have made some modifications of
the original model code. Here, the calculation of the Nusselt number in
the laminar region is based on the formulations from "Nellis G, Klein S.
Heat Transfer. New York: Cambridge University Press, 2009". For the
Nusselt number in the turbulent region, the more accurate equations from
"Gnielinski V. New equations for heat and mass transfer in turbulent
pipe and channel flow. International Chemical Engineering
1976;16:359-68" is followed. These formulations are used in EES in
calculating the heat transfer coefficient for internal pipe flow.
Attached is my revised source code for your reference. The original
equations are converted to comment statements so that you can easily
revert to the original model. Try and see if it is better.
Best Regards
Chun Kwong LEE
City University of Hong Kong
David BRADLEY wrote:
Chun Kwong,
I suspect that when the water flowrate drops below a certain point,
the heat transfer rate between the air and the water begins to
decrease. I modified your simulation so that there is an equation that
sets the pump control signal starting at a value of 1 and ending at 0
after 12 hours. That way, you can see the progression of the outlet
water temperature (I also used a shorter timestep). As you say, the
water outlet temperature starts out increasing but when the flow rate
through your coils gets to about 4000 kg/h (1.1 L/s) the outlet
temperature starts to go back down. At that point, the face velocity
of water through a single tube in the coil is only about 0.2 m/s,
which seems a bit slow. I don't think that you have reached a
laminar/turbulent flow point because there is not a sharp break in the
water outlet temperatures (which you usually see because the Reynolds
number drops very suddenly).
In any case, I am not sure that the model is wrong and I certainly
don't see anything incorrect in your implementation.
Best,
David
On 7/10/2013 02:30, a8304506@graduate.hku.hk wrote:
Hello.
I try to investigate the performance of a cooling coil at different chilled
water flow rates (see attached .tpf file). When I adjust the control signal of
the water pump from 1 down to around 0.4, the chilled water leaving
temperatures increase as expected. However, when I further reduce the control
signal, it appears that the leaving chilled water temperatures begin to
decrease. Is it the problem from the cooling coil model or any other issue?
Please advise. I am using TRNSYS16.1. Thank you!
Best Regards
Chun Kwong LEE
City University of Hong Kong
_______________________________________________
TRNSYS-users mailing list
TRNSYS-users@cae.wisc.edu
https://mailman.cae.wisc.edu/listinfo/trnsys-users
--
***************************
David BRADLEY
Principal
Thermal Energy Systems Specialists, LLC
22 North Carroll Street - suite 370
Madison, WI 53703 USA
P:+1.608.274.2577
F:+1.608.278.1475
d.bradley@tess-inc.com
http://www.tess-inc.com
http://www.trnsys.com
SUBROUTINE TYPE52(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
C-----------------------------------------------------------------------------------------------------------------------
C THIS ROUTINE DETERMINES THE THERMAL PERFORMANCE OF A COOLING COIL *
C USING AN EFFECTIVENESS MODEL AS DESCRIBED IN "METHODOLOGIES FOR *
C THE DESIGN AND CONTROL OF CHILLED WATER SYSTEMS", PHD THESIS, *
C SOLAR ENERGY LABORATORY, UNIVERSITY OF WISCONSIN, J.E. BRAUN, 1988 *
C**************************** PARAMETERS ***************************
C
C MODE = SIMPLE OR DETAILED CALCULATIONS
C IU = UNITS: 1-SI, 2-ENGLISH
C NROWS = NUMBER OF ROWS OF TUBES
C NTUBES = NUMBER TUBES PER ROW
C LT = LENGTH OF TUBES IN EACH ROW
C WD = WIDTH OF DUCT PERPENDICULAR TO THE TUBES
C DO = TUBE OUTER DIAMETER
C DI = TUBE INNER DIAMETER
C KT = TUBE THERMAL CONDUCTIVITY
C FT = FIN THICKNESS
C FS = FIN SPACING
C NFINS = NUMBER OF FINS ON ONE TUBE
C KF = FIN THERMAL CONDUCTIVITY
C FMODE = CONTINUOUS FLAT PLATE OR ANNULAR FINS
C ADIST = DISTANCE BETWEEN TUBE CENTERS ACROSS ROW
C DE = DIAMETER OF ANNULAR FIN
C CDIST = DISTANCE BETWEEN TUBE ROW CENTERLINES (FLOW DIRECTION)
C-----------------------------------------------------------------------------------------------------------------------
! Copyright © 2004 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved.
C-----------------------------------------------------------------------------------------------------------------------
C USE STATEMENTS
USE TrnsysFunctions
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS
!DEC$ATTRIBUTES DLLEXPORT :: TYPE52 !SET THE CORRECT TYPE NUMBER HERE
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C TRNSYS DECLARATIONS
IMPLICIT NONE
DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,TIME0,TFINAL,DELT
INTEGER*4 INFO(15),NP,NI,NOUT,ND,NPAR,IUNIT,ITYPE,ICNTRL,
1 NSTORED
CHARACTER*3 YCHECK,OCHECK
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C USER DECLARATIONS - SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI),
C OUTPUTS (NOUT), AND DERIVATIVES (ND) THAT MAY BE SUPPLIED FOR THIS TYPE
PARAMETER (NP=16,NI=5,NOUT=12,ND=0,NSTORED=0)
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C REQUIRED TRNSYS DIMENSIONS
DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT),
1 T(ND),DTDT(ND)
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES
DOUBLE PRECISION MA,MW,MHUA,KW,MHUW,LT,KT,KF,JF,MCONV,NTUD,NTUW,
& NTUO,NTUI,MSTAR,LCOIL,LIMITW,LIMITD,RNUM,TCONV,A1,A0,HCONV,H,B0,
& B1,CTOF,C1,C0,TS,HS,TSAT,PRW,CPA,CPW,CPV,HFG,PRA,PATM,TOL,PI,HSAT
& ,WD,DO,DI,FT,FS,ADIST,DE,CDIST,PSYDAT,RI,RO,RE,AC,FH,AX,AI,ATO,AF
& ,AO,FRATIO,WA1,TW1,UM,TA1,TA2,WA2,TW2,QCOIL,QSENS,QLAT,RHA2,RHA1,
& QDAIR,FDRY,GW,TPROP,REI,HI,HIT,HIL,DH,CJ1,CJ2,GA,REO,HDO,ALPHA,
& BETA,EFF,FIN,EFFD,CPM,RA,TDP,HA1,HW1,CSTAR,TEMP,C,EPSD,TA2D,TS2D,
& TW2W,HW2,CS,HWO,EFFW,ERAT,EPSW,G,G1,TOLD,G2,DGDT,TS1W,CK,TW2D,
& TMAX,TMIN,GMAX,TWX,GMIN,SGNMIN,DMAX1,DMIN1,HA2,HSAVG,TSAVG,TAX,
& HAX2,HAX,HSW,TSW,
& GZ,NU_LAM,RE_TUR,FD,NU_TUR !Added by Avatar Lee on 17/07/2013
INTEGER MODE,FMODE,EMODE,ROUND,IU,IMAX,NROWS,NTUBES,NFINS,IWARN,
& ISTAT,ITER
LOGICAL DRY,WET,MAXDRY
DIMENSION CPA(2),MHUA(2),CPW(2),KW(2),CPV(2),PSYDAT(9),HFG(2),
&A0(2),A1(2),B0(2),B1(2),C0(2),C1(2),MCONV(2)
DATA CPA/1.005, 0.240/, MHUA/0.0659, 0.0443/, CPW/4.186, 1.000/,
. KW/2.043, 0.328/, CPV/1.884, 0.450/,
. HFG/2452., 1054./, PRA/0.72/, PATM/1.0/
DATA TOL/1.E-04/, IMAX/20/, PI/3.14159/, IUNIT/0/
DATA A0/0., 32./, A1/1., 1.8/, B0/0., 7.687/, B1/1.,0.43/
. C0/32., 0./, C1/1.8, 1./, MCONV/1.488, 1.0/
C-----------------------------------------------------------------------------------------------------------------------
C EMODE DETERMINES HOW ERRORS ARE HANDLED FROM CALLS TO THE PSYCH
C SUBROUTINE. IF EMODE IS: 0 - NO ERROR MESSAGES WILL BE PRINTED,
C 1 - ERROR MESSAGES WILL BE PRINTED ONLY ONCE PER SIMULATION,
C 2 - ERROR MESSAGES WILL BE PRINTED EVERY TIMESTEMP THAT THEY OCCUR.
C
DATA EMODE/1/
C
C*********************** STATEMENT FUNCTIONS **************************
C
C-- ROUND-OFF REAL NUMBERS TO NEAREST INTEGER
C
ROUND(RNUM) = JFIX(RNUM + SIGN(0.5,RNUM))
C
C-- UNIT CONVERSIONS FOR TEMPERATURE AND ENTHALPY
C
C TEMPERATURE: SI TO ENGLISH (IU=2), SI TO SI (IU=1)
C ENTHALPY: ENGLISH TO SI (IU=2), SI TO SI (IU=1)
C TEMPERATURE: ENGLISH TO ENGLISH (IU=2), SI TO ENGLISH (IU = 1)
C
TCONV(T) = T*A1(IU) + A0(IU)
HCONV(H) = (H - B0(IU))/B1(IU)
CTOF(T) = T*C1(IU) + C0(IU)
C
C-- CORRELATION FOR SATURATION TEMPERATURE IN TERMS OF SATURATION
C ENTHALPY: CORRELATION IS IN SI UNITS AND IS GOOD FROM 9.473 KJ/KG
C TO 355.137 KJ/KG (0 C TO 55 C).
C (CORRELATION IS FOR A TOTAL SYSTEM PRESSURE OF 1 ATMOSPHERE)
C
TS(HS) = -5.79013 + 6.64030E-01*HS - 5.07802E-03*HS**2 +
. 2.80381E-05*HS**3 - 9.47051E-08*HS**4 +
. 1.72758E-10*HS**5 - 1.29547E-13*HS**6
C
TSAT(HSAT) = TCONV(TS(HCONV(HSAT)))
C
C-- CORRELATION FOR KINEMATIC WATER VISCOSITY IN TERMS OF TEMPERATURE:
C CORRELATION IS IN ENGLISH UNITS (LBM/FT-HR AND DEG. F) AND
C IS GOOD FROM 32 F TO 200 F.
C
MHUW(T) = MCONV(IU)*(7.67201 - 1.39290E-01*T + 1.23322E-03*T**2
. - 5.32425E-06*T**3 + 8.87455E-09*T**4)
C
C-- CORRELATION FOR PRANTL NUMBER OF WATER IN TERMS OF TEMPERATURE:
C TEMPERATURE IS IN ENGLISH UNITS (DEG. F) AND THE CORRELATION IS
C GOOD FROM 32 F TO 200 F.
C
PRW(T) = 2.73327E01 - 6.30480E-01*T + 7.56704E-03*T**2 -
. 5.04074E-05*T**3 + 1.74546E-07*T**4 - 2.43923E-10*T**5
C-----------------------------------------------------------------------------------------------------------------------
C TRNSYS FUNCTIONS
TIME0=getSimulationStartTime()
TFINAL=getSimulationStopTime()
DELT=getSimulationTimeStep()
C SET THE VERSION INFORMATION FOR TRNSYS
IF(INFO(7).EQ.-2) THEN
INFO(12)=16
RETURN 1
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(8).EQ.-1) THEN
RETURN 1
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C PERFORM ANY "AFTER-ITERATION" MANIPULATIONS THAT ARE REQUIRED
IF(INFO(13).GT.0) THEN
RETURN 1
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(7).EQ.-1) THEN
C RETRIEVE THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT FROM THE INFO ARRAY
IUNIT=INFO(1)
ITYPE=INFO(2)
C SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK
INFO(6)=12
INFO(9)=1
INFO(10)=0
C SET THE REQUIRED NUMBER OF INPUTS, PARAMETERS AND DERIVATIVES THAT THE USER SHOULD SUPPLY IN THE INPUT FILE
NPAR=15
C CALL THE TYPE CHECK SUBROUTINE TO COMPARE WHAT THIS COMPONENT REQUIRES TO WHAT IS SUPPLIED
CALL TYPECK(1,INFO,NI,NPAR,ND)
C SET THE VARIABLE TYPES
DATA YCHECK/'TE1','DM1','MF1','TE1','MF1'/
DATA OCHECK/'TE1','DM1','MF1','TE1','MF1','PW1','PW1','PW1',
1 'DM1','PC1','NAV','NAV'/
C CALL THE INPUT-OUTPUT CHECK SUBROUTINE TO SET THE CORRECT INPUT AND OUTPUT UNITS
CALL RCHECK(INFO,YCHECK,OCHECK)
C RETURN TO THE CALLING PROGRAM
RETURN 1
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME
IF (TIME.LT.(TIME0+DELT/2.D0)) THEN
C SET THE UNIT NUMBER FOR FUTURE CALLS
IUNIT=INFO(1)
ITYPE=INFO(2)
C READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
MODE=ROUND(PAR(1))
IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
IU=1
NROWS = ROUND(PAR(2))
NTUBES = ROUND(PAR(3))
LT = PAR(4)
WD = PAR(5)
DO = PAR(6)
DI = PAR(7)
KT = PAR(8)
FT = PAR(9)
FS = PAR(10)
NFINS = ROUND(PAR(11))
KF = PAR(12)
FMODE=ROUND(PAR(13))
IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
ADIST = PAR(14)
IF (FMODE .EQ. 2) DE = ADIST
CDIST = PAR(15)
PSYDAT(1) = PATM
C** COIL AREAS **
C
C RI & RO - INSIDE AND OUTSIDE TUBE RADIUS
C RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR
C CONTINUOUS FLAT FINS
C FH - FIN HEIGHT
C AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS
C FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA
C LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION
C (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS")
C AC - FLOW CROSS SECTIONAL AREA
C AX - CROSS SECTIONAL AREA OF ONE TUBE
C
RI = DI/2.
RO = DO/2.
IF (FMODE .EQ. 1) THEN
RE = DSQRT(ADIST*CDIST/PI)
AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.)
. -(NTUBES*LT*DO-(NFINS*FT*DO))
ELSE
RE = DE/2.
AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT
END IF
LCOIL = NROWS * CDIST
FH = RE - RO
AX = PI*RI**2
AI = NTUBES*NROWS*PI*DI*LT
ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT)
AF = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2)
AO = ATO + AF
FRATIO = AF/AO
C
C** TUBE HEAT TRANSFER COEFFICIENT **
C BASED UPON INSIDE AREA
C
UM = KT/RI/DLOG(RO/RI)
C PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
OUT(1) = XIN(1)
OUT(2) = XIN(2)
OUT(3) = XIN(3)
OUT(4) = XIN(4)
OUT(5) = XIN(5)
OUT(6) = 0.
OUT(7) = 0.
OUT(8) = 0.
OUT(9) = 1.0
OUT(10) = 0.0
C RETURN TO THE CALLING PROGRAM
RETURN 1
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C *** ITS AN ITERATIVE CALL TO THIS COMPONENT ***
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED SINCE THE LAST TIME THE PARAMETERS
C WERE READ IN
IF(INFO(1).NE.IUNIT) THEN
C RESET THE UNIT NUMBER
IUNIT=INFO(1)
ITYPE=INFO(2)
C READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
MODE=ROUND(PAR(1))
IF(MODE.NE.1 .AND. MODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
IU=1
NROWS = ROUND(PAR(2))
NTUBES = ROUND(PAR(3))
LT = PAR(4)
WD = PAR(5)
DO = PAR(6)
DI = PAR(7)
KT = PAR(8)
FT = PAR(9)
FS = PAR(10)
NFINS = ROUND(PAR(11))
KF = PAR(12)
FMODE=ROUND(PAR(13))
IF(FMODE.NE.1 .AND. FMODE.NE.2) CALL TYPECK(-4,INFO,0,0,0)
ADIST = PAR(14)
IF (FMODE .EQ. 2) DE = ADIST
CDIST = PAR(15)
PSYDAT(1) = PATM
C** COIL AREAS **
C
C RI & RO - INSIDE AND OUTSIDE TUBE RADIUS
C RE - ANNULAR FIN RADIUS OF EQUIVALENT RADIUS FOR
C CONTINUOUS FLAT FINS
C FH - FIN HEIGHT
C AI & AO - INSIDE AND OUTSIDE COIL SURFACE AREAS
C FRATIO - RATIO OF FIN AREA TO OUTSIDE COIL SURFACE AREA
C LCOIL - LENGTH OF THE COIL SECTION IN THE AIR FLOW DIRECTION
C (SEE KAYS & LONDON, "COMPACT HEAT EXCHANGERS")
C AC - FLOW CROSS SECTIONAL AREA
C AX - CROSS SECTIONAL AREA OF ONE TUBE
C
RI = DI/2.
RO = DO/2.
IF (FMODE .EQ. 1) THEN
RE = DSQRT(ADIST*CDIST/PI)
AC = LT*WD-NFINS*FT*(NTUBES*ADIST+ADIST/2.)
. -(NTUBES*LT*DO-(NFINS*FT*DO))
ELSE
RE = DE/2.
AC = LT*WD-NTUBES*LT*DO-NTUBES*NFINS*(DE-DO)*FT
END IF
LCOIL = NROWS * CDIST
FH = RE - RO
AX = PI*RI**2
AI = NTUBES*NROWS*PI*DI*LT
ATO = NTUBES*NROWS*PI*DO*(LT-NFINS*FT)
AF = NFINS*NTUBES*NROWS*2.*PI*(RE**2-RO**2)
AO = ATO + AF
FRATIO = AF/AO
C
C** TUBE HEAT TRANSFER COEFFICIENT **
C BASED UPON INSIDE AREA
C
UM = KT/RI/DLOG(RO/RI)
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C****************************** INPUTS ****************************
C
C TA1 = ENTERING AIR DRY BULB TEMPERATURE (C, F)
C WA1 = ENTERING AIR HUMIDITY RATIO
C MA = MASS FLOW RATE OF AIR (KG/HR, LBM/HR)
C TW1 = ENTERING WATER TEMPERATURE (C, F)
C MW = MASS FLOW RATE OF WATER (KG/HR LBM/HR)
C
TA1 = XIN(1)
WA1 = XIN(2)
MA = XIN(3)
TW1 = XIN(4)
MW = XIN(5)
C CHECK THE INLET AIR STATE
PSYDAT(1) = PATM
PSYDAT(2) = TA1
PSYDAT(6) = WA1
CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*99)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
99 TA1 = PSYDAT(2)
WA1 = PSYDAT(6)
RHA1 = PSYDAT(4)
C-----------------------------------------------------------------------------------------------------------------------
C************************* COIL ANALYSIS *************************
C
C TA2 = DRY-BULB TEMPERATURE OF AIR LEAVING COIL
C WA2 = ABSOLUTE HUMIDITY OF AIR LEAVING COIL
C RHA2 = RELATIVE HUMIDITY OF AIR LEAVING COIL
C TW2 = DRY-BULB TEMPERATURE OF WATER LEAVING COIL
C QCOIL = TOTAL HEAT TRANSFER TO COIL
C QSENS = SENSIBLE HEAT REMOVED FROM MOIST AIR
C QLAT = LATENT HEAT REMOVED FROM MOIST AIR
C QDAIR = SENSIBLE HEAT REMOVED FROM AIR ONLY (DRY AIR)
C FDRY = FRACTION OF COIL WHICH IS DRY
C GW = WATER MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA
C REI = INSIDE REYNOLDS NUMBER
C HI = INSIDE COEFFICIENT OF CONVECTION
C GA = AIR MASS FLOW RATE DIVIDED BY CROSS-SECTIONAL AREA
C REO = OUTSIDE REYNOLDS NUMBER
C CPM = MOIST AIR HEAT CAPACITY
C NTUO = NTU'S OUTSIDE
C NTUI = NTU'S INSIDE
C TDP = AIR DEW-POINT
C HA1 = ENTHALPY OF ENTERING AIR
C HA2 = ENTHALPY OF EXITING AIR
C *** NO FLOW ***
IF (MA .LE. 1.0E-06 .OR. MW .LE. 1.0E-06) THEN
TA2 = TA1
IF (MA .LE. 1.0E-02 .AND. MW .GT. 1.0E-02) TA2 = TW1
WA2 = WA1
RHA2 = RHA1
TW2 = TW1
QCOIL = 0.0
QSENS = 0.0
QLAT = 0.0
QDAIR = 0.0
FDRY = 1.0
GOTO 500
END IF
C *** FLOW ***
C ** INSIDE FLUID COEFFICIENT **
C FIND THE REYNOLDS NUMBER
GW = MW/NTUBES/AX
TPROP = CTOF(TW1)
IF ((32.0.GT.TPROP .OR. TPROP.GT.200.) .AND.
. (TIME-OUT(12)) .GT. DELT/2.) THEN
CALL MESSAGES(-1,'THE CORRELATIONS FOR THE VISCOSITY AND PRANDT
&L # OF WATER USED IN THE COOLING COIL MODEL WERE USED WITH A VALUE
& OUTSIDE OF THEIR INTENDED RANGE','WARNING',INFO(1),INFO(2))
IWARN=IWARN+1
OUT(12) = TIME
END IF
REI = GW*DI/MHUW(TPROP)
C Added by Avatar Lee on 17/07/2013
C Calculation of laminar Nusselt number based on Nellis G, Klein S. Heat Transfer. New York: Cambridge University Press, 2009.
GZ = DI*MIN(REI,2300.0D0)*PRW(TPROP)/LT
NU_LAM = 3.66+(0.049+0.02/PRW(TPROP))*GZ**1.12/(1.0+0.065*GZ**0.7)
C Calculation of turbulent Nusselt number based on Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow.
C International Chemical Engineering 1976;16:359-68.
RE_TUR = MAX(REI,3000.0D0)
FD = 4.0*(-0.0015702/LOG(RE_TUR)+0.3942031/(LOG(RE_TUR))**2.0
&+2.5341533/(Log(RE_TUR))**3.0)
NU_TUR = (FD/8.0)*(RE_TUR-1000.0)*PRW(TPROP)*(1.0+(DI/LT)**0.7)
&/(1.0+12.7*SQRT(FD/8.0)*(PRW(TPROP)**(2.0/3.0)-1.0))
C Revision ended...
C IF THE REYNOLDS NUMBER IS GREATER THEN 3000, USE THE CORRELATION
C FOR TURBULENT WATER FLOW IN A TUBE. IF THE REYNOLDS NUMBER IS
C LESS THAN 2000, USE THE CORRELATION FOR LAMINAR WATER FLOW IN A
C TUBE. FOR A REYNOLDS NUMBER BETWEEN 2000 AND 3000, A LINEAR
C RELATIONSHIP BETWEEN THE CONVECTION COEFFICIENTS FOR THE TURBULENT
C AND LAMINAR CASES IS USED. THE INLET TEMPERATURE IS USED FOR
C PROPERTY EVALUATIONS.
IF (REI .GE. 3000.) THEN
C HI = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4
HI = NU_TUR*KW(IU)/DI !Added by Avatar Lee on 17/07/2013
ELSE IF (REI .LE. 2300.) THEN
C HI = 4.36*KW(IU)/DI
HI = NU_LAM*KW(IU)/DI !Added by Avatar Lee on 17/07/2013
ELSE
C HIT = 0.023*KW(IU)/DI*REI**0.8*PRW(TPROP)**0.4
HIT = NU_TUR*KW(IU)/DI !Added by Avatar Lee on 17/07/2013
C HIL = 4.36*KW(IU)/DI
HIL = NU_LAM*KW(IU)/DI !Added by Avatar Lee on 17/07/2013
HI = HIL + (REI - 2000.0) * (HIT - HIL) / 1000.0
END IF
C ** OUTSIDE COEFFICIENT FOR DRY SURFACES **
C CORRELATION FOR OUTSIDE HEAT TRANSFER COEFFICIENT IS FROM
C "FINNED TUBE HEAT EXCHANGER: CORRELATION OF DRY SURFACE HEAT
C TRANSFER DATA," A.H. ELMAHDY, ASHRAE TRANSACTIONS. THE
C HYDRAULIC DIAMETER IS: 2 * THE FLOW LENGTH OF THE HEAT EXCHANGER *
C THE FLOW CROSS SECTIONAL AREA / TOTAL HEAT TRANSFER AREA.
C THE FIN EFFICIENCY EQUATION IS FOR ANNULAR FINS. STRAIGHT FINS
C ARE TREATED AS ANNULAR FINS BY FINDING AN EQUIVALENT ANNULAR FIN
C WITH THE SAME SURFACE AREA.
DH = 4.*LCOIL*AC/AO
CJ1 = 0.159*(FT/FH)**0.141*(DH/FT)**0.065
CJ2 = -0.323*(FT/FH)**0.049*(FS/FT)**0.077
GA = MA/AC
REO = GA*DH/MHUA(IU)
JF = CJ1*REO**CJ2
HDO = GA*CPA(IU)*JF/PRA**0.6667
C FIN EFFICIENCY CALCULATIONS
ALPHA = RO/RE
BETA = RE*DSQRT(2.*HDO/KF/FT)
EFF = FIN(ALPHA,BETA)
EFFD = 1. - FRATIO*(1.-EFF)
C ** DETERMINE NTU'S **
C THE AIR SPECIFIC HEAT IS FOR MOIST AIR
CPM = CPA(IU) + WA1*CPV(IU)
NTUO = EFFD*HDO*AO/(MA*CPM)
NTUI = AI/(1./HI + 1./UM)/(MW*CPW(IU))
RA = MW/MA
C ** ENTHALPIES FOR ENTERING CONDITIONS **
PSYDAT(2) = TA1
PSYDAT(6) = WA1
CALL PSYCHROMETRICS(TIME,INFO,IU,4,0,PSYDAT,EMODE,ISTAT,*101)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
101 TDP = PSYDAT(5)
HA1 = PSYDAT(7)
PSYDAT(2) = TW1
PSYDAT(3) = TW1
CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*102)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
102 HW1 = PSYDAT(7)
C ** DRY ANALYSIS **
C DETERMINES THE AIR EFFECTIVENESS (EPSD) AND OUTLET WATER TEMPERATURE
C (TW2D) ASSUMING THE COIL IS DRY
C IF THE SURFACE TEMPERATURE AT THE OUTLET IS GREATER THAN THE AIR
C ENTERING DEWPOINT THEN COIL IS COMPLETELY DRY, OTHERWISE DO THE WET
C ANALYSIS
CSTAR = CPM/RA/CPW(IU)
NTUD = NTUO/(1. + NTUO/NTUI*CSTAR)
TEMP = -NTUD*(1.-CSTAR)
IF (TEMP .LT. 50.0) THEN
C = DEXP(TEMP)
EPSD = (1. - C)/(1. - CSTAR*C)
ELSE
EPSD = 1/CSTAR
ENDIF
TW2D = TW1 + EPSD*CSTAR*(TA1 - TW1)
TA2D = TA1 - EPSD*(TA1 - TW1)
TS2D = TW1 + CSTAR*NTUD/NTUI*(TA2D - TW1)
C
DRY = TS2D .GT. TDP
FDRY = 1.0
IF(.NOT. DRY) THEN
FDRY = 0.0
C ** WET COIL **
C DETERMINES THE AIR EFFECTIVENESS (EPSW) AND OUTLET WATER TEMPERATURE
C (TW2W) ASSUMING THE COIL IS WET.
C THE AVERAGE SATURATION SPECIFIC HEAT (CS) DEPENDS UPON THE OUTLET
C CONDITIONS; THEREFORE AN ITERATIVE SOLUTION.
C THE WET SURFACE HEAT TRANSFER COEFFICIENT (HWO) AND THEREFORE THE FIN
C EFFICIENCY (EFFW) DEPEND UPON THE SATURATION SPECIFIC HEAT.
C IF THE SURFACE TEMPERATURE AT WATER OUTLET IS LESS THAN THE DEWPOINT
C OF THE ENTERING AIR THEN THE COIL IS COMPLETELY WET.
TW2W = TW2D
ITER = 0
10 ITER = ITER + 1
PSYDAT(2) = TW2W
PSYDAT(3) = TW2W
CALL PSYCHROMETRICS(TIME,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*104)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
104 HW2 = PSYDAT(7)
CS = (HW2 - HW1)/(TW2W - TW1)
MSTAR = CS/RA/CPW(IU)
HWO = CS*HDO/CPM
ALPHA = RO/RE
BETA = RE*DSQRT(2.*HWO/KF/FT)
EFF = FIN(ALPHA,BETA)
EFFW = 1. - FRATIO*(1.-EFF)
ERAT = EFFW/EFFD
NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR)
TEMP = -NTUW*(1.-MSTAR)
IF (TEMP .LT. 50.0) THEN
C = DEXP(TEMP)
EPSW = (1. - C)/(1. - MSTAR*C)
ELSE
EPSW = 1/MSTAR
ENDIF
C USE A SECANT METHOD TO CONVERGE ON TW2
G=TW2W-(TW1 + EPSW*(HA1 - HW1)/RA/CPW(IU))
IF (DABS(G) .GT. TOL .AND. ITER .LT. IMAX) THEN
IF (ITER .EQ. 1) THEN
G1 = G
TOLD = TW2W
TW2W = TOLD - G1
ELSE
G2 = G
DGDT = (G2 - G1)/(TW2W - TOLD)
DGDT = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT)
TOLD = TW2W
TW2W = TOLD - G2/DGDT
IF (DABS(TW2W-TOLD) .LT. TOL) TW2W = TOLD - G2
G1 = G2
END IF
GOTO 10
END IF
TS1W = TW2W + CSTAR/CPM*NTUW/NTUI*(HA1 - HW2)
ENDIF
C CHECK TO SEE IF COIL COMPLETELY WET
WET = ((.NOT. DRY) .AND. TS1W .LT. TDP)
IF(.NOT.(DRY) .AND. .NOT.(WET)) THEN
IF (MODE .EQ. 1) THEN
DRY = TW2D .GT. TW2W
FDRY = -1.0
ELSE IF (MODE .EQ. 2) THEN
C ** COMBINED WET AND DRY COIL PERFORMANCE **
C
C AN ITERATIVE BY-SECTION METHOD IS USEDD TO FIND TW2. A FUNCTION
C G, WHICH EQUALS ZERO AT THE SOLUTION, IS DEFINED FROM AN ENERGY
C BALANCE. THE VALUE OF G IS FOUND FOR TW2 CORRESPONDING TO FDRY=0
C AND FDRY=1. IF BOTH VALUES OF G ARE POSITIVE OR BOTH ARE NEGATIVE,
C THE LOGICAL VARIABLES "WET" AND "DRY" ARE APPROPRIATELY SET. (IE.
C THE COIL WILL BE CONSIDERED WET IF THE WET LIMIT GIVES THE G VALUE
C CLOSEST TO 0, OR THE COIL WILL BE CONSIDERED DRY IF THE DRY LIMIT
C GIVES THE G VALUE CLOSEST TO 0.)
CK = NTUD*(1. - CSTAR)
LIMITW=(TA1*(CSTAR-1.+CK/NTUO)+TDP*(1-CSTAR))/(CK/NTUO)
LIMITD=(TA1*(CSTAR-(1.-CK/NTUO)*DEXP(-CK))+TDP*(1.-CSTAR))/
& (1.-(1.-CK/NTUO)*DEXP(-CK))
IF(LIMITD .GT. LIMITW) THEN
TMAX = 0.999999 * LIMITD
TMIN = 1.000001 * LIMITW
MAXDRY = .TRUE.
ELSE
TMAX = 0.999999 * LIMITW
TMIN = 1.000001 * LIMITD
MAXDRY = .FALSE.
ENDIF
CALL CALCG(TMAX,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
& NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO,
& IU,GMAX,FDRY,EPSD,TWX,ERAT,EMODE)
CALL CALCG(TMIN,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
& NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,INFO,
& IU,GMIN,FDRY,EPSD,TWX,ERAT,EMODE)
IF ((GMAX .GT. 0.0 .AND. GMIN .GT. 0.0) .OR.
. (GMAX .LT. 0.0 .AND. GMIN .LT. 0.0)) THEN
IF (DABS(GMAX) .LT. DABS(GMIN) ) THEN
DRY = MAXDRY
WET = .NOT. MAXDRY
ELSE
WET = MAXDRY
DRY = .NOT. MAXDRY
ENDIF
ELSE
TW2 = TMIN
ITER = -1
C ITERATION LOOP TO DETERMINE TW2.
15 ITER = ITER + 1
CALL CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
& NTUI,NTUD,HDO,CPM,CSTAR,CPW(IU),EFFD,FRATIO,
& INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE)
IF ((DABS(G) .GT. TOL .AND. ITER .LT. IMAX) .AND.
. (DABS(TMIN - TMAX) .GT. TOL)) THEN
IF (ITER .EQ. 0) THEN
SGNMIN = SIGN(1.,G)
C USE THE MAXIMUM TW2 FROM THE WET AND DRY ANALYSIS AS AN INITIAL GUESS
TW2 = DMAX1(TMIN,DMIN1(DMAX1(TW2W,TW2D),TMAX))
GOTO 15
ELSE IF (ITER .EQ. 1) THEN
G1 = G
TOLD = TW2
TW2 = TOLD - G1
ELSE
G2 = G
DGDT = (G2 - G1)/(TW2 - TOLD)
DGDT = SIGN((DMAX1(1.0E-06,DABS(DGDT))),DGDT)
TOLD = TW2
TW2 = TOLD - G2/DGDT
G1 = G2
END IF
C NARROW THE BOUNDS TO THE SOLUTION
IF (G1 .LT. 0.) THEN
IF (SGNMIN .LT. 0.) THEN
TMIN = TOLD
ELSE
TMAX = TOLD
END IF
ELSE
IF (SGNMIN .GT. 0.) THEN
TMIN = TOLD
ELSE
TMAX = TOLD
END IF
END IF
C USE A BI-SECTION IF THE NEW VALUE OF TW2 IS OUTSIDE THE BOUNDS
IF (TW2.GE.TMAX.OR.TW2.LE.TMIN) TW2 = (TMAX + TMIN)/2.
GOTO 15
END IF
END IF
C * END OF COMBINED COIL ANALYSIS
END IF
END IF
C ** END OF COIL ANALYSIS
C
C *** NET COIL PERFORMANCE ***
C FOR MODE 1, UTILIZE THE RESULTS OF THE ANALYSIS THAT GIVES THE
C GREATEST HEAT TRANSFER (I.E. HIGHEST WATER RETURN TEMPERATURE)
C IF COIL IS NOT COMPLETELY WET OR DRY
C ** DRY **
IF(DRY) THEN
TW2 = TW2D
TA2 = TA2D
HA2 = HA1 - CPM*(TA1 - TA2)
WA2 = WA1
RHA2 = RHA1
ELSE IF ((WET) .OR. MODE .EQ. 1) THEN
C ** WET **
TW2 = TW2W
HA2 = HA1 - EPSW*(HA1 - HW1)
HSAVG = HA1 - (HA1 - HA2)/(1. - DEXP(-ERAT*NTUO))
IF ((9.473.GT.HCONV(HSAVG) .OR. HCONV(HSAVG).GT.355.137)
. .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN
CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA
&TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS
&IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2))
IWARN=IWARN+1
OUT(11) = TIME
END IF
TSAVG = TSAT(HSAVG)
TA2 = TSAVG + (TA1 - TSAVG)*DEXP(-ERAT*NTUO)
PSYDAT(2) = TA2
PSYDAT(7) = HA2
CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*105)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
105 WA2 = PSYDAT(6)
RHA2 = PSYDAT(4)
ELSE
C ** COMBINED **
HA2 = HA1 - RA*CPW(IU)*(TW2 - TW1)
TAX = TA1 - EPSD*(TA1 - TWX)
HAX2 = HA1 - RA*CPW(IU)*(TW2 - TWX)
HAX = HA1 - CPM*(TA1 - TAX)
HSW = HAX + (HA2 - HAX)/(1. - DEXP(-(1. - FDRY)*ERAT*NTUO))
IF ((9.473.GT.HCONV(HSW) .OR. HCONV(HSW).GT.355.137)
. .AND. (TIME-OUT(11)) .GT. DELT/2.) THEN
CALL MESSAGES(-1,'THE CORRELATION FOR THE SATURATION TEMPERA
&TURE USED IN THE COOLING COIL MODEL WAS USED WITH AN ENTHALPY OUTS
&IDE OF ITS INTENDED RANGE','WARNING',INFO(1),INFO(2))
IWARN=IWARN+1
OUT(11) = TIME
END IF
TSW = TSAT(HSW)
TA2 = TSW + (TAX - TSW)*DEXP(-(1.-FDRY)*ERAT*NTUO)
PSYDAT(2) = TA2
PSYDAT(7) = HA2
CALL PSYCHROMETRICS(TIME,INFO,IU,5,0,PSYDAT,EMODE,ISTAT,*106)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
106 WA2 = PSYDAT(6)
RHA2 = PSYDAT(4)
ENDIF
C-----------------------------------------------------------------------------------------------------------------------
C SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT
C****************************** OUTPUTS *****************************
C
QCOIL = MW*CPW(IU)*(TW2 - TW1)
QLAT = MA*(WA1 - WA2)*HFG(IU)
QSENS = QCOIL - QLAT
C
500 OUT(1) = TA2
OUT(2) = WA2
OUT(3) = MA
OUT(4) = TW2
OUT(5) = MW
OUT(6) = QCOIL
OUT(7) = QSENS
OUT(8) = QLAT
OUT(9) = FDRY
OUT(10) = 100.*DMIN1(1.,(DMAX1(0.,RHA2)))
C-----------------------------------------------------------------------------------------------------------------------
C EVERYTHING IS DONE - RETURN FROM THIS SUBROUTINE AND MOVE ON
RETURN 1
END
C-----------------------------------------------------------------------------------------------------------------------
C THIS SUBROUTINE CALCULATES G (WHICH IS BASED ON AN ENERGY BALANCE ON THE
C WET SIDE OF THE COIL), FDRY AND EPSD.
C SEE VARIABLE DEFINITIONS IN THE MAIN PROGRAM. CPW IS ASSUMED TO BE
C PASSED IN THE CORRECT UNITS.
C
SUBROUTINE CALCG(TW2,TDP,TA1,TW1,HA1,HW1,RO,RA,RE,KF,FT,NTUO,
. NTUI,NTUD,HDO,CPM,CSTAR,CPW,EFFD,FRATIO,
. INFO,IU,G,FDRY,EPSD,TWX,ERAT,EMODE)
IMPLICIT NONE
INTEGER*4 INFO
DIMENSION PSYDAT(9),INFO(15)
DOUBLE PRECISION KF,NTUD,NTUW,NTUO,NTUI,MSTAR,PATM,EXPK,TDP,TW2,
& CSTAR,TA1,FDRY,TEMP,C,EPSD,TWX,TW1,PSYDAT,HWX,CS,HW1,RA,HWO,HDO,
& CPM,ALPHA,RO,RE,BETA,FT,EFF,FIN,EFFW,FRATIO,ERAT,EFFD,CPW,EPSW,G,
& HA1
INTEGER EMODE,IU,ISTAT
DATA PATM/1.0/
EXPK = ((TDP-TW2)+CSTAR*(TA1-TDP))/
& (1.-NTUD*(1.-CSTAR)/NTUO)/(TA1-TW2)
EXPK = DMAX1(1.0E-6,EXPK)
FDRY = -1./NTUD/(1.-CSTAR)*DLOG(EXPK)
FDRY = DMAX1(0.D0,DMIN1(1.D0,FDRY))
TEMP = -FDRY*NTUD*(1. - CSTAR)
IF (TEMP .LT. 50.0) THEN
C = DEXP(TEMP)
EPSD = (1. - C)/(1. - CSTAR*C)
ELSE
EPSD = 1/CSTAR
ENDIF
TWX = (TW2 - CSTAR*EPSD*TA1)/(1 - CSTAR*EPSD)
TWX = DMAX1(TWX,TW1+1.E-04)
PSYDAT(1) = PATM
PSYDAT(2) = TWX
PSYDAT(3) = TWX
CALL PSYCHROMETRICS(1.0D0,INFO,IU,1,0,PSYDAT,EMODE,ISTAT,*107)
CALL LINKCK('TYPE52','PSYCHROMETRICS ',1,99)
107 HWX = PSYDAT(7)
CS = (HWX - HW1)/(TWX - TW1)
MSTAR= CS/RA/CPW
HWO = CS*HDO/CPM
ALPHA = RO/RE
BETA = RE*DSQRT(2.*HWO/KF/FT)
EFF = FIN(ALPHA,BETA)
EFFW = 1. - FRATIO*(1.-EFF)
ERAT = EFFW/EFFD
NTUW = ERAT*NTUO/(1. + ERAT*NTUO/NTUI*MSTAR)
TEMP = -(1. - FDRY)*NTUW*(1. - MSTAR)
IF (TEMP .LT. 50.0) THEN
C = DEXP(TEMP)
EPSW = (1. - C)/(1. - MSTAR*C)
ELSE
EPSW = 1/MSTAR
ENDIF
G=TWX*(1.-CSTAR*EPSW*EPSD)-TW1-CSTAR/CPM*EPSW*(HA1-EPSD*CPM*TA1
& -HW1)
RETURN
END
C******************************************************************************
C THIS FUNCTION CALCULATES THE FIN EFFICIENCY (EFFECTIVENESS)
C OF AN ANNULAR FIN OF CONSTANT THICKNESS.
C
C ALPHA = RADIUS AT FIN BASE / RADIUS AT FIN TIP
C BETA = RADIUS AT FIN TIP *
C (SQRT (2 * CONVECTION COEFFICIENT /
C FIN CONDUCTIVITY * FIN THICKNESS))
FUNCTION FIN(ALPHA,BETA)
IMPLICIT NONE
DOUBLE PRECISION FIN,I0,I1,K0,K1,ALPBET,ALPHA,BETA,XI0,XI1,XK0,
& XK1,YI0,YI1,YK0,YK1
ALPBET = ALPHA * BETA
CALL BESSEL(ALPBET,I0,I1,K0,K1)
XI0 = I0
XI1 = I1
XK0 = K0
XK1 = K1
CALL BESSEL(BETA,I0,I1,K0,K1)
YI0 = I0
YI1 = I1
YK0 = K0
YK1 = K1
FIN=2.*ALPHA/BETA/(1.- ALPHA**2)*(XK1*YI1-XI1*YK1)/(XK0*YI1
& +XI0*YK1)
RETURN
END
C****************************************************************************
C THIS SUBROUTINE USES POLYNOMIAL APPROXIMATIONS TO EVALUATE
C THE BESSEL FUNCTIONS. THE APPROXIMATIONS ARE FROM ABRAMOWITZ
C AND STEGUN, HANDBOOD OF MATHEMATICAL FUNCTIONS, DOVER
C PUBLICATIONS, INC., NEW YORK, NY.
C
SUBROUTINE BESSEL(X,I0,I1,K0,K1)
IMPLICIT NONE
DOUBLE PRECISION X,I0,I1,K0,K1,IT,A0,A1,A2,A3,A4,A5,A6,B0,B1,B2,
& B3,B4,B5,B6,B7,B8,C0,C1,C2,C3,C4,C5,C6,D0,D1,D2,D3,D4,D5,D6,D7,
& D8,E0,E1,E2,E3,E4,E5,E6,F0,F1,F2,F3,F4,F5,F6,G0,G1,G2,G3,G4,G5,
& G6,H0,H1,H2,H3,H4,H5,H6,T,TT,X1,X2
C THE FOLLOWING DATA STATEMENTS CONTAIN THE COEFFICIENTS TO
C THE POLYNOMIALS.
C I0
DATA A0/1.0/,A1/3.5156229/,A2/3.0899424/,A3/1.2067492/
DATA A4/0.2659732/,A5/0.0360768/,A6/0.0045813/
C I0
DATA B0/0.39894228/,B1/0.01328592/,B2/0.00225319/
DATA B3/-0.00157565/,B4/0.00916281/,B5/-0.02057706/
DATA B6/0.02635537/,B7/-0.01647633/,B8/0.00392377/
C I1
DATA C0/0.5/,C1/0.87890594/,C2/0.51498869/,C3/0.15084934/
DATA C4/0.02658733/,C5/0.00301532/,C6/0.00032411/
C I1
DATA D0/0.39894228/,D1/-0.03988024/,D2/-0.00362018/
DATA D3/0.00163801/,D4/-0.01031555/,D5/0.02282967/
DATA D6/-0.02895312/,D7/0.01787654/,D8/-0.00420059/
C K0
DATA E0/-0.57721566/,E1/0.4227842/,E2/0.23069756/
DATA E3/0.0348859/,E4/0.00262698/,E5/0.0001075/,E6/0.0000074/
C K0
DATA F0/1.25331414/,F1/-0.07832358/,F2/0.02189568/
DATA F3/-0.01062446/,F4/0.00587872/,F5/-0.0025154/
DATA F6/0.00053208/
C K1
DATA G0/1.0/,G1/0.15443144/,G2/-0.67278579/,G3/-0.18156897/
DATA G4/-0.01919402/,G5/-0.00110404/,G6/-0.00004686/
C K1
DATA H0/1.25331414/,H1/0.23498619/,H2/-0.0365562/
DATA H3/0.01504268/,H4/-0.00780353/,H5/0.00325614/
DATA H6/-0.00068245/
C
IF (X .LT. -3.75) THEN
CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C
&OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL',
&-1,-1)
RETURN
END IF
T=X/3.75
TT=T*T
C I0
IF (X .LE. 3.75) THEN
I0=A0+TT*(A1+TT*(A2+TT*(A3+TT*(A4+TT*(A5+TT*A6)))))
ELSE
IT=1/T
I0=(B0+IT*(B1+IT*(B2+IT*(B3+IT*(B4+IT*(B5+IT*(B6+IT*
. (B7+IT*B8))))))))/(DSQRT(X)*DEXP(-X))
END IF
C I1
IF (X .LE. 3.75) THEN
I1=(C0+TT*(C1+TT*(C2+TT*(C3+TT*(C4+TT*(C5+TT*C6))))))*X
ELSE
IT=1/T
I1=(D0+IT*(D1+IT*(D2+IT*(D3+IT*(D4+IT*(D5+IT*(D6+IT*
. (D7+IT*D8))))))))/(DSQRT(X)*DEXP(-X))
END IF
C K0
IF (X .LE. 0.0) THEN
CALL MESSAGES(-1,'THE BESSEL FUNCTION CALLED FROM THE COOLING C
&OIL SUBROUTINE COULD NOT BE EVALUATED AT THE GIVEN VALUE','FATAL',
&-1,-1)
RETURN
END IF
X1 = (X/2.)**2
X2 = 2./X
IF (X .LE. 2.0) THEN
K0=-DLOG(X/2)*I0+E0+X1*(E1+X1*(E2+X1*(E3+X1*(E4+X1*
. (E5+X1*E6)))))
ELSE
K0=(F0+X2*(F1+X2*(F2+X2*(F3+X2*(F4+X2*(F5+X2*F6))))))
. /(DSQRT(X)*DEXP(X))
END IF
C K1
IF (X .LE. 2.0) THEN
K1=(X*DLOG(X/2.)*I1+G0+X1*(G1+X1*(G2+X1*(G3+X1*(G4+X1*
. (G5+X1*G6))))))/X
ELSE
K1=(H0+X2*(H1+X2*(H2+X2*(H3+X2*(H4+X2*(H5+X2*H6))))))
. /(DSQRT(X)*DEXP(X))
END IF
RETURN
END