SUBROUTINE INITAL IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NR=21) COMMON/T/ Z, NSTOP, NORUN + /Y/ U(NR) + /F/ UZ(NR) + /S/ UR(NR), URR(NR), V(NR), PE, + R0, DR, DRS, R(NR), + ZL + /I/ IP C... C... C... PROBLEM PARAMETERS C... C... PECLET NUMBER IF(NORUN.EQ.1)PE= 60.0D0 IF(NORUN.EQ.2)PE=120.0D0 IF(NORUN.EQ.3)PE=180.0D0 C... C... RADIUS R0=1.0D0 C... C... LENGTH ZL=30.0D0 C... C... RADIAL GRID, VELOCITY PROFILE DR=R0/DFLOAT(NR-1) DRS=DR**2 DO 1 I=1,NR R(I)=DR*DFLOAT(I-1) V(I)=(1.0D0-R(I)**2) 1 CONTINUE C... C... INITIAL CONDITION DO 2 I=1,NR U(I)=0.0D0 2 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/ Z, NSTOP, NORUN + /Y/ U(NR) + /F/ UZ(NR) + /S/ UR(NR), URR(NR), V(NR), PE, + R0, DR, DRS, R(NR), + ZL + /I/ IP C... C... BOUNDARY CONDITION AT R = R0 = 1 U(NR)=1.0D0 UZ(NR)=0.0D0 NU=1 C... C... UR CALL DSS004(0.0D0,R0,NR,U,UR) C... C... BOUNDARY CONDITION AT R = 0 UR(1)=0.0D0 NL=2 C... C... URR CALL DSS044(0.0D0,R0,NR,U,UR,URR,NL,NU) C... C... PDE DO 1 I=1,NR-1 C... C... R = 0 IF(I.EQ.1)THEN UZ(I)=1.0D0/(PE*2.0D0*V(I))*2.0D0*URR(I) C... C... R NE 0 ELSE IF(I.NE.1)THEN UZ(I)=1.0D0/(PE*2.0D0*V(I))*(URR(I)+(1.0D0/R(I))*UR(I)) END IF 1 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NR=21) COMMON/T/ Z, NSTOP, NORUN + /Y/ U(NR) + /F/ UZ(NR) + /S/ UR(NR), URR(NR), V(NR), PE, + R0, DR, DRS, R(NR), + ZL + /I/ IP C... C... MONITOR THE OUTPUT WRITE(*,*)NORUN,Z C... C... OPEN A MATLAB FILE FOR PLOTTING IP=IP+1 IF(IP.EQ.1)THEN OPEN(1,FILE='h5a.out') OPEN(2,FILE='h5b.out') END IF C... C... WRITE THE NUMERICAL AND AND ANALYTICAL SOLUTIONS AT SELECTED Z IF(((IP-1)*(IP-51)*(IP-101)).EQ.0)THEN C... C... NUMERICAL SOLUTION FOR 0 LE R LE R0 NRI=4 WRITE(NO,1)Z,(U(I),I=1,NR,NRI), + (UZ(I),I=1,NR,NRI) 1 FORMAT(//,' z = ',F4.1,/, + 12X,'r=0',4X,'r=0.2',4X,'r=0.4', + 4X,'r=0.6',4X,'r=0.8',4X,'r=1.0', + /,' u(r,z)',6F8.5, + /,' uz(r,z)',6F8.5,/) C... C... SERIES SOLUTION AT R = 0 CALL SERIES(PE,Z,UA) WRITE(NO,5)PE,Z,UA 5 FORMAT(/,' Series solution, Ua(0,z)',//, + 2X,' Pe = ',F5.0 ,/, + 2X,' z = ',F5.1 ,/, + 2X,' Ua = ',F8.5 ,//) END IF C... C... WRITE THE SOLUTION FOR MATLAB PLOTTING WRITE(1,2)Z,(U(I),I=1,NR) 2 FORMAT(22F8.4) C... C... COMPUTE THE SERIES SOLUTION TO BE SUPERIMPOSED ON THE NUMERICAL C... SOLUTION (AT THE END OF THE THIRD SOLUTION) IF((IP.EQ.101).AND.(NORUN.EQ.3))THEN DO 3 NCASE=1,3 IF(NCASE.EQ.1)PE= 60.0D0 IF(NCASE.EQ.2)PE=120.0D0 IF(NCASE.EQ.3)PE=180.0D0 DO 4 I=1,101 ZA=ZL*DFLOAT(I-1)/DFLOAT(101-1) CALL SERIES(PE,ZA,UA) WRITE(2,2)ZA,UA 4 CONTINUE 3 CONTINUE END IF RETURN END SUBROUTINE SERIES(PE,Z,U) C... C... The Graetz problem stated by Jakob (1949) is C... C... 2*vavg*v(r)*ua = D*(ua + (1/r)*ua ) (1) C... z rr r C... C... v(r) = 1 - r**(2) C... C... where ua = (T - Ts)/(To - Ts) C... C... For the numerical solution, u = (T - To)/(Ts - To). u can be C... expressed in terms of ua. Thus, C... C... T = Ts + ua*(To - Ts) C... C... u = (Ts + ua*(To - Ts) - To)/(Ts - To) = 1 - ua C... C... The analytical (series) solution given by Jakob, eq. (22-45), C... is C... C... ua = 1.477*exp(-3.658*(1/Pe)*z)*R0(r) C... C... - 0.810*exp(22.178*(1/Pe)*z)*R1(r) C... C... + 0.385*exp(-53.05*(1/Pe)*z)*R2(r) C... C... where z = z/ro, r = r/ro, (the LHS z and r are dimensionaless), C... R0(0) = R1(0) = R2(0) = 1 (from Jakob, Table,22-1, page 455), C... Pe = ro*vavg/D (D = k/(rho*Cp)) C... C... Jakob, M., Heat Transfer, John Wiley and Sons, Inc. New York, C... Vol. 1, 1949 C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... RATIO OF Z/PE USED IN SERIES SOLUTION ZPE=Z/PE C... C... UA(0,Z) TERM1= 1.477D0*DEXP( -3.658D0*ZPE) TERM2=-0.810D0*DEXP(-22.178D0*ZPE) TERM3= 0.385D0*DEXP( -53.05D0*ZPE) UA=TERM1+TERM2+TERM3 C... C... FOR COMPARISON WITH THE NUMERICAL SOLUTION, U = 1 - UA U=1.0D0-UA C... C... END OF CALCULATION RETURN END