SUBROUTINE INITAL C... C... H4 - DYNAMIC GRAETZ PROBLEM WITH CONSTANT WALL HEAT FLUX C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... SELECT APPROXIMATION OF SPATIAL DERIVATIVES (SEE THE COMMENTS C... IN SUBROUTINE DERV FOR AN EXPLANATION OF THE TWO CASES) NCASE=2 C... C... PROBLEM PARAMETERS C... C... PECLET NUMBER PE=60.0D0 C... C... QW*R0/K = TN (NORMALIZATION TEMPERATURE) TN=1.0D0 C... C... OUTSIDE RADIUS R0=1.0D0 C... C... AXIAL LENGTH IF(NORUN.EQ.1)ZL=30.0D0 IF(NORUN.EQ.2)ZL=90.0D0 C... C... RADIAL GRID, VELOCITY PROFILE DR=R0/DFLOAT(NR-1) DO 1 IR=1,NR R(IR)=DFLOAT(IR-1)*DR V(IR)=2.0D0*(1.0D0-R(IR)**2) 1 CONTINUE DRS=DR**2 C... C... AXIAL GRID DZ=ZL/DFLOAT(NZ-1) DO 2 JZ=1,NZ Z(JZ)=DFLOAT(JZ-1)*DZ 2 CONTINUE DZS=DZ**2 C... C... INITIAL CONDITION DO 3 IR=1,NR DO 3 JZ=1,NZ U(IR,JZ)=1.0D0 3 CONTINUE C... C... GIVE THE TIME DERIVATIVES A VALUE WHICH CAN THEN BE USED TO C... TEST FOR UNDEFINED DERIVATIVES BY SETTING THEM TO PI, THEN C... CALLING DERV, AND CHECKING IN SUBROUTINE PRINT IF THEY HAVE C... ALL BEEN SET TO VALUES OTHER THAN PI. THIS IS A SIMPLE PRO- C... CEDURE TO ENSURE THAT ALL NR X NZ DERIVATIVES HAVE BEEN SET C... IN DERV PI=3.1415927D0 DO 4 IR=1,NR DO 4 JZ=1,NZ UT(IR,JZ)=PI 4 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/I/ IP, NCASE C... C... APPROXIMATIONS OF SPATIAL DERIVATIVES C... C... NCASE = 1 - EXPLICIT FINITE DIFFERENCES (FDM) IF(NCASE.EQ.1)THEN CALL DERV1 C... C... NCASE = 2 - FINITE DIFFERENCES FROM DSS034 ELSE + IF(NCASE.EQ.2)THEN CALL DERV2 END IF RETURN END SUBROUTINE DERV1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=1.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ UNRP1=U(NR-1,JZ)+(2.0D0*DR) UT(NR,JZ)=-PE*0.0D0 + +(UNRP1-2.0D0*U(NR,JZ)+U(NR-1,JZ))/DRS + +(UNRP1-U(NR-1,JZ))/(2.0D0*DR) 2 CONTINUE C... C... R = 0, Z NE 0 DO 3 JZ=2,NZ UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ + +4.0D0*(U(2,JZ)-U(1,JZ))/DRS 3 CONTINUE C... C... R NE 0, R0, Z NE 0 DO 4 JZ=2,NZ DO 4 IR=2,NR-1 UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ + +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS + +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR) 4 CONTINUE RETURN END SUBROUTINE DERV2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR UR, URR, UZ DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ) C... C... Z = 0 DO 1 IR=1,NR U(IR,1)=1.0D0 UT(IR,1)=0.0D0 1 CONTINUE C... C... UR CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0) C... C... R = R0, Z NE 0 DO 2 JZ=2,NZ UR(NR,JZ)=1.0D0 2 CONTINUE C... C... R = 0, Z NE 0 DO 3 JZ=2,NZ UR(1,JZ)=0.0D0 3 CONTINUE C... C... URR CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0) C... C... UZ (FOR CONVECTION) CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0) C... C... Z NE 0 DO 4 JZ=2,NZ DO 4 IR=1,NR C... C... R = 0 IF(IR.EQ.1)THEN UT(1,JZ)=-PE*V(1)*UZ(1,JZ) + +2.0D0*URR(1,JZ) C... C... R NE 0 ELSE + IF(IR.NE.1)THEN UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ) + +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ) END IF 4 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... ARRAYS FOR THE ANALYTICAL SOLUTION DIMENSION UZ(NR), UA(NR) C... C... MONITOR SOLUTION WRITE(*,*)NCASE,NORUN,T C... C... PRINT THE INITIAL CONDITION AND INITIAL TEMPORAL DERIVATIVE IF(IP.EQ.0)THEN WRITE(NO,10)T,((U(IR,JZ),IR=1,NR),JZ=1,NZ) 10 FORMAT(' t = ',F6.2,//,' u(r,z,t)',/,(11F7.3)) WRITE(NO,11)(( UT(IR,JZ),IR=1,NR),JZ=1,NZ) 11 FORMAT(//,' ut(r,z,t)',/,(11F7.1)) C... C... CHECK FOR UNDEFINED DERIVATIVES DO 12 IR=1,NR DO 12 JZ=1,NZ IF((UT(IR,JZ).GT.3.14D0).AND.(UT(IR,JZ).LT.3.15D0))THEN WRITE (NO,13)IR,JZ,U(IR,JZ),UT(IR,JZ) 13 FORMAT(' IR = ',I3,' JZ = ',I3, + ' U(IR,JZ) = ',F7.3,' UT(IR,JZ) = ',F7.1) END IF 12 CONTINUE END IF C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,1)T 1 FORMAT(//,' t = ',F6.2,//,29X,'u(r,z,t)',//, + 3X,'z',4X,'r=0.0',4X,'r=0.2',4X,'r=0.4', + 4X,'r=0.6',4X,'r=0.8',4X,'r=1.0',/) DO 3 JZ=1,NZ NRI=2 WRITE(NO,2)Z(JZ),(U(IR,JZ),IR=1,NR,NRI) 2 FORMAT(F4.1,6F9.4) 3 CONTINUE C... C... AT THE END OF THE RUN, CALCULATE AND PRINT THE ANALYTICAL, C... STEADY STATE TEMPERATURE IP=IP+1 IF(IP.EQ.11)THEN DO 14 IR=1,NR UZ(IR)=(U(IR,NZ)-U(IR,NZ-1))/DZ UA(IR)=U(NR,NZ)+2.0D0*PE*TN*UZ(IR)* + ((1.0D0/ 4.0D0)*(R(IR)**2) + -(1.0D0/16.0D0)*(R(IR)**4)-3.0D0/16.0D0) 14 CONTINUE WRITE(NO,15)(UZ(IR),IR=1,NR,NRI),(UA(IR),IR=1,NR,NRI) 15 FORMAT(/,10X,'Axial gradient; analytical, steady state,', + /,16X,'thermally developed solution',//, + ' uz',6F9.4,/,' ua',6F9.4,/) END IF C... C... STORE THE SOLUTION AND PLOT THE SOLUTION AT THE END OF THE C... SECOND RUN C... C... PLOT RADIAL PROFILES CALL PLOTDR C... C... PLOT CENTERLINE AXIAL PROFILES CALL PLOTDZ RETURN END SUBROUTINE PLOTDR IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... OPEN A FILE FOR MATLAB PLOTTING AT THE BEGINNING OF THE C... FIRST RUN IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN OPEN(1,FILE='h4r.out') END IF C... C... WRITE THE TEMPERATURES, U(R,ZL,T), FOR SUBSEQUENT MATLAB C... PLOTTING DO 2 IR=1,NR WRITE(1,3)R(IR),U(IR,NZ) 3 FORMAT(F10.3,F10.5) 2 CONTINUE RETURN END SUBROUTINE PLOTDZ IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=11,NZ=31) COMMON/T/ T, NSTOP, NORUN + /Y/ U(NR,NZ) + /F/ UT(NR,NZ) + /G/ R(NR), Z(NZ), V(NR) + /C/ PE, TN, R0, ZL, + DR, DRS, DZ, DZS + /I/ IP, NCASE C... C... OPEN A FILE FOR MATLAB PLOTTING AT THE BEGINNING OF THE C... FIRST RUN IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN OPEN(2,FILE='h4z.out') END IF C... C... WRITE THE CENTERLINE TEMPERATURE, U(0,Z,T), FOR SUBSEQUENT C... MATLAB PLOTTING DO 2 JZ=1,NZ WRITE(2,3)Z(JZ),U(1,JZ) 3 FORMAT(F10.3,F10.5) 2 CONTINUE RETURN END