SUBROUTINE INITAL C... C... H3 - HEAT CONDUCTION IN A CIRCULAR FIN C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(NR=21) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR) + /F/ UT(NR) + /P/ ALPHA, BETA, UA, U0, DR, DRS, + R(NR) + /I/ IP 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=1.0D0 C... C... RHO*CP RHOCP=2.0D0 C... C... THERMAL DIFFUSIVITY ALPHA=0.01D0 C... C... INITIAL TEMPERATURE UI=25.0D0 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, SQUARE OF INCREMENT DR=(R1-R0)/DFLOAT(NR-1) DRS=DR**2 C... C... RADIAL GRID AND INITIAL CONDITION DO 1 I=1,NR R(I)=R0+DFLOAT(I-1)*DR U(I)=UI 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(NR=21) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR) + /F/ UT(NR) + /P/ ALPHA, BETA, UA, U0, DR, DRS, + R(NR) + /I/ IP C... C... BOUNDARY CONDITION AT R = R0 U(1)=U0 UT(1)=0.0D0 C... C... BOUNDARY CONDITION AT R = R1 C... C... RADIAL GROUP NOT INCLUDED IF(NORUN.EQ.1)THEN UT(NR)=-BETA*(U(NR)-UA) C... C... RADIAL GROUP INCLUDED ELSE IF(NORUN.EQ.2)THEN UT(NR)=ALPHA*2.0D0*(U(NR-1)-U(NR))/DRS-BETA*(U(NR)-UA) END IF C... C... INTERIOR POINTS DO 1 I=2,NR-1 UT(I)=ALPHA*(U(I+1)-2.0D0*U(I)+U(I-1))/DRS + +ALPHA*(1.0D0/R(I))*(U(I+1)-U(I-1))/(2.0D0*DR) + -BETA*(U(NR)-UA) 1 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(NR=21) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR) + /F/ UT(NR) + /P/ ALPHA, BETA, UA, U0, DR, DRS, + R(NR) + /I/ IP C... C... MONITOR OUTPUT WRITE(*,*)' NORUN = ',NORUN,' T = ',T C... C... IF THE SOLUTION HAS A PHYSICALLY UNREALISTIC VALUE, TERMINATE C... EXECUTION OF CURRENT RUN DO 2 I=1,NR IF((U(I).GT.U0).OR.(U(I).LT.UA))THEN WRITE(NO,1)NORUN,T 1 FORMAT(//,' RUN NO = ',I2,' AT T = ',F7.2,' TERMINATED',//) NSTOP=1 END IF 2 CONTINUE C... C... PRINT THE SOLUTION WRITE(NO,3)T,(R(I),I=1,NR,4),(U(I),I=1,NR,4) 3 FORMAT(' t = ',F6.2,/,' r ',6F9.2,/,' u(r,t) ',6F9.3,/) RETURN END