PROGRAM H3SS1 C... C... H3 - HEAT CONDUCTION IN A CIRCULAR FIN, STEADY STATE C... SOLUTION, VERSION 1 C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION(A-H,O-Z) C... C... ADDITIONAL VARIABLES ARE DECLARED DOUBLE PRECISION (TO C... ACCOMMODATE THE MODIFIED BESSEL FUNCTIONS) DOUBLE PRECISION I00, I01, I11, + K00, K01, K11 C... C... COMMON AREA WITH PROBLEM PARAMETERS PARAMETER(NR=6) COMMON/P/ ALPHA, BETA, UA, U0, U(NR), + DR, R0, R1, R(NR), B, + RHOCP C... C... OPEN AN OUTPUT FILE NO=6 OPEN(NO,FILE='OUTPUT', STATUS='UNKNOWN') C... C... CALL SUBROUTINE BFTEST TO CHECK THE PERFORMANCE OF FUNCTIONS C... BESSI0, BESSI1, BESSK0 AND BESSK1 CALL BFTEST(NO) C... C... PROBLEM PARAMETERS C... C... FIN INSIDE RADIUS R0=0.5D0 C... C... FIN OUTSIDE RADIUS R1=2.0D0 C... C... FIN THICKNESS B=0.2D0 C... C... HEAT TRANSFER COEFFICIENT H=0.01D0 C... C... RHO*CP RHOCP=2.0D0 C... C... THERMAL DIFFUSIVITY ALPHA=0.01D0 C... C... AMBIENT TEMPERATURE UA=25.0D0 C... C... INSIDE TEMPERATURE U0=150.0D0 C... C... RECIPROCAL TIME FOR HEAT TRANSFER BETA=(2.0D0*H)/(B*RHOCP) C... C... INCREMENT FOR RADIAL GRID DR=(R1-R0)/DFLOAT(NR-1) C... C... RADIAL GRID DO 1 I=1,NR R(I)=R0+DFLOAT(I-1)*DR 1 CONTINUE C... C... STEADY STATE SOLUTION WITH DIRICHLET BOUNDARY CONDITIONS C... C... CONSTANTS C1, C2 ARG0=DSQRT(BETA/ALPHA)*R0 ARG1=DSQRT(BETA/ALPHA)*R1 I00=BESSI0(ARG0) I01=BESSI0(ARG1) K00=BESSK0(ARG0) K01=BESSK0(ARG1) DEN=I00*K01-I01*K00 C1= K01/DEN C2=-I01/DEN C... C... EVALUATE THE SOLUTION AT A SERIES OF RADIAL POSITIONS DO 2 I=1,NR ARG=DSQRT(BETA/ALPHA)*R(I) U(I)=C1*BESSI0(ARG)+C2*BESSK0(ARG) C... C... CONVERSION FROM DIMENSIONLESS TO DIMENSIONAL TEMPERATURE U(I)=UA+(U0-UA)*U(I) 2 CONTINUE C... C... PRINT THE SOLUTION WRITE(NO,3)(R(I),I=1,NR),(U(I),I=1,NR) 3 FORMAT(/,' Dirichlet boundary conditions',/, + ' r ',6F9.2,/,' u(r) ',6F9.3,/) C... C... COMPUTE AND PRINT THE TOTAL HEAT TRANSFER TO THE FIN DTDR=(C1*BESSI1(ARG0)-C2*BESSK1(ARG0))*DSQRT(BETA/ALPHA)*(U0-UA) PI=4.0D0*DATAN(1.0D0) TOTAL=-2.0D0*PI*R0*B*ALPHA*RHOCP*DTDR WRITE(NO,4)TOTAL 4 FORMAT(' Total heat transfer rate = ',F7.3,//) C... C... STEADY STATE SOLUTION WITH DIRICHLET-NEUMANN BOUNDARY C... CONDITIONS C... C... CONSTANTS C1, C2 ARG0=DSQRT(BETA/ALPHA)*R0 ARG1=DSQRT(BETA/ALPHA)*R1 I00=BESSI0(ARG0) I11=BESSI1(ARG1) K00=BESSK0(ARG0) K11=BESSK1(ARG1) DEN=I00*K11+I11*K00 C1=K11/DEN C2=I11/DEN C... C... EVALUATE THE SOLUTION AT A SERIES OF RADIAL POSITIONS DO 12 I=1,NR ARG=DSQRT(BETA/ALPHA)*R(I) U(I)=C1*BESSI0(ARG)+C2*BESSK0(ARG) C... C... CONVERSION FROM DIMENSIONLESS TO DIMENSIONAL TEMPERATURE U(I)=UA+(U0-UA)*U(I) 12 CONTINUE C... C... PRINT THE SOLUTION WRITE(NO,13)(R(I),I=1,NR),(U(I),I=1,NR) 13 FORMAT(/,' Dirichlet-Neumann boundary conditions',/, + ' r ',6F9.2,/,' u(r) ',6F9.3,/) C... C... COMPUTE AND PRINT THE TOTAL HEAT TRANSFER TO THE FIN DTDR=(C1*BESSI1(ARG0)-C2*BESSK1(ARG0))*DSQRT(BETA/ALPHA)*(U0-UA) PI=4.0D0*DATAN(1.0D0) TOTAL=-2.0D0*PI*R0*B*ALPHA*RHOCP*DTDR WRITE(NO,14)TOTAL 14 FORMAT(' Total heat transfer rate = ',F7.3,//) STOP END SUBROUTINE BFTEST(NO) C... C... PROGRAM BFTEST CALLS FUNCTIONS BESSI0, BESSI1, BESSK0 AND C... BESSK1 TO CHECK THE CALCULATION OF I0(X), I1(X), K0(X) AND C... K1(X). C... C... THE FOLLOWING TABULATED VALUES CAN BE USED TO CHECK THE C... PERFORMANCE OF BESSI0 TO BESSK1 [SPIEGEL (1991), PP. 246- C... 247] C... C... I0(1) = 1.266, I0(5) = 27.24 C... C... I1(1) = 0.5652, I1(5) = 24.34 C... C... K0(1) = 0.4210, K0(5) = 0.003691 C... C... K1(1) = 0.6019, K1(5) = 0.004045 C... C... SPIEGEL, M. R., (1991), MATHEMATICAL HANDBOOK OF FORMULAS AND C... TABLES, MCGRAW-HILL, NEW YORK C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... STEP THROUGH A SERIES OF VALUES OF X DO 1 I=1,10 X=0.5D0*DFLOAT(I) YI0=BESSI0(X) YI1=BESSI1(X) YK0=BESSK0(X) YK1=BESSK1(X) C... C... WRITE THE COMPUTED BESSEL FUNCTIONS WRITE(NO,2)X, YI0, YI1, YK0, YK1 2 FORMAT(' X = ',F5.2,' YI0 = ',F8.4,' YI1 = ',F8.4, + ' YK0 = ',F8.4,' YK1 = ',F8.4,/) C... C... NEXT X 1 CONTINUE C... C... END OF CALCULATION END DOUBLE PRECISION FUNCTION BESSI0(X) C... C... FUNCTION BESSI0 COMPUTES THE BESSEL FUNCTION I0(X). BESSI0 IS C... TAKEN FROM PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY AND C... W. T. VETTERLING (1985), NUMERICAL RECIPIES, CAMBRIDGE UNIVERSITY C... PRESS, CAMBRIDGE, CHAPTER 6, WITH CONVERSION TO DOUBLE PRECISION C... AND EDITING BY W. E. SCHIESSER. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA P1, P2, P3, + P4, P5, P6, + P7 + /1.0D0, 3.5156229D0, 3.0899424D0, + 1.2067492D0, 0.2659732D0, 0.360768D-1, + 0.45813D-2/ DATA Q1, Q2, Q3, + Q4, Q5, Q6, + Q7, Q8, Q9 + /0.39894228D0, 0.1328592D-1, 0.225319D-2, + -0.157565D-2, 0.916281D-2, -0.2057706D-1, + 0.2635537D-1, -0.1647633D-1, 0.392377D-2/ IF(DABS(X).LT.3.75D0)THEN Y=(X/3.75D0)**2 BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))) ELSE AX=DABS(X) Y=3.75D0/AX BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 + +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) END IF RETURN END DOUBLE PRECISION FUNCTION BESSI1(X) C... C... FUNCTION BESSI1 COMPUTES THE BESSEL FUNCTION I1(X). BESSI1 IS C... TAKEN FROM PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY AND C... W. T. VETTERLING (1985), NUMERICAL RECIPIES, CAMBRIDGE UNIVERSITY C... PRESS, CAMBRIDGE, CHAPTER 6, WITH CONVERSION TO DOUBLE PRECISION C... AND EDITING BY W. E. SCHIESSER. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA P1, P2, P3, + P4, P5, P6, + P7 + /0.5D0, 0.87890594D0, 0.51498869D0, + 0.15084934D0, 0.2658733D-1, 0.301532D-2, + 0.32411D-3/ DATA Q1, Q2, Q3, + Q4, Q5, Q6, + Q7, Q8, Q9 + /0.39894228D0, -0.3988024D-1, -0.362018D-2, + 0.163801D-2, -0.1031555D-1, 0.2282967D-1, + -0.2895312D-1, 0.1787654D-1, -0.420059D-2/ IF(DABS(X).LT.3.75D0)THEN Y=(X/3.75D0)**2 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE AX=DABS(X) Y=3.75D0/AX BESSI1=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+ + Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) END IF RETURN END DOUBLE PRECISION FUNCTION BESSK0(X) C... C... FUNCTION BESSK0 COMPUTES THE BESSEL FUNCTION K0(X). BESSK0 IS C... TAKEN FROM PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY AND C... W. T. VETTERLING (1985), NUMERICAL RECIPIES, CAMBRIDGE UNIVERSITY C... PRESS, CAMBRIDGE, CHAPTER 6, WITH CONVERSION TO DOUBLE PRECISION C... AND EDITING BY W. E. SCHIESSER. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA P1, P2, P3, + P4, P5, P6, + P7 + /-0.57721566D0, 0.42278420D0, 0.23069756D0, + 0.3488590D-1, 0.262698D-2, 0.10750D-3, + 0.74D-5/ DATA Q1, Q2, Q3, + Q4, Q5, Q6, + Q7 + /1.25331414D0, -0.7832358D-1, 0.2189568D-1, + -0.1062446D-1, 0.587872D-2, -0.251540D-2, + 0.53208D-3/ IF(X.LE.2.0D0)THEN Y=X*X/4.0D0 BESSK0=(-DLOG(X/2.0D0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+ + Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=(2.0D0/X) BESSK0=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ + Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) END IF RETURN END DOUBLE PRECISION FUNCTION BESSK1(X) C... C... FUNCTION BESSK1 COMPUTES THE BESSEL FUNCTION K1(X). BESSK1 IS C... TAKEN FROM PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY AND C... W. T. VETTERLING (1985), NUMERICAL RECIPIES, CAMBRIDGE UNIVERSITY C... PRESS, CAMBRIDGE, CHAPTER 6, WITH CONVERSION TO DOUBLE PRECISION C... AND EDITING BY W. E. SCHIESSER. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA P1, P2, P3, + P4, P5, P6, + P7 + /1.0D0, 0.15443144D0, -0.67278579D0, + -0.18156897D0, -0.1919402D-1, -0.110404D-2, + -0.4686D-4/ DATA Q1, Q2, Q3, + Q4, Q5, Q6, + Q7 + /1.25331414D0, 0.23498619D0, -0.3655620D-1, + 0.1504268D-1, -0.780353D-2, 0.325614D-2, + -0.68245D-3/ IF(X.LE.2.0D0)THEN Y=X*X/4.0D0 BESSK1=(DLOG(X/2.0D0)*BESSI1(X))+(1.0/X)*(P1+Y*(P2+ + Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=2.0D0/X BESSK1=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ + Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))) END IF RETURN END DOUBLE PRECISION FUNCTION DFLOAT(I) C... C... FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO A C... DOUBLE PRECISION FLOATING POINT. THIS FUNCTION IS PROVIDED C... IN CASE THE USER'S FORTRAN COMPILER DOES NOT INCLUDE DFLOAT. C... DFLOAT=DBLE(FLOAT(I)) RETURN END