* This is a listing of all the Fortran programs discussed in the * textbook "Introduction to Convective Heat Transfer Analysis" ************************************************************************ ***** ***** ***** DEVPIPE ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR LAMINAR THERMALLY DEVELOPING FLOW ***** ***** ***** ***** IN A PIPE USING THE FINITE DIFFERENCE TECHNIQUE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ***** __ _ ***** ************************************************************************ * DIMENSION T(2,150),A(150),B(150) DIMENSION C(150),D(150),R(150),H(150),U(150) REAL ND,NDM * ********************************************************************** * OPEN(1,FILE='DEVPIPR.DAT') OPEN(2,FILE='DEVPIPL.DAT') * ********************************************************************** * * ************************ SPECIFY KNOWN VALUES *********************** * WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) IWALL * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) IF (IWALL.EQ.1) THEN WRITE (1,1002) ELSE WRITE (1,1003) ENDIF * ************************** ASSIGN KNOWN VALUES *********************** * Z = 0.0 ZMAX=0.7 N = 91 N1=N-1 R(1) = 0.0 U(1)=2.0 DR =0.5/(N-1) * DO 1 J = 2,N R(J) = R(J-1) + DR U(J)=2.0*(1.0-R(J)*R(J)/0.25) 1 CONTINUE * *********************** ASSIGN INITIAL VALUES *********************** * IF (IWALL.EQ.1) THEN DO 2 J = 1,N T(1,J) = 1.0 2 CONTINUE T(1,N) = 0.0 ELSE DO 3 J = 1,N T(1,J) = 0.0 3 CONTINUE T(1,N) = DR ENDIF DZ = 0.0001 DZMAX=0.005 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * IF (DZ.LT.DZMAX) THEN DZ = 1.05*DZ ELSE DZ = DZMAX ENDIF * Z = Z + DZ WRITE (1,2000) Z * * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * DO 12 J = 2,N1 A(J)= U(J)/DZ+(R(J+1)+R(J))/(2.0*R(J)*DR*DR)+ & (R(J)+R(J-1))/(2.0*R(J)*DR*DR) B(J) = -(R(J+1)+R(J))/(2.0*R(J)*DR*DR) C(J) = -(R(J)+R(J-1))/(2.0*R(J)*DR*DR) D(J)= (U(J)/DZ)*T(1,J) 12 CONTINUE * A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 IF (IWALL.EQ.1) THEN A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=0.0 ELSE A(N)=1.0 B(N)=0.0 C(N)=-1.0 D(N)=DR ENDIF * CALL TRISOL(N,A,B,C,D,H) * * DO 13 J = 1,N T(2,J) = H(J) 13 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N T(1,J)=T(2,J) 30 CONTINUE * ************************** WRITE THE RESULTS ************************** * IF (IWALL.EQ.1) THEN ND = (T(2,N-1)-T(2,N))/DR ELSE ND=1.0/T(2,N) ENDIF SUM=0.0 DO 35 J=2,N SUM=0.5*(U(J)*T(2,J)*R(J)+U(J-1)*T(2,J-1)*R(J-1))*DR+SUM 35 CONTINUE IF (IWALL.EQ.1) THEN NDM=ND/(8.0*SUM) ELSE NDM=1.0/(T(2,N)-8.0*SUM) ENDIF WRITE(1,4000) Z,ND,NDM,T(2,1) WRITE(6,4200) Z,ND,NDM,T(2,1) WRITE(2,4100) Z,ND,NDM,T(2,1) WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,R(J),U(J),T(2,J) 40 CONTINUE * ************************** CHECK FOR MAXIMUM Z ************************ * IF(Z.GE.ZMAX) GO TO 102 GO TO 100 102 WRITE(1,5001) WRITE(6,5001) CLOSE(1) CLOSE(2) STOP * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'THERMALLY DEVELOPING FLOW IN A PIPE') 1001 FORMAT (10X,'-----------------------------------',//) 1002 FORMAT (14X,'UNIFORM WALL TEMPERATURE',//) 1003 FORMAT (15X,'UNIFORM WALL HEAT FLUX',//) 2000 FORMAT (/,5X,'SOLUTION FOR Z =',F12.4,/ C ,5X,'___________________________',/) 4000 FORMAT (5X,'Z =',F10.6,'ND =',F10.5,'NDM =',F10.5,'TCEN =',F10.5) 4001 FORMAT (2X,'J',7X,'R',12X,'U',12X,'T') 4002 FORMAT (I4,3F12.4) 4100 FORMAT (F10.5,',',F10.5,',',F10.5,',',F10.5) 4200 FORMAT (5X,'Z =',F10.6,'ND =',F10.5,'NDM =',F10.5,'TCEN =',F10.5) 5001 FORMAT (/,5X,'********** ZMAX REACHED **********') 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUE OF THE PARAMETER THAT DETERMINES', & /,' THE WALL THERMAL BOUNDARY CONDITION', & /,' - 1 FOR A UNIFORM TEMPERATURE, 2 FOR A UNIFORM HEAT FLUX') END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(150),B(150),C(150),D(150),H(150),W(150),Q(150),G(150) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END ************************************************************************ ***** ***** ***** DEVDUCT ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR LAMINAR THERMALLY DEVELOPING FLOW ***** ***** ***** ***** IN A PLANE DUCT USING THE FINITE DIFFERENCE TECHNIQUE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ***** __ _ ***** ************************************************************************ * DIMENSION T(2,150),A(150),B(150) DIMENSION C(150),D(150),Y(150),H(150),U(150) REAL NW,NWM * ********************************************************************** * OPEN(1,FILE='DEVDUPR.DAT') OPEN(2,FILE='DEVDUPL.DAT') OPEN(3,FILE='DEVDUIN.DAT') * ********************************************************************** * *********************** SPECIFY KNOWN VALUES *********************** * * WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) IWALL * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) IF (IWALL.EQ.1) THEN WRITE (1,1002) ELSE WRITE (1,1003) ENDIF * ************************** ASSIGN KNOWN VALUES *********************** * Z = 0.0 ZMAX=0.7 N = 91 N1=N-1 Y(1) = 0.0 U(1)=1.5 DY =0.5/(N-1) * DO 1 J = 2,N Y(J) = Y(J-1) + DY U(J)=1.5*(1.0-Y(J)*Y(J)/0.25) 1 CONTINUE * *********************** ASSIGN INITIAL VALUES *********************** * IF (IWALL.EQ.1) THEN DO 2 J = 1,N T(1,J) = 1.0 2 CONTINUE T(1,N) = 0.0 ELSE DO 3 J = 1,N T(1,J) = 0.0 3 CONTINUE T(1,N) = DY ENDIF DZ = 0.0001 DZMAX=0.005 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * IF (DZ.LT.DZMAX) THEN DZ = 1.05*DZ ELSE DZ = DZMAX ENDIF * Z = Z + DZ WRITE (1,2000) Z * * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * DO 12 J = 2,N1 A(J)= U(J)/DZ+2.0/(DY*DY) B(J) = -1.0/(DY*DY) C(J) = -1.0/(DY*DY) D(J)= (U(J)/DZ)*T(1,J) 12 CONTINUE * A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 IF (IWALL.EQ.1) THEN A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=0.0 ELSE A(N)=1.0 B(N)=0.0 C(N)=-1.0 D(N)=DY ENDIF * CALL TRISOL(N,A,B,C,D,H) * * DO 13 J = 1,N T(2,J) = H(J) 13 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N T(1,J)=T(2,J) 30 CONTINUE * ************************** WRITE THE RESULTS ************************** * IF (IWALL.EQ.1) THEN NW = (T(2,N-1)-T(2,N))/DY ELSE NW=1.0/T(2,N) ENDIF SUM=0.0 DO 35 J=2,N SUM=0.5*(U(J)*T(2,J)+U(J-1)*T(2,J-1))*DY+SUM 35 CONTINUE IF (IWALL.EQ.1) THEN NWM=NW/(2.0*SUM) ELSE NWM=1.0/(T(2,N)-2.0*SUM) ENDIF WRITE(1,4000) Z,NW,NWM,T(2,1) WRITE(6,4200) Z,NW,NWM,T(2,1) WRITE(2,4100) Z,NW,NWM,T(2,1) WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,Y(J),U(J),T(2,J) 40 CONTINUE * ************************** CHECK FOR MAXIMUM Z ************************ * IF(Z.GE.ZMAX) GO TO 102 GO TO 100 102 WRITE(1,5001) WRITE(6,5001) CLOSE(1) CLOSE(2) CLOSE(3) STOP * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'THERMALLY DEVELOPING FLOW IN A PLANE DUCT') 1001 FORMAT (10X,'-----------------------------------------',//) 1002 FORMAT (14X,'UNIFORM WALL TEMPERATURE',//) 1003 FORMAT (15X,'UNIFORM WALL HEAT FLUX',//) 2000 FORMAT (/,5X,'SOLUTION FOR Z =',F12.4,/ C ,5X,'___________________________',/) 4000 FORMAT (5X,'Z =',F10.6,'NW =',F10.5,'NWM =',F10.5,'TCEN =',F10.5) 4001 FORMAT (2X,'J',7X,'Y',12X,'U',12X,'T') 4002 FORMAT (I4,3F12.4) 4100 FORMAT (F10.5,',',F10.5,',',F10.5,',',F10.5) 4200 FORMAT (5X,'Z =',F10.6,'NW =',F10.5,'NWM =',F10.5,'TCEN =',F10.5) 5001 FORMAT (/,5X,'********** ZMAX REACHED **********') 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUE OF THE PARAMETER THAT DETERMINES', & /,' THE WALL THERMAL BOUNDARY CONDITION', & /,' - 1 FOR A UNIFORM TEMPERATURE, 2 FOR A UNIFORM HEAT FLUX') END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(150),B(150),C(150),D(150),H(150),W(150),Q(150),G(150) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END ************************************************************************ ***** ***** ***** DUCTSYM ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR DEVELOPING LAMINAR PLANE DUCT FLOW ***** ***** ***** ***** THE VELOCITY AND TEMP. FIELDS ARE SIMULTANEOUSLY DEVELOPING ***** ***** ***** ***** THE FLOW IS ASSUMED TO BE SYMMETRICAL ABOUT THE CENTRE-LINE ***** ***** ***** ***** AND THE FLOW IN HALF OF THE DUCT IS CALCULATED ***** ***** ***** ***** AN EXPLICIT FINITE DIFFERENCE TECHNIQUE IS USED ***** ***** ***** ***** EITHER THE WALL TEMP. OR THE WALL HEAT FLUX IS UNIFORM ***** ***** ***** ************************************************************************ * REAL U(2,100),V(2,100),T(2,100),A(100),B(100),C(100) REAL PR,ZMAX,DY,DZ * *********************************************************************** * OPEN(UNIT=1,FILE='DUSTPR.DAT') OPEN(UNIT=2,FILE='DUSTPL.DAT') * *********************************************************************** * WRITE(1,760) WRITE(6,760) * *********************** ENTER SPECIFIED VALUES ********************** * write(6,1400) write(6,1070) read(5,*) ZMAX write(6,1060) read(5,*) PR write(6,1050) read(5,*) N write(6,1600) read(5,*) IQT * ********************************************************************* IPRF=50 IPPF=100 * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * IPRF=PRINT FREQUENCY * IPPF=PROFILE PRINT FREQUENCY * ZMAX=MAXIMUM Z TO WHICH SOLUTION IS TO BE UNDER TAKEN * IQT=BOUNDARY CONDITION IDENTIFIER (1=UNIFORM TEMP., 2=UNIFORM FLUX) * IF (IQT.EQ.1) THEN WRITE (1,1100) ELSE WRITE (1,1200) ENDIF WRITE (1,1000) N,PR,IPRF,IPPF,ZMAX * ********************************************************************** * IPRT=0 IPPT=0 N1=N-1 DY=0.5/N1 ZFAC=10.0*PR * ************************** ASSIGN INITIAL VALUES ********************* * Z=0.0 P1=-0.5 DO 10 J=1,N U(1,J)=1.0/(1.0-DY) V(1,J)=0.0 T(1,J)=0.0 10 CONTINUE U(1,N)=0.0 IF (IQT.EQ.1) THEN T(1,N)=1.0 ELSE T(1,N)=DY ENDIF * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * ***************** FIND Z-STEP TO ENSURE STABILITY ******************* * DZ=DY*DY*U(1,N1)/ZFAC * ****************** FIND COEFFICIENTS ******************* * A(1)=P1/U(1,1)+(PR*DZ/U(1,1))*(2.0*U(1,2)-2.0*U(1,1))/(DY*DY) DO 20 J=2,N1 A(J)=P1/U(1,J)+(PR*DZ/U(1,J))* & (-V(1,J)*(U(1,J+1)-U(1,J-1))/(2.0*DY*PR) & +(U(1,J+1)+U(1,J-1)-2.0*U(1,J))/(DY*DY)) 20 CONTINUE DO 25 J=2,N1 B(J)=(DY/(2.0*DZ))*(1.0/U(1,J)+1.0/U(1,J-1)) C(J)=(DY/(2.0*DZ))*(A(J)+A(J-1)) 25 CONTINUE B(N)=(DY/(2.0*DZ))*(1.0/U(1,N-1)) C(N)=(DY/(2.0*DZ))*A(N-1) * **************** CALCULATE THE PRESSURE ***************************** * BSUM=0.0 CSUM=0.0 DO 30 J=2,N BSUM=BSUM+B(J) CSUM=CSUM+C(J) 30 CONTINUE * P2=CSUM/BSUM * ********* CALCULATE STREAM WISE VELOCITY COMPONENT ******************* * DO 40 J=2,N1 U(2,J)=U(1,J)+A(J)-P2/U(1,J) 40 CONTINUE U(2,1)=U(2,2) U(2,N)=0.0 * ********* SOLVE FOR THE CROSS-STREAM VELOCITY COMPONENT ************ * V(2,1)=0.0 DO 50 J=2,N1 V(2,J)=V(2,J-1)+B(J)*P2-C(J) 50 CONTINUE V(2,N)=0.0 * ************* CALCULATE TEMPERATURE FROM ENERGY EQUATION ************ * DO 60 J=2,N1 T(2,J)=T(1,J)+(DZ/U(1,J))*(-V(1,J)*(T(1,J+1)-T(1,J-1))/(2.0*DY) $ +(T(1,J+1)+T(1,J-1)-2.0*T(1,J))/(DY*DY)) 60 CONTINUE T(2,1)=T(2,2) IF (IQT.EQ.1) THEN T(2,N)=1.0 ELSE T(2,N)=T(2,N1)+DY ENDIF * *********************** CALCULATE NUSSELT NUMBER ********************** * IF (IQT.EQ.1) THEN XNUW=(T(2,N)-T(2,N1))/DY ELSE XNUW=1.0/T(2,N) ENDIF XNUMW=0.0 XMASS=0.0 DO 70 J=2,N XNUMW=XNUMW+0.5*(U(2,J)*T(2,J)+U(2,J-1)*T(2,J-1))*DY XMASS=XMASS+0.5*(U(2,J)+U(2,J-1))*DY 70 CONTINUE IF (IQT.EQ.1) THEN XNUDM=XNUW/(1.0-XNUMW/XMASS) ELSE XNUDM=1.0/(T(2,N)-XNUMW/XMASS) ENDIF Z=Z+DZ XNUMW=XNUMW/Z * ********************* WRITE THE RESULTS ****************************** * IPRT=IPRT+1 IF(IPRT.LT.IPRF) GO TO 300 IPRT=0 IF (IQT.EQ.1) THEN WRITE(1,2000) Z,XNUW,XNUDM,XNUMW WRITE(2,2200) Z,XNUW,XNUDM,XNUMW WRITE(6,2000) Z,XNUW,XNUDM,XNUMW ELSE WRITE(1,2100) Z,XNUW,XNUDM WRITE(2,2300) Z,XNUW,XNUDM WRITE(6,2100) Z,XNUW,XNUDM ENDIF IPPT=IPPT+1 IF(IPPT.LT.IPPF) GO TO 300 IPPT=0 WRITE(1,3000) Y=0.0 DO 80 J=1,N WRITE(1,4000) Y,U(2,J),T(2,J),V(2,J) Y=Y+DY 80 CONTINUE 300 CONTINUE * ************************** CHECK FOR MAXIMUM X ************************ * IF(Z.GT.ZMAX) GO TO 200 * ***************************** RETURN VALUES *************************** * DO 90 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 90 CONTINUE P1=P2 GO TO 100 200 CONTINUE * ************************** FORMAT STATEMENTS ************************** * 760 FORMAT(//,3X,'DEVELOPING FLOW THROUGH IN A PLANE DUCT',//) 1000 FORMAT(6X,'NUMBER OF GRID POINTS = ',I4,/, $6X,'PRANDTL NUMBER = ' ,F10.5,/, $6X,'PRINT FREQUENCY = ',I3,/, $6X,'PROFILE FREQUENCY = ',I3,/, $6X,'MAXIMUM Z = ',F10.5,//) 1050 format(//,5X,'Input the Number of Nodal Points ',/) 1060 format(//,5X,'Input the Prandtl Number ',/) 1070 format(//,5X,'Input the Maximum Z-value to be Considered ',/) 1100 FORMAT(6X,'UNIFORM WALL TEMPERATURE') 1200 FORMAT(6X,'UNIFORM WALL HEAT FLUX') 1400 format(////////,5X,'DEVELOPING LAMINAR DUCT FLOW', & /,5X,'--------------------------',///) 1600 format(//,5X,'Uniform Temp. (1) or Uniform Heat Flux (2) ',/) 2000 FORMAT(2X,'Z =',F6.4,2X, 'NU (TW-TIN)=',F7.4, $ 2X,'NU (TW-TMEAN)=',F7.4,2X,'NU (AVG. TO Z)=',F7.4) 2100 FORMAT(2X,'Z =',F8.6,3X, 'NU (TW-TIN)=' ,F9.4, $ 3X, 'NU (TW-TMEAN)=' ,F9.4) 2200 FORMAT(F6.4,',',F7.4,',',F7.4,',',F7.4) 2300 FORMAT(F6.4,',',F7.4,',',F7.4) 3000 FORMAT(' Y U T V' ) 4000 FORMAT(4F12.4) * ************************** ENDING ************************************* * STOP CLOSE(1) CLOSE(2) END ************************************************************************ ***** ***** ***** ENCLREC ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR NATURAL CONVECTIVE ***** ***** ***** ***** LAMINAR FLOW IN A RECTANGULAR ***** ***** ***** ***** ENCLOSURE USING A FINITE DIFFERENCE TECHNIQUE WITH ***** ***** ***** ***** AN UNDER-RELAXATION ITERATIVE PROCEDURE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ************************************************************************ C DIMENSION PSI(61,61),T(61,61),OME(61,61) DIMENSION TOLD(61,61),NU(2,61) INTEGER CNT REAL LEN,NU C C *************************************************************** C OPEN (2,FILE='ENCLCON.DAT') OPEN (1,FILE='ENCLREC.DAT') C C ******************************************************************** C C **************SPECIFIED VALUES************************** C WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) AR,RA,PR WRITE(6,6002) READ(5,*) IEND C 1000 FORMAT (10X,'FLOW IN A RECTANGULAR ENCLOSURE') 1001 FORMAT (10X,'-------------------------------',//) 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUES OF THE ASPECT RATIO,', & /,' THE RAYLEIGH NUMBER AND THE PRANDTL NUMBER - A,RA,PR',///) 6002 FORMAT(/,' INPUT THE VALUE OF THE PARAMETER THAT DETERMINES', & /,' THE END-WALL THERMAL BOUNDARY CONDITION', & /,' 1 FOR ADIABATIC OR 2 FOR PERFECTLY CONDUCTING END-WALLS',//) C NX=37 NY=37 TOL=1.0E-5 IMAX=50000 IMIN=25000 IAVG=IMIN-5000 REX=0.10 REXB=0.15 REXT=0.05 C C ************INITIAL CONDITIONS************************* C HSUM=0.0 ITER=0 DO 100 I=1,NX DO 110 J=1,NY PSI(I,J)=0.0 T(I,J)=1.0-(I-1)/(NX-1) TOLD(I,J)=T(I,J) 110 OME(I,J)=0.0 100 CONTINUE C C ************BOUNDARY CONDITIONS************************* C DO 120 J=1,NY T(1,J)=1.0 T(NX,J)=0.0 TOLD(1,J)=1.0 TOLD(NX,J)=0.0 120 CONTINUE C C ***********CALCULATION BEGINS************************* C NX1=NX-1 NY1=NY-1 NY2=NY-2 DX=1.0/NX1 DY=AR/NY1 FAC=2*(1/DX**2+1/DY**2) 199 ITER=ITER+1 C C ************CALCULATION OF STREAM FUNCTION "PSI"******** C DO 200 IREP=1,2 DO 210 I=2,NX1 DO 220 J=2,NY1 F1=1/FAC*(OME(I,J)+(PSI(I+1,J)+PSI(I-1,J))/DX**2+ * (PSI(I,J+1)+PSI(I,J-1))/DY**2) 220 PSI(I,J)=PSI(I,J)+REX*(F1-PSI(I,J)) 210 CONTINUE 200 CONTINUE C C ************CALCULATION OF TEMPERATURE FIELD "T"******** C DO 300 IREP=1,2 DO 310 I=2,NX1 DO 320 J=2,NY1 F2=1/FAC*((T(I+1,J)+T(I-1,J))/DX**2+ * (T(I,J+1)+T(I,J-1))/DY**2+ * (PSI(I+1,J)-PSI(I-1,J))*(T(I,J+1)-T(I,J-1))/(4*DX*DY)- * (PSI(I,J+1)-PSI(I,J-1))*(T(I+1,J)-T(I-1,J))/(4*DX*DY)) 320 T(I,J)=T(I,J)+REX*(F2-T(I,J)) IF(IEND.EQ.1) THEN F2=(4.0/3.0)*(T(I,2)-T(I,3)/4.0) T(I,1)=T(I,1)+REXT*(F2-T(I,1)) F2=(4.0/3.0)*(T(I,NY1)-T(I,NY2)/4.0) T(I,NY)=T(I,NY)+REXT*(F2-T(I,NY)) ELSE T(I,1)=1.0-((I-1.0)/(NX-1.0)) T(I,NY)=1.0-((I-1.0)/(NX-1.0)) ENDIF 310 CONTINUE 300 CONTINUE C C *******CALCULATION OF VORTICITY "OME" FOR INTERNAL NODES***** C DO 400 IREP=1,2 DO 410 I=2,NX1 DO 420 J=2,NY1 F3=1/FAC*((OME(I+1,J)+OME(I-1,J))/DX**2+ * (OME(I,J+1)+OME(I,J-1))/DY**2+ * RA*(T(I+1,J)-T(I-1,J))/(2*DX)- * (PSI(I,J+1)-PSI(I,J-1))*(OME(I+1,J)-OME(I-1,J))/(4*PR*DX*DY)+ * (PSI(I+1,J)-PSI(I-1,J))*(OME(I,J+1)-OME(I,J-1))/(4*PR*DX*DY)) 420 OME(I,J)=OME(I,J)+REX*(F3-OME(I,J)) 410 CONTINUE 400 CONTINUE C C *******CALCULATION OF VORTICITY "OME" FOR BOUNDARY NODES***** C DO 510 J=1,NY F41=-2*PSI(2,J)/DX**2 OME(1,J)=OME(1,J)+REXB*(F41-OME(1,J)) F42=-2*PSI(NX1,J)/DX**2 510 OME(NX,J)=OME(NX,J)+REXB*(F42-OME(NX,J)) DO 520 I=1,NX F43=-2*PSI(I,2)/DY**2 OME(I,1)=OME(I,1)+REXB*(F43-OME(I,1)) F44=-2*PSI(I,NY1)/DY**2 520 OME(I,NY)=OME(I,NY)+REXB*(F44-OME(I,NY)) C C *************CALCULATION OF NUSSELT NUMBER***************** C ANU1=0.0 ANU2=0.0 DO 450 J=1,NY NU(1,J)=(T(1,J)-T(2,J))/DX NU(2,J)=(T(NX1,J)-T(NX,J))/DX 450 CONTINUE DO 460 J=2,NY1 ANU1=ANU1+NU(1,J) ANU2=ANU2+NU(2,J) 460 CONTINUE ANU=(ANU1+ANU2+(NU(1,1)+NU(1,NY)+NU(2,1)+NU(2,NY))*0.5) * /(2*NY1) WRITE(6,470) ANU,RA,ITER 470 FORMAT(2X,'Average Nusselt No.=',F8.5,3X,'Ra=',F9.1, * 3X,'No.of iterations=',I5) C C *******CHECK FOR CONVERGENCE FOR TEMPERATURE FIELD********** C IF(ITER.GE.IAVG) HSUM=HSUM+ANU IF(ITER .LT. IMIN) GOTO 800 IF(ITER .GT. IMAX) GOTO 840 DET=0.0 DO 610 I=2,NX1 DO 620 J=2,NY1 DET=ABS((TOLD(I,J)-T(I,J))/T(I,J)) 620 IF(DET .GT. TOL) GOTO 800 610 CONTINUE GO TO 860 C C ************OUTPUT OF DATA************************************ C 840 WRITE(6,850) WRITE(1,850) 850 FORMAT(' NO CONVERGENCE ') 860 CONTINUE WRITE(6,700) ITER,RA 700 FORMAT(4X,'No. OF ITERATIONS=',I5,5X,'Ra=',F9.1) WRITE(1,705) RA 705 FORMAT(6X,'Ra=',F9.1) WRITE(1,710) 710 FORMAT(4X,'X',6X,'Y',8X,'PSI',12X,'T',12X,'OME') DO 720 I=1,NX XC=DX*(I-1) DO 730 J=1,NY YC=DY*(J-1) WRITE(1,740) XC,YC,PSI(I,J),T(I,J),OME(I,J) 740 FORMAT(1X,F5.2,2X,F5.2,4X,E10.3,4X,F9.4,5X,E10.3) WRITE(2,750) XC,YC,PSI(I,J),T(I,J),OME(I,J) 750 FORMAT(F9.4,',',F9.4,',',F10.4,',',F10.4,',',F10.3) 730 CONTINUE 720 CONTINUE ANUA=HSUM/(ITER-IAVG+1) WRITE(1,470) ANUA,RA,ITER GOTO 999 C C ************FAILURE TO CONVERGE,GO BACK TO LOOP************** C 800 CONTINUE DO 810 I=1,NX DO 820 J=1,NY 820 TOLD(I,J)=T(I,J) 810 CONTINUE IF(ITER .GT. IMAX) GOTO 840 GOTO 199 C 999 CONTINUE CLOSE(1) STOP END ************************************************************************ ***** ***** ***** ENCLPOR ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR NATURAL CONVECTIVE ***** ***** ***** ***** FLOW IN A POROUS MEDIUM FILLED RECTANGULAR ENCLOSURE ***** ***** ***** ***** USING THE FINITE DIFFERENCE TECHNIQUE WITH ***** ***** ***** ***** AN UNDER-RELAXATION ITERATIVE PROCEDURE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ************************************************************************ C DIMENSION PSI(60,60),T(60,60) DIMENSION TOLD(60,60),XNU(2,60) INTEGER CNT C C *************************************************************** C OPEN (1,FILE='ENCLPOR.DAT') C C ******************************************************************** C C **************INPUT OF DATA************************** C NX=33 NY=41 TOL=1.0E-5 IMAX=6000 REX=0.4 C ********************************************************************* C write(6,1400) write(6,1040) read(5,*) RA write(6,1041) read(5,*) AR 1040 format(/////,5X,'Input the Value of the Rayleigh Number ',//) 1041 format(//,5X,'Input the Value of the Aspect Ratio ',//) 1400 format(//////,5X,'FLOW IN A POROUS MEDIA FILLED ENCLOSURE', & /,5X,'-------------------------------------------',///, & /,5X,' Output is written to a file named ENCLPOR.DAT',///) C ********************************************************************* C C ITER=0 CNT=0 DO 100 I=1,NX DO 100 J=1,NY PSI(I,J)=0.0 T(I,J)=1.0-REAL(I-1)/REAL(NX-1) TOLD(I,J)=T(I,J) 100 CONTINUE C C ************BOUNDARY CONDITIONS************************* C DO 120 J=1,NY T(1,J)=1.0 T(NX,J)=0.0 TOLD(1,J)=1.0 TOLD(NX,J)=0.0 120 CONTINUE C C ***********CALCULATION BEGINS************************* C NX1=NX-1 NY1=NY-1 DX=1.0/NX1 DY=AR/NY1 FAC=2*(1/DX**2+1/DY**2) 199 ITER=ITER+1 C C ************CALCULATION OF STREAM FUNCTION "PSI"******** C CNT=0 200 DO 210 I=2,NX1 DO 220 J=2,NY1 F1=1/FAC*(RA*(T(I+1,J)-T(I-1,J))/(2*DX)+ + (PSI(I+1,J)+PSI(I-1,J))/DX**2+(PSI(I,J+1)+PSI(I,J-1))/DY**2) 220 PSI(I,J)=PSI(I,J)+REX*(F1-PSI(I,J)) 210 CONTINUE CNT=CNT+1 IF (CNT .LT. 2) GOTO 200 C C ************CALCULATION OF TEMPERATURE FIELD "T"******** C CNT=0 300 DO 310 I=2,NX1 DO 320 J=2,NY1 F2=1/FAC*((T(I+1,J)+T(I-1,J))/DX**2+ + (T(I,J+1)+T(I,J-1))/DY**2+ + (PSI(I+1,J)-PSI(I-1,J))*(T(I,J+1)-T(I,J-1))/(4*DX*DY)- + (PSI(I,J+1)-PSI(I,J-1))*(T(I+1,J)-T(I-1,J))/(4*DX*DY)) 320 T(I,J)=T(I,J)+REX*(F2-T(I,J)) T(I,1)=T(I,2) T(I,NY)=T(I,NY1) 310 CONTINUE CNT=CNT+1 IF(CNT .LT. 2) GOTO 300 C C *************CALCULATION OF NUSSELT NUMBER***************** C ANU1=0.0 ANU2=0.0 DO 450 J=1,NY XNU(1,J)=(T(1,J)-T(2,J))/DX XNU(2,J)=(T(NX1,J)-T(NX,J))/DX 450 CONTINUE DO 460 J=2,NY1 ANU1=ANU1+XNU(1,J) ANU2=ANU2+XNU(2,J) 460 CONTINUE ANU1=(ANU1+(XNU(1,1)+XNU(1,NY))*0.5)/NY1 ANU2=(ANU2+(XNU(2,1)+XNU(2,NY))*0.5)/NY1 ANU=(ANU1+ANU2)/2.0 C C *******CHECK FOR CONVERGENCE FOR TEMPERATURE FIELD********** C IF(ITER.le.1000) goto 800 DET=0.0 DO 610 I=2,NX1 DO 620 J=2,NY1 DET=ABS((TOLD(I,J)-T(I,J))/T(I,J)) 620 IF(DET .GT. TOL) GOTO 800 610 CONTINUE C C ************OUTPUT OF DATA************************************ C WRITE(1,705) RA 705 FORMAT(6X,'Ra=',F9.1) WRITE(1,710) 710 FORMAT(4X,'X',6X,'Y',8X,'PSI',12X,'T') DO 720 I=1,NX XC=DX*(I-1) DO 730 J=1,NY YC=DY*(J-1) WRITE(1,740) XC,YC,PSI(I,J),T(I,J) 740 FORMAT(1X,F5.2,2X,F5.2,4X,E10.3,4X,F9.4) 730 CONTINUE 720 CONTINUE WRITE(6,470) ANU,RA,ITER 470 FORMAT(2X,'Average Nusselt No.=',F8.5,3X,'Ra=',F9.1, + 3X,'No.of iterations=',I5) GOTO 900 C C ************FAILURE TO CONVERGE,GO BACK TO LOOP************** C 800 CONTINUE DO 810 I=1,NX DO 820 J=1,NY 820 TOLD(I,J)=T(I,J) 810 CONTINUE WRITE(6,830) ITER,ANU1,ANU2,ANU 830 FORMAT(4X,'NO. OF ITERATIONS.=',I5,3X,'NU1=',F9.4, + 3X,'NU2',F9.4,3X,'NUA',F9.4) IF(ITER .GT. IMAX) GOTO 840 GOTO 199 840 WRITE(6,850) WRITE(1,850) 850 FORMAT(' NO CONVERGENCE ') 900 CONTINUE CLOSE(1) STOP END C C ********************************************************************* C C PROGRAM "EXTSQCYL" C C THIS PROGRAM SOLVES THE FLOW OVER A SQUARE CYLINDER. STEADY SYMMETRICAL C FLOW IS ASSUMED. THE FULL NAVIER-STOKES AND ENERGY EQUATIONS C ARE USED IN OBTAINING THE SOLUTION. THE PROGRAM IS INTENDED TO C BE USED TO STUDY LOW REYNOLDS NUMBER FLOWS. C C ********************************************************************* C REAL LU,LD,LN,S(2,91,65),V(2,91,65),T(2,91,65),A(91),B(91),C(91), $ D(91),Q(91),H(91),QL(91) C C ******************************************************************** C OPEN (1,FILE='EXSQCYPR.DAT') OPEN (2,FILE='EXSQCYPL.DAT') C C ********************************************************************* C INPUT C ********************************************************************* C C RE=REYNOLDS NUMBER C PR=PRANDTL NUMBER C LU=DIMENSIONLESS UPSTREAM DISTANCE C LD=DIMENSIONLESS DOWNSTREAM DISTANCE C LN=DIMENSIONLESS CROSS-STREAM DISTANCE C C ********************************************************************* C WRITE(6,3141) 3141 FORMAT(//,' LAMINAR FLOW OVER A SQUARE CYLINDER',/) WRITE(6,3142) 3142 FORMAT(' -----------------------------------',///////) WRITE(6,3143) 3143 FORMAT(' INPUT REYNOLDS AND PRANDTL NUMBERS (RE,PR)',//) READ(5,*) RE,PR WRITE(6,3144) 3144 FORMAT(////,' INPUT DIMENSIONLESS UPSTREAM, DOWNSTREAM',/) WRITE(6,3145) 3145 FORMAT(' AND CROSS-STREAM DISTANCES (LU,LD,LN)',//) READ(5,*) LU,LD,LN WRITE(6,3146) 3146 FORMAT(/////) IF (RE.GT.100.0) THEN WRITE(6,6996) 6996 FORMAT(' ******** RE100 *********') STOP ENDIF C C ********************************************************************* WRITE(1,5237) RE,PR 5237 FORMAT(//,' REYNOLDS NUMBER= ', F12.1,/,' PRANDTL NUMBER= ',F12.3, $//) C ********************************************************************* C REX=0.2 RES=0.7 DX=0.050 NT=15 DY=0.5/(NT-1) C C ********************************************************************* C DX2=DX*DX DY2=DY*DY DX2I=1.0/DX2 DY2I=1.0/DY2 RX2=1.0/(RE*DX2) RY2=1.0/(RE*DY2) PX2=RX2/PR PY2=RY2/PR DXI=1.0/(2.0*DX) DYI=1.0/(2.0*DY) NU=INT(LU/DX)+1 NP=INT(1.0/DX)+NU ND=INT(LD/DX)+NP NN=INT(LN/DY)+1 NN1=NN-1 ND1=ND-1 NU1=NU-1 NP1=NP-1 NT1=NT-1 NUP=NU+1 NTP=NT+1 NPP=NP+1 NTN=NN-NT+1 NDN=ND-NP+1 IW=NP+(ND-NP)/2 ITM=NU+NP/2 C C ********************************************************************* WRITE(2,7237) RE,PR,NU,NP,ND,NT,NN,DX,DY 7237 FORMAT(F12.4,','F12.4,','5(I4,',',)F12.6,',',F12.6) C ********************************************************************* C C ********************************************************************* C INITIAL GUESSED VALUES C ********************************************************************* C DO 100 I=1,ND DO 100 J=1,NN S(1,I,J)=0.0 IF(J.EQ.NN) S(1,I,J)=LN IF(I.EQ.1) S(1,I,J)=(J-1)*LN/(NN-1) IF(I.EQ.ND) S(1,I,J)=(J-1)*LN/(NN-1) IF(I.EQ.NU.AND.J.LE.NT) S(1,I,J)=0.0 IF(I.EQ.NP.AND.J.LE.NT) S(1,I,J)=0.0 V(1,I,J)=0.0 T(1,I,J)=0.0 100 CONTINUE DO 176 I=NU,NP S(1,I,NT)=0.0 176 CONTINUE C C ********************************************************************* C WALL TEMPERATURE VARIATION C ********************************************************************* C DO 200 I=NU,NP T(1,I,NT)=1.0 200 CONTINUE DO 253 J=1,NT T(1,NU,J)=1.0 T(1,NP,J)=1.0 253 CONTINUE C C ********************************************************************* C SOLVE INVISCID FLOW EQUATIONS TO GET INITIAL STREAM FUNCTION VALUES C ********************************************************************* C INTR=0 3100 CONTINUE C C ********************************************************************* C INTR=INTR+1 C C C ********************************************************************* C FIND STREAM FUNCTION ON Y-DIRECTION LINES C ********************************************************************* C DO 5507 J=1,NN S(2,1,J)=S(1,1,J) 5507 CONTINUE DO 5700 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 5600 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 5600 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN CALL TRISOL(NN,A,B,C,D,H) DO 5800 J=1,NN S(2,I,J)=S(1,I,J)+RES*(H(J)-S(1,I,J)) 5800 CONTINUE 5700 CONTINUE DO 6791 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 6691 K=NTP,NN1 J=K-NT+1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,K)+S(1,I-1,K)) 6691 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=LN CALL TRISOL(NTN,A,B,C,D,H) DO 6891 K=1,NTN J=NT+K-1 S(2,I,J)=S(1,I,J)+RES*(H(K)-S(1,I,J)) 6891 CONTINUE 6791 CONTINUE DO 6781 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 9600 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 9600 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=S(1,I,NN) CALL TRISOL(NN,A,B,C,D,H) DO 9800 J=1,NN S(2,I,J)=S(1,I,J)+RES*(H(J)-S(1,I,J)) 9800 CONTINUE 6781 CONTINUE DO 5900 J=1,NN S(2,ND,J)=S(1,ND,J) 5900 CONTINUE C C ********************************************************************* C DO 6207 I=1,ND DO 6207 J=1,NN S(1,I,J)=S(2,I,J) 6207 CONTINUE C C ********************************************************************* C FIND STREAM FUNCTION ON X-DIRECTION LINES C ********************************************************************* C DO 6510 I=1,ND S(2,I,1)=S(1,I,1) 6510 CONTINUE DO 6110 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 7001 I=2,NU1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 7001 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=S(1,NU,J) CALL TRISOL(NU,A,B,C,D,H) DIFFX=0.0 DO 6211 I=1,NU DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+RES*(H(I)-S(1,I,J)) 6211 CONTINUE 6110 CONTINUE DO 5910 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,NP,J) DO 8182 K=NPP,ND1 I=K-NP+1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,K,J+1)+S(1,K,J-1)) 8182 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=0.0 D(NDN)=S(1,ND,J) CALL TRISOL(NDN,A,B,C,D,H) DIFFX=0.0 DO 6011 K=1,NDN I=NP+K-1 DIFF=ABS(H(K)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+RES*(H(K)-S(1,I,J)) 6011 CONTINUE 5910 CONTINUE DO 8197 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 9091 I=2,ND1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 9091 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=0.0 D(ND)=S(1,ND,J) CALL TRISOL(ND,A,B,C,D,H) DO 6097 I=1,ND DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+RES*(H(I)-S(1,I,J)) 6097 CONTINUE 8197 CONTINUE DO 6250 I=1,ND DO 6250 J=1,NN S(1,I,J)=S(2,I,J) 6250 CONTINUE C C ********************************************************************* C BOTTOM OF INVISCID FLOW STREAM FUNCTION LOOP C ********************************************************************* C WRITE(6,5584) INTR 5584 FORMAT(' INVISCID FLOW LOOP, ITERATION NUMBER ',I6) C C ********************************************************************* C IF(INTR.LT.300) GO TO 3100 IF(DIFFX.LT.0.0001) GO TO 4107 IF(INTR.GT.1500) GO TO 4107 GO TO 3100 4107 CONTINUE C C ********************************************************************* C INTR=0 C C ********************************************************************* C TOP OF VORTICITY - STREAM FUNCTION LOOP C ********************************************************************* C 1000 CONTINUE C C ********************************************************************* C INTR=INTR+1 C C ********************************************************************* C CALCULATION OF WALL VORTICITY VALUES C ********************************************************************* C DO 300 I=1,ND IF(I.GE.NU.AND.I.LE.NP) THEN V(2,I,NT)=V(1,I,NT)+REX*(-2.0*(S(1,I,NT+1)- $ S(1,I,NT))/DY2-V(1,I,NT)) ENDIF IF(I.LT.NU) V(2,I,1)=0.0 IF(I.GT.NP) V(2,I,1)=0.0 V(2,I,NN)=0.0 300 CONTINUE DO 353 J=1,NN V(2,1,J)=0.0 IF(J.LE.NT) THEN V(2,NU,J)=V(1,NU,J)+REX*(-2.0*(S(1,NU-1,J)- $ S(1,NU,1))/DY2-V(1,NU,J)) V(2,NP,J)=V(1,NP,J)+REX*(-2.0*(S(1,NP+1,J)- $ S(1,NP,1))/DY2-V(1,NP,J)) ENDIF 353 CONTINUE V(2,NU,NT)=V(1,NU,NT)+REX*(-2.0*(S(1,NU-1,NT)- $ S(1,NU,NT))/DX2-2.0*(S(1,NU,NT+1)- $ S(1,NU,NT))/DY2-V(1,NU,NT)) V(2,NP,NT)=V(1,NP,NT)+REX*(-2.0*(S(1,NP+1,NT)- $ S(1,NP,NT))/DX2-2.0*(S(1,NP,NT+1)- $ S(1,NP,NT))/DY2-V(1,NP,NT)) C C ********************************************************************* C FIND VORTICITY ON Y-DIRECTION LINES C ********************************************************************* C DO 600 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,1) DO 500 J=2,NN1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(V(1,I+1,J)-V(1,I-1,J)) $+RX2*(V(1,I+1,J)+V(1,I-1,J)) 500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 700 J=1,NN V(2,I,J)=V(1,I,J)+REX*(H(J)-V(1,I,J)) 700 CONTINUE 600 CONTINUE DO 691 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,NT) DO 591 K=NTP,NN1 J=K-NT+1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-RY2 C(J)=+DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-RY2 D(J)=-DYI*(S(1,I,K+1)-S(1,I,K-1))*DXI*(V(1,I+1,K)-V(1,I-1,K)) $+RX2*(V(1,I+1,K)+V(1,I-1,K)) 591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=0.0 CALL TRISOL(NTN,A,B,C,D,H) DO 791 K=1,NTN J=NT+K-1 V(2,I,J)=V(1,I,J)+REX*(H(K)-V(1,I,J)) 791 CONTINUE 691 CONTINUE DO 681 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,1) DO 581 J=2,NN1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(V(1,I+1,J)-V(1,I-1,J)) $+RX2*(V(1,I+1,J)+V(1,I-1,J)) 581 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 781 J=1,NN V(2,I,J)=V(1,I,J)+REX*(H(J)-V(1,I,J)) 781 CONTINUE 681 CONTINUE DO 800 J=1,NN V(2,ND,J)=V(2,ND-1,J) 800 CONTINUE DO 2100 I=1,ND DO 2100 J=1,NN V(1,I,J)=V(2,I,J) 2100 CONTINUE C C ********************************************************************* C DO 2410 I=1,ND V(2,I,1)=V(1,I,1) 2410 CONTINUE C C ********************************************************************* C FIND VORTICITY ON X-DIRECTION LINES C ********************************************************************* C DO 1101 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 991 I=2,NU1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(V(1,I,J+1)-V(1,I,J-1)) $+RY2*(V(1,I,J+1)+V(1,I,J-1)) 991 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=V(2,NU,J) CALL TRISOL(NU,A,B,C,D,H) DO 1191 I=1,NU V(2,I,J)=V(1,I,J)+REX*(H(I)-V(1,I,J)) 1191 CONTINUE 1101 CONTINUE DO 1192 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,NP,J) DO 1082 K=NPP,ND1 I=K-NP+1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-RX2 C(I)=-DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-RX2 D(I)=+DXI*(S(1,K+1,J)-S(1,K-1,J))*DYI*(V(1,K,J+1)-V(1,K,J-1)) $+RY2*(V(1,K,J+1)+V(1,K,J-1)) 1082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DO 1374 K=1,NDN I=NP+K-1 V(2,I,J)=V(1,I,J)+REX*(H(K)-V(1,I,J)) 1374 CONTINUE 1192 CONTINUE DO 1010 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 900 I=2,ND1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(V(1,I,J+1)-V(1,I,J-1)) $+RY2*(V(1,I,J+1)+V(1,I,J-1)) 900 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 1100 I=1,ND V(2,I,J)=V(1,I,J)+REX*(H(I)-V(1,I,J)) 1100 CONTINUE 1010 CONTINUE DO 1200 I=1,ND V(2,I,NN)=0.0 1200 CONTINUE DO 2150 I=1,ND DO 2150 J=1,NN V(1,I,J)=V(2,I,J) 2150 CONTINUE C C ********************************************************************* C FIND STREAM FUNCTION ON Y-DIRECTION LINES C ********************************************************************* C DO 3400 J=1,NN S(2,1,J)=S(1,1,J) 3400 CONTINUE DO 3600 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 3500 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,J)-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 3500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN CALL TRISOL(NN,A,B,C,D,H) DO 3700 J=1,NN S(2,I,J)=S(1,I,J)+REX*(H(J)-S(1,I,J)) 3700 CONTINUE 3600 CONTINUE DO 4691 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 4591 K=NTP,NN1 J=K-NT+1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,K)-DX2I*(S(1,I+1,K)+S(1,I-1,K)) 4591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=LN CALL TRISOL(NTN,A,B,C,D,H) DO 4791 K=1,NTN J=NT+K-1 S(2,I,J)=S(1,I,J)+REX*(H(K)-S(1,I,J)) 4791 CONTINUE 4691 CONTINUE DO 4681 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 7500 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,J)-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 7500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN CALL TRISOL(NN,A,B,C,D,H) DO 7700 J=1,NN S(2,I,J)=S(1,I,J)+REX*(H(J)-S(1,I,J)) 7700 CONTINUE 4681 CONTINUE DO 3800 J=1,NN S(2,ND,J)=S(2,ND-1,J) 3800 CONTINUE C C ********************************************************************* C DO 4100 I=1,ND DO 4100 J=1,NN S(1,I,J)=S(2,I,J) 4100 CONTINUE C C ********************************************************************* C FIND STREAM FUNCTION ON X-DIRECTION LINES C ********************************************************************* C DO 4410 I=1,ND S(2,I,1)=0.0 4410 CONTINUE DO 4010 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 4900 I=2,NU1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,I,J)-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 4900 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=S(1,NU,J) CALL TRISOL(NU,A,B,C,D,H) DIFFX=0.0 DO 4111 I=1,NU DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(I)-S(1,I,J)) 4111 CONTINUE 4010 CONTINUE DO 8010 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,NP,J) DO 6082 K=NPP,ND1 I=K-NP+1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,K,J)-DY2I*(S(1,K,J+1)+S(1,K,J-1)) 6082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DO 8111 K=1,NDN I=NP+K-1 DIFF=ABS(H(K)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(K)-S(1,I,J)) 8111 CONTINUE 8010 CONTINUE DO 6091 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 6991 I=2,ND1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,I,J)-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 6991 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 8191 I=1,ND DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(I)-S(1,I,J)) 8191 CONTINUE 6091 CONTINUE DO 4200 I=1,ND S(2,I,NN)=LN 4200 CONTINUE DO 4150 I=1,ND DO 4150 J=1,NN S(1,I,J)=S(2,I,J) 4150 CONTINUE C C ********************************************************************* C BOTTOM OF VORTICITY - STREAM FUNCTION LOOP C ********************************************************************* C WRITE(6,5500) INTR 5500 FORMAT(' VORTICITY-STREAM FUNCTION LOOP, ITERATION NUMBER ',I6) C C ********************************************************************* C IF(INTR.LT.300) GO TO 1000 IF(DIFFX.LT.0.001) GO TO 2000 IF(INTR.GT.1500) GO TO 2000 GO TO 1000 2000 CONTINUE C C ********************************************************************* C DO 6150 I=1,ND DO 6150 J=1,NN T(2,I,J)=T(1,I,J) 6150 CONTINUE INTR=0 C C ********************************************************************* C TOP OF TEMPERATURE LOOP C ********************************************************************* C 6000 CONTINUE C C ********************************************************************* C INTR=INTR+1 C C ********************************************************************* C FIND TEMPERATURE ON Y-DIRECTION LINES C ********************************************************************* C DO 6400 J=1,NN T(2,1,J)=0.0 6400 CONTINUE DO 6600 I=2,NU1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 6500 J=2,NN1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(T(1,I+1,J)-T(1,I-1,J)) $+PX2*(T(1,I+1,J)+T(1,I-1,J)) 6500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 6700 J=1,NN T(2,I,J)=T(1,I,J)+REX*(H(J)-T(1,I,J)) 6700 CONTINUE 6600 CONTINUE DO 6696 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,I,NT) DO 6591 K=NTP,NN1 J=K-NT+1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-PY2 C(J)=+DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-PY2 D(J)=-DYI*(S(1,I,K+1)-S(1,I,K-1))*DXI*(T(1,I+1,K)-T(1,I-1,K)) $+PX2*(T(1,I+1,K)+T(1,I-1,K)) 6591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=0.0 CALL TRISOL(NTN,A,B,C,D,H) DO 6799 K=1,NTN J=NT+K-1 T(2,I,J)=T(1,I,J)+REX*(H(K)-T(1,I,J)) 6799 CONTINUE 6696 CONTINUE DO 6681 I=NPP,ND1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 6581 J=2,NN1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(T(1,I+1,J)-T(1,I-1,J)) $+PX2*(T(1,I+1,J)+T(1,I-1,J)) 6581 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 6783 J=1,NN T(2,I,J)=T(1,I,J)+REX*(H(J)-T(1,I,J)) 6783 CONTINUE 6681 CONTINUE DO 6800 J=1,NN T(2,ND,J)=T(2,ND-1,J) 6800 CONTINUE C C ********************************************************************* C DO 6109 I=1,ND DO 6109 J=1,NN T(1,I,J)=T(2,I,J) 6109 CONTINUE C C ********************************************************************* C FIND TEMPERATURE ON X-DIRECTION LINES C ********************************************************************* C DO 6410 I=1,ND T(2,I,1)=T(1,I,1) 6410 CONTINUE DO 7101 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 7991 I=2,NU1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(T(1,I,J+1)-T(1,I,J-1)) $+PY2*(T(1,I,J+1)+T(1,I,J-1)) 7991 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=T(2,NU,J) CALL TRISOL(NU,A,B,C,D,H) DO 8127 I=1,NU T(2,I,J)=T(1,I,J)+REX*(H(I)-T(1,I,J)) DIFF=ABS(H(I)-T(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF 8127 CONTINUE 7101 CONTINUE DO 7192 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,NP,J) DO 7082 K=NPP,ND1 I=K-NP+1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-PX2 C(I)=-DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-PX2 D(I)=+DXI*(S(1,K+1,J)-S(1,K-1,J))*DYI*(T(1,K,J+1)-T(1,K,J-1)) $+PY2*(T(1,K,J+1)+T(1,K,J-1)) 7082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DIFFX=0.0 DO 7374 K=1,NDN I=NP+K-1 T(2,I,J)=T(1,I,J)+REX*(H(K)-T(1,I,J)) DIFF=ABS(H(K)-T(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF 7374 CONTINUE 7192 CONTINUE DO 7010 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 6900 I=2,ND1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(T(1,I,J+1)-T(1,I,J-1)) $+PY2*(T(1,I,J+1)+T(1,I,J-1)) 6900 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 7100 I=1,ND T(2,I,J)=T(1,I,J)+REX*(H(I)-T(1,I,J)) DIFF=ABS(H(I)-T(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF 7100 CONTINUE 7010 CONTINUE DO 7200 I=1,ND T(2,I,NN)=0.0 T(2,I,1)=T(2,I,2) 7200 CONTINUE DO 6159 I=1,ND DO 6159 J=1,NN T(1,I,J)=T(2,I,J) 6159 CONTINUE C C ********************************************************************* C BOTTOM OF TEMPERATURE LOOP C ********************************************************************* C WRITE(6,5531) INTR 5531 FORMAT(' TEMPERATURE LOOP, ITERATION NUMBER ',I6) C C ********************************************************************* C IF(INTR.LT.300) GO TO 6000 IF(DIFFX.LT.0.001) GO TO 7000 IF(INTR.GT.2500) GO TO 7000 GO TO 6000 7000 CONTINUE C C ********************************************************************* C WRITE OUT STREAM FUNCTION, VORTICITY AND TEMPERATURE VALUES C ********************************************************************* C WRITE(1,5931) 5931 FORMAT(/,' X Y S V $ T',/) DO 6149 I=1,ND DO 6149 J=1,NN X=(I-1)*DX Y=(J-1)*DY WRITE(1,5287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) 5287 FORMAT(5F12.5) IF(I.LE.NU) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF IF(I.GE.NP) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF IF(I.GT.NU.AND.I.LT.NP) THEN IF(J.GE.NT) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF ENDIF 7287 FORMAT(4(F12.5,',',)F12.5) 6149 CONTINUE C C ********************************************************************* C WRITE OUT SURFACE VORTICITY AND HEAT TRANSFER RATE VALUES C ********************************************************************* C WRITE(1,5991) 5991 FORMAT(/,' POINT DISTANCE VORTICITY HEAT FLUX',/) IL=0 DO 5103 J=1,NT IL=IL+1 QL(IL)=(1.0-T(2,NU1,J))/DX XREL=0.0 WRITE(1,5200) IL,XREL,V(2,NU,J),(1.0-T(2,NU1,J))/DX 5103 CONTINUE DO 5100 I=NU,NP IL=IL+1 QL(IL)=(1.0-T(2,I,NTP))/DY XREL=REAL(I-NU)/REAL(NP-NU) WRITE(1,5200) IL,XREL,V(2,I,NT),(1.0-T(2,I,NTP))/DY 5200 FORMAT(I4,3F12.5) 5100 CONTINUE DO 5106 J=1,NT IL=IL+1 K=NT-J+1 QL(IL)=(1.0-T(2,NPP,K))/DX XREL=1.0 WRITE(1,5200) IL,XREL,V(2,NP,K),(1.0-T(2,NPP,K))/DX 5106 CONTINUE C C ********************************************************************* C DETERMINE AND WRITE OUT MEAN HEAT TRANSFER RATE C ********************************************************************* C QM=QL(1)*DY/2.0 DO 4928 J=2,NT1 QM=QM+QL(J)*DY 4928 CONTINUE QM=QM+QL(NT)*DY/2.0 QM=QM+QL(NTP)*DX/2.0 DO 4927 I=NUP,NP1 K=I-NU+1+NT QM=QM+QL(K)*DX 4927 CONTINUE QM=QM+QL(NP-NU+1+NT)*DX/2.0 QM=QM+QL(NP-NU+2+NT)*DY/2.0 DO 4926 J=2,NT1 K=NT+(NP-NU+1)+J QM=QM+QL(K)*DY 4926 CONTINUE QM=QM+QL(2*NT+NP-NU+1)*DY/2.0 WRITE(1,5283) QM/2.0 5283 FORMAT(/,' MEAN NUSSELT NUMBER = ',F12.5) WRITE(6,5247) RE,PR 5247 FORMAT(/,' REYNOLDS NUMBER= ', F12.1,/,' PRANDTL NUMBER= ',F12.3) WRITE(6,5283) QM/2.0 CLOSE(1) CLOSE(2) STOP END C C ********************************************************************* C END OF MAIN PROGRAM C ********************************************************************* C C C ********************************************************************* C SUBROUTINE TRISOL(N,A,B,C,D,H) C C ********************************************************************* C TRI-DIAGONAL MATRIX SOLVER C ********************************************************************* C REAL A(91),B(91),C(91),D(91),H(91),W(91),R(91),G(91) W(1)=A(1) G(1)=D(1)/W(1) DO 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 CONTINUE H(N)=G(N) N1=N-1 DO 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 CONTINUE RETURN END ************************************************************************ ***** ***** ***** LAMBMIX ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** This Program Solves for Mixed (Combined) Convective ***** ***** ***** ***** Laminar Boundary Layer Flow Over a Vertical ***** ***** ***** ***** Flat Plate Using a Finite Difference Technique. ***** ***** ***** ***** A Constant (i.e. Uniform) Cross-stream Grid Spacing is Used.***** ***** ***** ***** There is either a Specified Surface Temperature ***** ***** ***** ***** or a Specified Surface Heat Flux. ***** ***** ***** ************************************************************************ * * * *************************************************** * * ** ** * * ** THIS PROGRAM IS COPYRIGHT BY ** * * ** P. H. OOSTHUIZEN ** * * ** HEAT TRANSFER LABORATORY ( Q'HEAT) ** * * ** DEPARTMENT OF MECHANICAL ENGINEERING ** * * ** QUEEN'S UNIVERSITY ** * * ** KINGSTON, ONTARIO, CANADA ** * * ** TELEPHONE 613-545-2573 ** * * ** FAX 613-545-6489 ** * * ** ** * * *************************************************** * * * ************************************************************************ * dimension U(2,100),V(2,100),T(2,100),A(100),B(100) dimension C(100),D(100),Y(100),H(100) real NX,NL,NNX * *********************** OPEN OUTPUT FILE ***************************** * open(1,file='LAMBMIX.DAT') open(2,file='LAMIXPL.DAT') * ********************************************************************** * write (1,1000) write (1,1001) * ******************** SET INPUT PARAMETERS ************************ * * * GBUYL is the Buoyancy Parameter * * PR is the Prandtl Number * * ISUR Determines the Thermal Boundary Condition at the Wall * * Uniform Wall Temp., ISUR=1, Uniform Wall Heat Flux, ISUR=2 * * IDIR Determines the Direction of the Buoyancy Forces * * Assisting Flow, IDIR=+1, Opposing Flow, IDIR=-1 * * N is Initial Number of Grid Points in Cross-Stream Direction * * DY is the Grid Size in the Cross-Stream Direction * * * ****************************************************************** * write(6,1400) write(6,1050) read(5,*) PR write(6,1100) read(5,*) GBUYL write(6,1500) read(5,*) IDIR write(6,1600) read(5,*) ISUR N = 45 DY = 0.06 * ****************************************************************** * write (1,1002) PR,GBUYL,IDIR,ISUR * ****************************************************************** * X = 0.0 XMAX = 1.00 * Y(1) = 0.0 DY = 0.075 * do 1 J = 2,N Y(J) = Y(J-1) + DY 1 continue * *********************** ASSIGN INITIAL VALUES *********************** * U(1,1) = 0.0 V(1,1) = 0.0 do 2 J = 2,N U(1,J) = 1.0 T(1,J) = 0.0 V(1,J) = 0.0 2 continue * if(ISUR.EQ.1) then T(1,1)=1.0 else T(1,1)=T(1,2)+DY endif * DX = 0.002 DXMAX=0.005 * *************************** SOLUTION BEGINS ************************** * 100 continue VCHX=0.0 * if (DX.LT.DXMAX) then DX = 1.1*DX else DX = DXMAX endif * X = X + DX write (1,2000) X write (6,2001) X * V(2,1) = 0.0 U(2,1) = 0.0 T(2,N) = 0.0 * N1 = N-1 N2 = N-2 * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************* * A(1)=1.0 C(1)=0.0 if(ISUR.EQ.1) then B(1)=0.0 D(1)=1.0 else B(1)=-1.0 D(1)=DY endif do 12 J = 2,N1 A(J) = (2.0/(DY*DY*PR))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) D(J) = U(1,J)*T(1,J)/DX 12 continue * A(N)=1.0 C(N)=0.0 D(N)=0.0 * call TRISOL(N,A,B,C,D,H) * * do 13 J = 1,N T(2,J) = H(J) 13 continue * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ****************** * A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 do 10 J = 2,N1 A(J) = (2.0/(DY*DY))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY) D(J) = U(1,J)*U(1,J)/DX+IDIR*GBUYL*T(2,J) 10 continue * A(N)=1.0 C(N)=0.0 D(N)=1.0 * call TRISOL(N,A,B,C,D,H) * do 11 J = 1,N U(2,J) = H(J) 11 continue * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** * do 14 J = 2, N V(2,J) = V(2,J-1)-(DY/(2.0*DX))*(U(2,J)-U(1,J)+ C U(2,J-1)-U(1,J-1)) 14 continue * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UU = 0.99 TT = 0.01*T(2,1) do 20 J=10,N1 if (U(2,J).GT.UU) go to 21 20 continue 21 DU =Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) do 22 J=1,N1 if (T(2,J).le.TT) go to 23 22 continue 23 DT = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * if (DU.LE.Y(N-3).and.DT.LE.Y(N-3)) go to 24 NMIN = N+1 NMAX = N+5 if (NMAX.GT.100) go to 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) do 25 I=NMIN,NMAX Y(I) = Y(I-1) + DY U(2,I) = 1.0 T(2,I) = 0.0 V(2,I)=V(2,I-1)+DVY*DY 25 continue write (1,3000) write (6,3000) N = NMAX 24 continue * * ************************* TRANSFER VALUES *************************** * do 30 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 30 continue * ************************ WRITE OUT THE RESULTS ************************ * DTY = (T(2,1)-T(2,2))/(Y(2)-Y(1)) NL = DTY/T(2,1) NX = DTY*(X**0.50)/T(2,1) GBUYX=GBUYL*X NNX=NX/GBUYX**0.25 write(2,4999)X,GBUYX,DU,DT,NL,NX,NNX write(1,4000)DTY,NL,NX,DU,DT,GBUYX write(1,4001) DO 40 J=1,N write(1,4002)J,Y(J),U(2,J),V(2,J),T(2,J) 40 continue * ************************* CHECK FOR MAXIMUM X ************************ * if(X.ge.XMAX) go to 102 if(U(2,2).le.0.0) go to 101 go to 100 101 write (1,5000) write (6,5000) go to 103 102 write(1,5001) write(6,5001) 103 continue * stop close(1) * ************************** FORMAT STATEMENTS ************************** * 1000 format (10X,'MIXED CONVECTIVE LAMINAR BOUNDARY FLOW') 1001 format (10X,'--------------------------------------',//) 1002 format (5X,'PRANDTL NUMBER=',F10.4,/, + 5X,'BOUYANCY FORCE PARAMETER =',E12.4,/, + 5X,'BUOYANCY FORCE DIRECTION PARAMETER =',I4,/, + 5X,'WALL THERMAL CONDITION PARAMETER =',I4,///) 1050 format(/////,5X,'Input the Value of the Prandtl Number ',/) 1100 format(//,5X,'Input the Value of the Buoyancy Parameter ',/) 1400 format(////////,5X,'MIXED CONVECTION OVER A VERTICAL FLAT PLATE', & /,5X,'-------------------------------------------',///, & /,5X,' Output is written to a file named LAMBMIX.DAT',///) 1500 format(//,5X,'Assisting Flow (+1) or Opposing Flow (-1) ',/) 1600 format(//,5X,'Uniform Temp. (1) or Uniform Heat Flux (2) ',/) 2000 format (/,5X,'SOLUTION FOR X =',F12.4,/ + ,5X,'___________________________',/) 2001 format (5X,'SOLUTION FOR X =',F12.4,' COMPLETED') 3000 format (5X,'****** NUMBER OF GRID POINTS INCREASED ******') 4000 format (5X,'DT/DY Y=0 =',E12.4,/,5X,'NL =',E12.4,/, + 5X,'NX/RX^1/2 =',E12.4,/,5X, + 'DELTU=',E12.4,/,5X,'DELTT=',E12.4,/,5X,'GX/RX^2=',F12.4,/) 4001 format (2X,'J',7X,'Y',12X,'U',12X,'V',12X,'T') 4002 format (I4,4E12.4) 4999 format (7F11.5) 5000 format (/,5X,'********* BUOYANCY INDUCED SEPARATION **********') 5001 format (/,5X,'********** XMAX REACHED **********') END * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * ******************************************************************** * dimension A(100),B(100),C(100),D(100),H(100),W(100),Q(100),G(100) W(1)=A(1) G(1)=D(1)/W(1) do 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 continue H(NN)=G(NN) N1=NN-1 do 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 continue return end ************************************************************************ ***** ***** ***** LAMBNAT ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR NATURAL CONVECTIVE ***** ***** ***** ***** LAMINAR BOUNDARY LAYER FLOW OVER A VERTICAL ***** ***** ***** ***** SURFACE USING THE FINITE DIFFERENCE TECHNIQUE WITH ***** ***** ***** ***** AN UNDER-RELAXATION ITERATIVE PROCEDURE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ************************************************************************ * DIMENSION U(2,100),V(2,100),T(2,100),A(100),B(100) DIMENSION C(100),D(100),Y(100),H(100) REAL NX * ********************************************************************** * OPEN(1,FILE='C:LAMBNAT.DAT') * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) * ************************** ASSIGN KNOWN VALUES *********************** * WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) XMAX,PR X = 0.0 N = 35 Y(1) = 0.0 DY = 0.075 REX=0.2 * DO 1 J = 2,N Y(J) = Y(J-1) + DY 1 CONTINUE * *********************** ASSIGN INITIAL VALUES *********************** * U(1,1) = 0.0 * V(1,1) = 0.0 * CALL TEMP (X,T(1,1)) * DO 2 J = 2,N U(1,J) = 0.0 T(1,J) = 0.0 V(1,J) = 0.0 2 CONTINUE DX = 0.005 DXMAX=0.05 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE ITER=0 VCHX=0.0 DO 120 J = 1,N U(2,J) = U(1,J) T(2,J) = T(1,J) V(2,J) = V(1,J) 120 CONTINUE * IF (DX.LT.DXMAX) THEN DX = 1.1*DX ELSE DX = DXMAX ENDIF * X = X + DX WRITE (1,2000) X WRITE (6,2000) X * V(2,1) = 0.0 U(2,1) = 0.0 T(2,N) = 0.0 * CALL TEMP (X,T(2,1)) * 500 ITER=ITER+1 * N1 = N-1 N2 = N-2 * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,1) DO 12 J = 2,N1 A(J) = (2.0/(DY*DY*PR))+(U(2,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) D(J) = U(2,J)*T(1,J)/DX 12 CONTINUE * A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=0.0 * CALL TRISOL(N,A,B,C,D,H) * * DO 13 J = 1,N T(2,J) = T(2,J)+REX*(H(J)-T(2,J)) 13 CONTINUE * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ******************* * A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 10 J = 2,N1 A(J) = (2.0/(DY*DY))+(U(2,J)/DX) B(J) = V(2,J)/(2.0*DY)-1.0/(DY*DY) C(J) = -V(2,J)/(2.0*DY)-1.0/(DY*DY) D(J) = U(2,J)*U(1,J)/DX+T(2,J) 10 CONTINUE * A(N)=1.0 C(N)=0.0 D(N)=0.0 * CALL TRISOL(N,A,B,C,D,H) * DO 11 J = 1,N U(2,J) = U(2,J)+REX*(H(J)-U(2,J)) 11 CONTINUE * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** * DO 14 J = 2, N V(2,J) = V(2,J-1)-(DY/(2.0*DX))*(U(2,J)-U(1,J)+ C U(2,J-1)-U(1,J-1)) 14 CONTINUE IF(ITER.LT.50) GO TO 500 IF(ITER.GT.75) GO TO 101 VDIFF=ABS(V(2,N)-VCHX) IF(VDIFF.LT.0.01) GO TO 300 VCHX=V(2,N) GO TO 500 300 CONTINUE * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UX=0.0 DO 26 J=1,N1 IF (U(2,J).GT.UX) UX=U(2,J) 26 CONTINUE UU = 0.01*UX TT = 0.01*T(2,1) DO 20 J=5,N1 IF (U(2,J).LT.UU) GO TO 21 20 CONTINUE 21 DU =Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) DO 22 J=1,N1 IF (T(2,J).LE.TT) GO TO 23 22 CONTINUE 23 DT = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * IF (DU.LE.Y(N-3).AND.DT.LE.Y(N-3)) GO TO 24 NMIN = N+1 NMAX = N+5 IF (NMAX.GT.100) GO TO 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) DO 25 I=NMIN,NMAX Y(I) = Y(I-1) + DY U(2,I) = 0.0 T(2,I) = 0.0 V(2,I)=V(2,I-1)+DVY*DY 25 CONTINUE WRITE (1,3000) WRITE (6,3000) N = NMAX 24 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 30 CONTINUE * ************************** WRITE THE RESULTS ************************** * DTY = (T(2,1)-T(2,2))/(Y(2)-Y(1)) NX = DTY*(X**0.25)/T(2,1) WRITE(1,4000)DTY,NX,DU,DT WRITE(6,4500)ITER,DTY,NX,DU,DT WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,Y(J),U(2,J),V(2,J),T(2,J) 40 CONTINUE * ************************** CHECK FOR MAXIMUM X ************************ * IF(X.GE.XMAX) GO TO 102 GO TO 100 101 WRITE (1,5000) WRITE (6,5000) GO TO 103 102 WRITE(1,5001) WRITE(6,5001) 103 CONTINUE STOP CLOSE(1) * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'LAMINAR BOUNDARY FLOW') 1001 FORMAT (10X,'---------------------',//) 2000 FORMAT (/,5X,'SOLUTION FOR X =',F12.4,/ C ,5X,'___________________________',/) 3000 FORMAT (5X,'****** NUMBER OF GRID POINTS INCREASED ******') 4000 FORMAT (5X,'DT/DY Y=0 =',E12.4,/,5X,'NX/GX 1/4 =',E12.4,/,5X, + 'DELTU=',E12.4,/,5X,'DELTT=',E12.4,/) 4500 FORMAT (5X,'ITER =',I4,/,5X, + 'DT/DY Y=0 =',E12.4,/,5X,'NX/GX 1/4 =',E12.4,/,5X, + 'DELTU=',E12.4,/,5X,'DELTT=',E12.4,/) 4001 FORMAT (2X,'J',7X,'Y',12X,'U',12X,'V',12X,'T') 4002 FORMAT (I4,4E12.4) 5000 FORMAT (5X,'ITERATION NUMBER TOO LARGE OR N EXCEEDS 100') LAM02090 5001 FORMAT (/,5X,'********** XMAX REACHED **********') 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUES OF THE MAXIMUM X TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN AND', & /,' OF THE PRANDTL NUMBER - (XMAX,PR)',//) END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(100),B(100),C(100),D(100),H(100),W(100),Q(100),G(100) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END * * ******************************************************************** * SUBROUTINE TEMP(X,TW) * *********** THIS DETERMINES THE WALL TEMPERATURE ******************* * ---- * TW = 1.0 C TW = X RETURN END ************************************************************************ ***** ***** ***** LAMBOUN ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR LAMINAR BOUNDARY LAYER FLOW OVER ***** ***** --- ---- ***** ***** A SURFACE USING THE FINITE DIFFERENCE TECHNIQUE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ***** THE DIMENSIONLESS WALL TEMPERATURE AND FREESTREAM VELOCITY ***** ***** ***** ***** DISTRIBUTIONS ARE DESCRIBED BY 3RD-ORDER POLYNOMIALS. ***** ***** ***** ************************************************************************ C DIMENSION U(2,200),V(2,200),T(2,200),A(200),B(200) DIMENSION C(200),D(200),Y(200),H(200) REAL NX COMMON AT,BT,CT,DT,AV,BV,CV,DV * ********************************************************************** * OPEN(1,FILE='LAMBPRT.DAT') OPEN(2,FILE='LAMBPLT.DAT') * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) * ************************** ASSIGN KNOWN VALUES *********************** * WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) XMAX,PR WRITE(6,6002) READ(5,*) AT,BT,CT,DT WRITE(6,6003) READ(5,*) AV,BV,CV,DV WRITE (1,1002) XMAX,PR,AT,BT,CT,DT,AV,BV,CV,DV * X = 0.0 N = 49 Y(1) = 0.0 DY = 0.050 * DO 1 J = 2,N Y(J) = Y(J-1) + DY 1 CONTINUE * *********************** ASSIGN INITIAL VALUES *********************** * U(1,1) = 0.0 * CALL TEMP(X,TW) T(1,1)=TW * V(1,1) = 0.0 * DO 2 J = 2,N U(1,J) = 1.0 T(1,J) = 0.0 V(1,J) = 0.0 2 CONTINUE DX = 0.0015 DXMAX=0.01 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * IF (DX.LT.DXMAX) THEN DX = 1.05*DX ELSE DX = DXMAX ENDIF * X = X + DX WRITE (1,2000) X * V(2,1) = 0.0 U(2,1) = 0.0 T(2,N) = 0.0 * CALL TEMP (X,TW) CALL VELDP (X,UF,DPX) * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ******************* * N1 = N-1 N2 = N-2 DO 10 J = 2,N1 A(J) = (2.0/(DY*DY))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY) D(J) = U(1,J)*U(1,J)/DX-DPX 10 CONTINUE * A(1)=1.0 B(1)=0.0 D(1)=0.0 A(N)=1.0 C(N)=0.0 D(N)=UF * CALL TRISOL(N,A,B,C,D,H) * DO 11 J = 1,N U(2,J) = H(J) 11 CONTINUE * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * DO 12 J = 1,N A(J) = (2.0/(DY*DY*PR))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) D(J) = U(1,J)*T(1,J)/DX 12 CONTINUE * A(1)=1.0 D(1)=TW B(1)=0.0 A(N)=1.0 C(N)=0.0 D(N)=0.0 * CALL TRISOL(N,A,B,C,D,H) * * DO 13 J = 1,N T(2,J) = H(J) 13 CONTINUE * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** LAM01330 * DO 14 J = 2, N V(2,J) = V(2,J-1)-(DY/(2.0*DX))*(U(2,J)-U(1,J)+ C U(2,J-1)-U(1,J-1)) 14 CONTINUE * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UU = 0.99*U(2,N) TT = 0.01*T(2,1) DO 20 J=1,N1 IF (U(2,J).GE.UU) GO TO 21 20 CONTINUE 21 DU =Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) DO 22 J=1,N1 IF (T(2,J).LE.TT) GO TO 23 22 CONTINUE 23 DE = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * IF (DU.LE.Y(N-5).AND.DE.LE.Y(N-5)) GO TO 24 NMIN = N+1 NMAX = N+10 IF (NMAX.GT.200) GO TO 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) DO 25 I=NMIN,NMAX Y(I) = Y(I-1) + DY U(2,I) = U(2,N) T(2,I) = T(2,N) V(2,I)=V(2,I-1)+DVY*DY 25 CONTINUE WRITE (1,3000) WRITE (6,3000) N = NMAX 24 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 30 CONTINUE * ************************** WRITE THE RESULTS ************************** * DTY = (T(2,1)-T(2,2))/(Y(2)-Y(1)) NX = DTY*SQRT(X)/T(2,1) WRITE(1,4000)DTY,NX,T(2,1),DU,DE WRITE(6,4200) X WRITE(2,4100) X,DTY,NX,T(2,1),DU,DE WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,Y(J),U(2,J),V(2,J),T(2,J) 40 CONTINUE * ************************** CHECK FOR MAXIMUM X ************************ * IF(X.GE.XMAX) GO TO 102 GO TO 100 101 WRITE (1,5000) WRITE (6,5000) GO TO 103 102 WRITE(1,5001) WRITE(6,5001) 103 CONTINUE STOP CLOSE(1) * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'LAMINAR BOUNDARY FLOW') 1001 FORMAT (10X,'---------------------',//) 1002 FORMAT (5X,' MAXIMUM X =', F10.4,/,5X,' PRANDTL NO. =',F10.4, $ /,5X, 'TEMPERATURE COEFFS. =',4F10.4, $ /,5X, 'VELOCITY COEFFS. =',4F10.4,//) 2000 FORMAT (/,5X,'SOLUTION FOR X =',F12.4,/ & ,5X,'___________________________',/) 3000 FORMAT (5X,'****** NUMBER OF GRID POINTS INCREASED ******') 4000 FORMAT (5X,'DT/DY Y=0 =',E12.4,/,5X,'NX/RX 1/2 =',E12.4,/,5X, & 'TWALL=',E12.4,/,5X,'DELTU=',E12.4,/,5X,'DELTT=',E12.4,/) 4001 FORMAT (2X,'J',7X,'Y',12X,'U',12X,'V',12X,'T') 4002 FORMAT (I4,4E12.4) 4100 FORMAT (F10.5,',',F10.5,',',F10.5,',',F10.5,',',F10.5,',',F10.5) 4200 FORMAT (5X,' X =',F10.4) 5000 FORMAT (5X,'N EXCEEDS 200') 5001 FORMAT (/,5X,'********** XMAX REACHED **********') 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUES OF THE MAXIMUM X TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN AND', & /,' OF THE PRANDTL NUMBER - (XMAX,PR)',//) 6002 FORMAT(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE WALL TEMPERATURE VARIATION - & (AT,BT,CT,DT) ',//) 6003 FORMAT(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE FREESTEAM VELOCITY VARIATION - & (AV,BV,CV,DV) ',//) END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(200),B(200),C(200),D(200),H(200),W(200),Q(200),G(200) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END * * ******************************************************************** * SUBROUTINE TEMP(X,TW) COMMON AT,BT,CT,DT,AV,BV,CV,DV * *********** THIS DETERMINES THE WALL TEMPERATURE ******************* * ---- * TW = AT+BT*X+CT*X*X+DT*X*X*X RETURN END * * ******************************************************************** * SUBROUTINE VELDP (X,UF,DPX) COMMON AT,BT,CT,DT,AV,BV,CV,DV * ****** THIS DETERMINES FREESTREAM VELOCITY AND PRESSURE GRADIENT ****** * --- - * UF = AV+BV*X+CV*X*X+DV*X*X*X DPX = -(BV+2.0*CV*X+3.0*DV*X*X) RETURN END ************************************************************************ ***** ***** ***** LAMBOUQ ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR LAMINAR BOUNDARY LAYER FLOW OVER ***** ***** ***** ***** A SURFACE USING THE FINITE DIFFERENCE TECHNIQUE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) GRID SPACING IS USED. ***** ***** ***** ***** THE DIMENSIONLESS WALL HEAT FLUX AND FREESTREAM VELOCITY ***** ***** ***** ***** DISTRIBUTIONS ARE DESCRIBED BY 3RD-ORDER POLYNOMIALS. ***** ***** ***** ************************************************************************ C DIMENSION U(2,200),V(2,200),T(2,200),A(200),B(200) DIMENSION C(200),D(200),Y(200),H(200) REAL NX COMMON AQ,BQ,CQ,DQ,AV,BV,CV,DV * ********************************************************************** * OPEN(1,FILE='LAMQPRT.DAT') OPEN(2,FILE='LAMQPLT.DAT') * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) * ************************** ASSIGN KNOWN VALUES *********************** * WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6001) READ(5,*) XMAX,PR WRITE(6,6002) READ(5,*) AQ,BQ,CQ,DQ WRITE(6,6003) READ(5,*) AV,BV,CV,DV WRITE (1,1002) XMAX,PR,AQ,BQ,CQ,DQ,AV,BV,CV,DV * X = 0.0 N = 49 Y(1) = 0.0 DY = 0.050 * DO 1 J = 2,N Y(J) = Y(J-1) + DY 1 CONTINUE * *********************** ASSIGN INITIAL VALUES *********************** * U(1,1) = 0.0 * CALL FLUX(X,QW) T(1,1)=QW*DY * V(1,1) = 0.0 * DO 2 J = 2,N U(1,J) = 1.0 T(1,J) = 0.0 V(1,J) = 0.0 2 CONTINUE DX = 0.0015 DXMAX=0.01 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * IF (DX.LT.DXMAX) THEN DX = 1.05*DX ELSE DX = DXMAX ENDIF * X = X + DX WRITE (1,2000) X * V(2,1) = 0.0 U(2,1) = 0.0 T(2,N) = 0.0 * CALL FLUX (X,TW) CALL VELDP (X,UF,DPX) * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ******************* * N1 = N-1 N2 = N-2 DO 10 J = 2,N1 A(J) = (2.0/(DY*DY))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY) D(J) = U(1,J)*U(1,J)/DX-DPX 10 CONTINUE * A(1)=1.0 B(1)=0.0 D(1)=0.0 A(N)=1.0 C(N)=0.0 D(N)=UF * CALL TRISOL(N,A,B,C,D,H) * DO 11 J = 1,N U(2,J) = H(J) 11 CONTINUE * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * DO 12 J = 1,N A(J) = (2.0/(DY*DY*PR))+(U(1,J)/DX) B(J) = V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) C(J) = -V(1,J)/(2.0*DY)-1.0/(DY*DY*PR) D(J) = U(1,J)*T(1,J)/DX 12 CONTINUE * A(1)=1.0 B(1)=-1.0 D(1)=DY*QW A(N)=1.0 C(N)=0.0 D(N)=0.0 * CALL TRISOL(N,A,B,C,D,H) * * DO 13 J = 1,N T(2,J) = H(J) 13 CONTINUE * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** LAM01330 * DO 14 J = 2, N V(2,J) = V(2,J-1)-(DY/(2.0*DX))*(U(2,J)-U(1,J)+ C U(2,J-1)-U(1,J-1)) 14 CONTINUE * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UU = 0.99*U(2,N) TT = 0.01*T(2,1) DO 20 J=1,N1 IF (U(2,J).GE.UU) GO TO 21 20 CONTINUE 21 DU =Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) DO 22 J=1,N1 IF (T(2,J).LE.TT) GO TO 23 22 CONTINUE 23 DE = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * IF (DU.LE.Y(N-5).AND.DE.LE.Y(N-5)) GO TO 24 NMIN = N+1 NMAX = N+10 IF (NMAX.GT.200) GO TO 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) DO 25 I=NMIN,NMAX Y(I) = Y(I-1) + DY U(2,I) = U(2,N) T(2,I) = T(2,N) V(2,I)=V(2,I-1)+DVY*DY 25 CONTINUE WRITE (1,3000) WRITE (6,3000) N = NMAX 24 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 30 CONTINUE * ************************** WRITE THE RESULTS ************************** * QX = QW/T(2,1) NX = QW*SQRT(X)/T(2,1) WRITE(1,4000) QX,NX,T(2,1),DU,DE WRITE(6,4200) X,T(2,1) WRITE(2,4100) X,QX,NX,T(2,1),DU,DE WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,Y(J),U(2,J),V(2,J),T(2,J) 40 CONTINUE * ************************** CHECK FOR MAXIMUM X ************************ * IF(X.GE.XMAX) GO TO 102 GO TO 100 101 WRITE (1,5000) WRITE (6,5000) GO TO 103 102 WRITE(1,5001) WRITE(6,5001) 103 CONTINUE STOP CLOSE(1) * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'LAMINAR BOUNDARY FLOW') 1001 FORMAT (10X,'---------------------',//) 1002 FORMAT (5X,' MAXIMUM X =', F10.4,/,5X,' PRANDTL NO. =',F10.4, $ /,5X, 'WALL HEAT FLUX COEFFS. =',4F10.4, $ /,5X, 'VELOCITY COEFFS. =',4F10.4,//) 2000 FORMAT (/,5X,'SOLUTION FOR X =',F12.4,/ & ,5X,'___________________________',/) 3000 FORMAT (5X,'****** NUMBER OF GRID POINTS INCREASED ******') 4000 FORMAT (5X,'QW/TW =',E12.4,/,5X,'NX/RX 1/2 =',E12.4,/,5X, & 'TWALL=',E12.4,/,5X,'DELTU=',E12.4,/,5X,'DELTT=',E12.4,/) 4001 FORMAT (2X,'J',7X,'Y',12X,'U',12X,'V',12X,'T') 4002 FORMAT (I4,4E12.4) 4100 FORMAT (F10.5,',',F10.5,',',F10.5,',',F10.5,',',F10.5,',',F10.5) 4200 FORMAT (5X,' X =',F10.4,2X,' TW =',F10.4) 5000 FORMAT (5X,'N EXCEEDS 200') 5001 FORMAT (/,5X,'********** XMAX REACHED **********') 6000 FORMAT(////////) 6001 FORMAT(/,' INPUT THE VALUES OF THE MAXIMUM X TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN AND THE PRANDTL - (XMAX,PR)',//) 6002 FORMAT(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE WALL HEAT FLUX VARIATION - & (AQ,BQ,CQ,DQ) ',//) 6003 FORMAT(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE FREESTEAM VELOCITY VARIATION - & (AV,BV,CV,DV) ',//) END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(200),B(200),C(200),D(200),H(200),W(200),Q(200),G(200) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END * * ******************************************************************** * SUBROUTINE FLUX(X,QW) COMMON AQ,BQ,CQ,DQ,AV,BV,CV,DV * *********** THIS DETERMINES THE WALL HEAT FLUX ********************* * ---- * QW = AQ+BQ*X+CQ*X*X+DQ*X*X*X RETURN END * * ******************************************************************** * SUBROUTINE VELDP (X,UF,DPX) COMMON AQ,BQ,CQ,DQ,AV,BV,CV,DV * ****** THIS DETERMINES FREESTREAM VELOCITY AND PRESSURE GRADIENT ****** * --- - * UF = AV+BV*X+CV*X*X+DV*X*X*X DPX = -(BV+2.0*CV*X+3.0*DV*X*X) RETURN END ************************************************************************ ***** ***** ***** NATCHAN ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR NATURAL CONVECTIVE ***** ***** ***** ***** STEADY LAMINAR FLOW THROUGH A VERTICAL PLANE CHANNEL. ***** ***** ***** ***** THE FLOW IS ASSUMED TO BE SYMMETRICAL ABOUT THE ***** ***** ***** ***** CENTRE-LINE. A UNIFORM WALL TEMPERATURE IS ASSUMED. ***** ***** ***** ***** THE PARABOLIC FORM OF THE GOVERNING EQUATIONS ARE SOLVED ***** ***** ***** ***** USING A DIFFERENCE TECHNIQUE WITH A UNIFORM GRID ***** ***** ***** ***** SPACING IN THE CROSS-STREAM DIRECTION. ***** ***** ***** ************************************************************************ * PARAMETER (NM=41) DIMENSION U(NM), V(NM), T(NM),U0(NM),V0(NM),T0(NM) DIMENSION C0(NM),D0(NM),E0(NM) * ********************************************************************** * OPEN(1,FILE='C:NATCHAN.DAT') OPEN(2,FILE='C:NATCHAN.VAR') * ********************************************************************** * C ***************PRINT HEADING********************************** C WRITE (1,10) 10 FORMAT (4X,'Z',12X,'Y',11X,'U',8X,'V',8X,'T',10X,'P', $ 10X,'Q') WRITE (2,11) 11 FORMAT (5X,'Z',11X,'Q',10X,'P',11X,'UC',10X,'TC') C C **************SPECIFY INPUT DATA************************ C WRITE(6,1007) 1007 FORMAT(//,' NATURAL CONVECTION IN A VERTICAL PLANE CHANNEL', & /,'------------------------------------------------',///) WRITE(6,1005) 1005 FORMAT(' INPUT PRANDTL NUMBER') READ(5,*)PR WRITE(6,1006) 1006 FORMAT(' INPUT DIMENSIONLESS INLET VELOCITY') READ(5,*)UIN C PR=0.7 C UIN=0.005 C C ******************************************************************** C NY=NM NY1=NY-1 DY=0.5/NY1 DY2=2.0*DY DYSQ=DY**2 IPF=1E+3 IZM=1E+6 K=1 C C ***********INATIALIZE VALUES****************************** C P0=-(UIN**2)/2.0 DO 20 J=1,NY U0(J)=UIN T0(J)=0.0 V0(J)=0.0 20 CONTINUE Z=0.0 IZ=1 IP=1 IP2=0 C C ***********APPLY BOUNDARY CONDITIONS********************** C U0(NY)=0.0 T0(NY)=1.0 V0(NY)=0.0 U(NY)=0.0 T(NY)=1.0 V(NY)=0.0 V(1)=0.0 C C ***********SOLUTION BEGINS********************************* C 30 CONTINUE WRITE (6,1350) Z,Q,P 1350 FORMAT (2X,'Z=',F10.8,3X,'Q=',F10.8,3X,'P=',F10.8) C C ***********FIND Z-STEP************************************* C UMIN=U0(1) IF (U0(NY1) .LT. UMIN) UMIN=U0(NY1) DZ=DY*DY*UMIN/8.0 Z=Z+DZ C C ***********CALCULATE C,D,E AND P*************************** C DO 50 J=2,NY1 FC1=U0(J)*U0(J)/DZ FC2=V0(J)*(U0(J+1)-U0(J-1))/DY2 FC3=P0/DZ FC4=(U0(J+1)+U0(J-1)-2.0*U0(J))/DYSQ FC5=T0(J) C0(J)=(FC1-FC2+FC3+FC4+FC5)*DZ/U0(J) 50 CONTINUE C0(1)=(1.0/U0(1))*(U0(1)**2+P0+2.0*(U0(2)-U0(1)) $ *DZ/DY**2+T0(1)*DZ) DO 55 J=2,NY1 D0(J)=-(C0(J)+C0(J-1)-U0(J)-U0(J-1))*DY/(2.0*DZ) E0(J)=-(1.0/U0(J)+1.0/U0(J-1))*DY/(2.0*DZ) 55 CONTINUE D0(NY)=-(C0(NY1)-U0(NY1))*DY/(2.0*DZ) E0(NY)=-DY/(2.0*DZ*U0(NY1)) DSUM=0.0 ESUM=0.0 DO 60 J=2,NY DSUM=DSUM+D0(J) ESUM=ESUM+E0(J) 60 CONTINUE P=DSUM/ESUM C C ************CALCULATE U,V,T**************************** C DO 70 J=2, NY1 U(J)=C0(J)-P/U0(J) V(J)=V(J-1)+D0(J)-E0(J)*P FT1=T0(J) FT2=V0(J)*(T0(J+1)-T0(J-1))/DY2*DZ/U0(J) FT3=(T0(J+1)+T0(J-1)-2.0*T0(J))/PR/DYSQ*DZ/U0(J) T(J)=FT1-FT2+FT3 70 CONTINUE C C **********APPLY SYMMETRY AXIS CONDITION***************** C U(1)=4.0*U(2)/3.0-U(3)/3.0 T(1)=4.0*T(2)/3.0-T(3)/3.0 C C *************CALCULATE HEAT TRANSFER RATE************** C Q=0.0 DO 250 J=1,NY1 DIFQ=(U(J)*T(J)+U(J+1)*T(J+1))*DY/2.0 Q=Q+DIFQ 250 CONTINUE C C *************OUTPUT DATA******************************** C IP2=IP2+1 IF (IP2.EQ.5) THEN IP2=0 WRITE (2,1360) Z,Q,P,U(1),T(1) 1360 FORMAT (F11.8,',',F11.8,',',F11.8,',',F11.8,',',F11.8) ENDIF IF (IP .GT. IPF) THEN IP=1 Y=0.0 DO 350 J=1,NY,2 Y=DY*(J-1) WRITE (1,320) Z,Y,U(J),V(J),T(J),P,Q 320 FORMAT (1X,E12.5,1X,E10.3,3(1X,F8.4),2(1X,E10.3)) 350 CONTINUE WRITE (1,380) 380 FORMAT (/) ELSE IP=IP+1 ENDIF C C *************CHECK FOR P******************************* C IZ=IZ+1 IF (IZ .GT. IZM) GOTO 951 IF (P .GT. 0.0) GOTO 900 DO 500 J=1,NY U0(J)=U(J) V0(J)=V(J) T0(J)=T(J) 500 CONTINUE P0=P GOTO 30 C C *************PRINT CHANNEL LENGTH********************** C 900 CONTINUE Z=Z-DZ-P0*DZ/(P-P0) Q=Q-DIFQ-P0*DIFQ/(P-P0) QA=Q/Z WRITE (1,910) Z,UIN,Q,QA STOP 910 FORMAT (2X,'Z=',F10.8,3X,'UIN=',F10.8,3X,'Q=',F10.8, & 3X,'QMEAN=',F10.5) 951 CONTINUE WRITE (1,961) STOP 961 FORMAT (2X,'**** MAX. ALLOWABLE NUMBER OF STEPS EXCEEDED ****') END C C ********************************************************************* C C PROGRAM "PIPETEM" C C THIS PROGRAM SOLVES FOR THE HEAT TRANSFER RATE IN FULLY- C DEVELOPED LAMINAR FLOW IN A PIPE WITH A UNIFORM WALL TEMPERATURE. C C ********************************************************************* C REAL G(1500),F(1500),H(1500),R(1500) C C ******************************************************************** C INPUT C ********************************************************************* C WRITE(6,6000) 6000 FORMAT(////////) WRITE(6,6001) 6001 FORMAT (10X,'FULLY-DEV. FLOW IN A PIPE WITH A UNIF. WALL TEMP.',/) WRITE(6,6002) 6002 FORMAT (10X,'-------------------------------------------------',/) NR=1500 DR=1.0/(NR-1) C C ********************************************************************* C NR1=NR-1 C C ********************************************************************* C INITIAL GUESSED VALUES C ********************************************************************* C DO 100 I=1,NR R(I)=(I-1)*DR G(I)=1.0-4.0*R(I)*R(I)/3.0+R(I)*R(I)*R(I)*R(I)/3.0 100 CONTINUE INTR=0 XNUP=0.0 C C ********************************************************************* C BEGIN SOLUTION C ********************************************************************* C 1000 CONTINUE INTR=INTR+1 C C ********************************************************************* C DO 200 I=1,NR H(I)=G(I)*(R(I)-R(I)*R(I)*R(I)) 200 CONTINUE F(1)=0.0 DO 300 I=2,NR F(I)=(F(I-1)+(H(I)+H(I-1))*DR/2.0) 300 CONTINUE DO 350 I=2,NR F(I)=F(I)/R(I) 350 CONTINUE H(1)=0.0 DO 400 I=2,NR H(I)=H(I-1)+(F(I)+F(I-1))*DR/2.0 400 CONTINUE DO 500 I=1,NR G(I)=1.0-H(I)/H(NR) 500 CONTINUE XNU=G(NR1)/(2.0*DR*F(NR)) C C ********************************************************************* C WRITE(6,2500) INTR,XNU 2500 FORMAT(4X,'ITERATION NUMBER = ',I5,' NU = ',F12.4) C C ********************************************************************* C CHECK FOR ENDING C ********************************************************************* C IF(INTR.LT.5) GO TO 1000 DIFFX=ABS(XNU-XNUP) IF(DIFFX.LT.0.00000001) GO TO 4107 IF(INTR.GT.1000) GO TO 4107 XNUP=XNU GO TO 1000 4107 CONTINUE STOP END ************************************************************************ ***** ***** ***** PIPEFLOW ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR DEVELOPING LAMINAR PIPE FLOW ***** ***** ***** ***** THE VELOCITY AND TEMP. FIELDS ARE SIMULTANEOUSLY DEVELOPING. ***** ***** ***** ***** AN EXPLICIT FINITE DIFFERENCE TECHNIQUE IS USED ***** ***** ***** ***** EITHER THE WALL TEMP. OR THE WALL HEAT FLUX IS UNIFORM ***** ***** ***** *********************************************************************** * REAL U(2,100),V(2,100),T(2,100),A(100),B(100),C(100),R(100) * *********************************************************************** * OPEN(UNIT=1,FILE='PIPEFLPR.DAT') OPEN(UNIT=2,FILE='PIPEFLPL.DAT') * *********************************************************************** * WRITE(1,760) WRITE(6,760) * *********************** ENTER SPECIFIED VALUES ********************** * write(6,1400) write(6,1070) read(5,*) ZMAX write(6,1060) read(5,*) PR write(6,1050) read(5,*) N write(6,1600) read(5,*) IQT IPRF=150 IPPF=150 ********************************************************************* * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * IPRF=PRINT FREQUENCY * IPPF=PROFILE PRINT FREQUENCY * ZMAX=MAXIMUM Z TO WHICH SOLUTION IS TO BE UNDER TAKEN * IQT=BOUNDARY CONDITION IDENTIFIER (1=UNIFORM TEMP., 2=UNIFORM FLUX) ********************************************************************* * IF (IQT.EQ.1) THEN WRITE (1,1100) ELSE WRITE (1,1200) ENDIF WRITE (1,1000) N,PR,ZMAX * ********************************************************************** * IPRT=0 IPPT=0 N1=N-1 DR=0.5/N1 ZFAC=10.0*PR * ************************** ASSIGN INITIAL VALUES ********************* * Z=0.0 P1=-0.5 DO 10 J=1,N1 R(J)=(J-1)*DR U(1,J)=1.0/(1.0-2.0*DR) V(1,J)=0.0 T(1,J)=0.0 10 CONTINUE R(N)=0.5 U(1,N)=0.0 IF (IQT.EQ.1) THEN T(1,N)=1.0 ELSE T(1,N)=DY ENDIF * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * ***************** FIND Z-STEP TO ENSURE STABILITY ******************* * DZ=DR*DR*U(1,N1)/ZFAC * **************** CALCULATE THE PRESSURE ***************************** * A(1)=(P1/U(1,1)) & +(PR*DZ/U(1,1))*((2.0*U(1,2)-2.0*U(1,1))/(DR*DR)+GBUOY*T(1,1)) DO 20 J=2,N1 A(J)=(P1/U(1,J)) & +(PR*DZ/U(1,J))*(-V(1,J)*(U(1,J+1)-U(1,J-1))/(2.0*DR*PR) & +(U(1,J+1)+U(1,J-1)-2.0*U(1,J))/(DR*DR) & +(U(1,J+1)-U(1,J-1))/(2.0*R(J)*DR)+GBUOY*T(1,J)) 20 CONTINUE DO 25 J=2,N1 B(J)=(R(J)+R(J-1))*(DR/(4.0*DZ))*(1.0/U(1,J)+1.0/U(1,J-1)) C(J)=(R(J)+R(J-1))*(DR/(4.0*DZ))*(A(J)+A(J-1)) 25 CONTINUE B(N)=(R(N)+R(N1))*(DR/(4.0*DZ))*(1.0/U(1,N-1)) C(N)=(R(N)+R(N1))*(DR/(4.0*DZ))*A(N-1) BSUM=0.0 CSUM=0.0 DO 30 J=2,N BSUM=BSUM+B(J) CSUM=CSUM+C(J) 30 CONTINUE * P2=CSUM/BSUM * ********* CALCULATE STREAM WISE VELOCITY COMPONENT ******************* * DO 40 J=2,N1 U(2,J)=U(1,J)+A(J)-P2/U(1,J) 40 CONTINUE U(2,N)=0.0 U(2,1)=U(2,2) * **************** SOLVE FOR THE RADIAL VELOCITY COMPONENT *************** * DRDZ=DR/(2.0*DZ) V(2,1)=0.0 DO 50 J=2,N1 V(2,J)=(R(J-1)*V(2,J-1)+B(J)*P2-C(J))/R(J) 50 CONTINUE V(2,N)=0.0 * *********** SOLVE FOR TEMPERATURE FROM THE ENERGY EQUATION ************* * T(2,1)=0.0 DO 60 J=2,N1 T(2,J)=T(1,J)+(DZ/U(1,J))*(-V(1,J)*(T(1,J+1)-T(1,J-1))/(2.0*DR) * +(T(1,J+1)+T(1,J-1)-2.0*T(1,J))/(DR*DR) * +(T(1,J+1)-T(1,J-1))/(R(J)*2.0*DR)) 60 CONTINUE T(2,1)=T(2,2) IF (IQT.EQ.1) THEN T(2,N)=1.0 ELSE T(2,N)=T(2,N1)+DR ENDIF * *********************** CALCULATE NUSSELT NUMBER ********************** * IF (IQT.EQ.1) THEN XNUW=(T(2,N)-T(2,N1))/DR ELSE XNUW=1.0/T(2,N) ENDIF XNUMW=0.0 XMASS=0.0 DO 70 J=2,N XNUMW=XNUMW+(R(J)*U(2,J)*T(2,J)+R(J-1)*U(2,J-1)*T(2,J-1))*DR/2.0 XMASS=XMASS+(R(J)*U(2,J)+R(J-1)*U(2,J-1))*DR/2.0 70 CONTINUE IF (IQT.EQ.1) THEN XNUDM=XNUW/(1.0-XNUMW/XMASS) ELSE XNUDM=1.0/(T(2,N)-XNUMW/XMASS) ENDIF Z=Z+DZ XNUMW=XNUMW/Z * ********************* WRITE THE RESULTS ****************************** * IPRT=IPRT+1 IF(IPRT.LT.IPRF) GO TO 300 IPRT=0 IF (IQT.EQ.1) THEN WRITE(1,2000) Z,XNUW,XNUDM,XNUMW WRITE(2,2200) Z,XNUW,XNUDM,XNUMW,U(2,1),T(2,1) WRITE(6,2000) Z,XNUW,XNUDM,XNUMW ELSE WRITE(1,2100) Z,XNUW,XNUDM WRITE(2,2300) Z,XNUW,XNUDM,U(2,1),T(2,1) WRITE(6,2100) Z,XNUW,XNUDM ENDIF IPPT=IPPT+1 IF(IPPT.LT.IPPF) GO TO 300 IPPT=0 WRITE(1,3000) DO 80 J=1,N WRITE(1,4000) R(J),U(2,J),T(2,J),V(2,J) 80 CONTINUE 300 CONTINUE * ************************** CHECK FOR MAXIMUM X ************************ * IF(Z.GT.ZMAX) GO TO 200 * *************************** RETURN VALUES ****************************** * DO 90 J=1,N U(1,J)=U(2,J) V(1,J)=V(2,J) T(1,J)=T(2,J) 90 CONTINUE P1=P2 GO TO 100 * ****************************** ENDING ********************************* * 200 CONTINUE WRITE(1,3000) DO 86 J=1,N WRITE(1,4000) R(J),U(2,J),T(2,J),V(2,J) 86 CONTINUE STOP CLOSE(1) CLOSE(2) * ************************** FORMAT STATEMENTS ************************** * 760 FORMAT(//,3X,'LAMINAR DEVELOPING FLOW IN A PIPE', &/,3X,'_______________________________',//) 1000 FORMAT(6X,'NUMBER OF GRID POINTS = ',I4,/, $6X,'PRANDTL NUMBER = ' ,F10.5,/, $6X,'MAXIMUM Z = ',F10.6,//) 1050 format(//,5X,'Input the Number of Nodal Points ',/) 1060 format(//,5X,'Input the Prandtl Number ',/) 1070 format(//,5X,'Input the Maximum Z-value to be Considered ',/) 1100 FORMAT(6X,'UNIFORM WALL TEMPERATURE') 1200 FORMAT(6X,'UNIFORM WALL HEAT FLUX') 1400 format(////////,5X,'DEVELOPING LAMINAR PIPE FLOW', & /,5X,'--------------------------',///) 1600 format(//,5X,'Uniform Temp. (1) or Uniform Heat Flux (2) ',/) 2000 FORMAT(2X,'Z =',F8.6,2X, 'NU (TW-TIN)=',F7.4, $ 2X,'NU (TW-TMEAN)=',F7.4,2X,'NU (AVG. TO Z)=',F7.4) 2100 FORMAT(2X,'Z =',F8.6,3X, 'NU (TW-TIN)=' ,F9.4, $ 3X, 'NU (TW-TMEAN)=' ,F9.4) 2200 FORMAT(F8.6,',',F7.4,',',F7.4,',',F7.4,',',F7.4,',',F7.4) 2300 FORMAT(F8.6,',',F7.4,',',F7.4,',',F7.4,',',F7.4) 3000 FORMAT(' R U T V' ) 4000 FORMAT(4F12.5) * ********************************************************************* * END ************************************************************************ ***** ***** ***** PORCYL ***** ***** ______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR BOUNDARY LAYER FLOW OVER A CYLINDER ***** ***** ***** ***** IN A POROUS MEDIUM USING THE FINITE DIFFERENCE TECHNIQUE. ***** ***** ***** ***** A CONSTANT (I.E. UNIFORM) Y-GRID SPACING IS USED. ***** ***** ***** ***** THE WALL TEMPERATURE IS ASSUMED TO BE UNIFORM ***** ***** ***** ************************************************************************ C DIMENSION V(500),T(2,500),E(500),F(500) DIMENSION G(500),H(500),Y(500),A(500) REAL ND,NDM * ********************************************************************** * OPEN(1,FILE='PORCYPRT.DAT') OPEN(2,FILE='PORCYPLT.DAT') * ********************************************************************** * WRITE (1,1000) WRITE (1,1001) WRITE(6,6000) WRITE(6,1000) WRITE(6,1001) WRITE(6,6000) * SND=0.0 X = 0.0 N = 50 N1=N-1 Y(1) = 0.0 DY = 0.050 * DO 1 J = 2,N Y(J) = Y(J-1) + DY 1 CONTINUE DX = 0.001 DXMAX=0.005 XMAX=1.5707 * *************************** SOLUTION BEGINS ************************** * 100 CONTINUE * WRITE (1,2000) X * * CALL VELDV (X,U,DUX) * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * DO 12 J = 1,N V(J)= -Y(J)*DUX E(J) = (2.0/(DY*DY))+(U/DX) F(J) = V(J)/(2.0*DY)-1.0/(DY*DY) G(J) = -V(J)/(2.0*DY)-1.0/(DY*DY) H(J) = U*T(1,J)/DX 12 CONTINUE * E(1)=1.0 F(1)=0.0 G(1)=0.0 H(1)=1.0 E(N)=1.0 F(N)=0.0 G(N)=0.0 H(N)=0.0 * CALL TRISOL(N,E,F,G,H,A) * * DO 13 J = 1,N T(2,J) = A(J) 13 CONTINUE * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * TT = 0.01*T(2,1) DO 22 J=1,N1 IF (T(2,J).LE.TT) GO TO 23 22 CONTINUE 23 DE = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * IF (DE.LE.Y(N-5)) GO TO 24 NMIN = N+1 NMAX = N+10 IF (NMAX.GT.500) GO TO 101 DO 25 I=NMIN,NMAX Y(I) = Y(I-1) + DY T(2,I) = 0.0 25 CONTINUE WRITE (1,3000) WRITE (6,3000) N = NMAX 24 CONTINUE * * *************************** RETURN VALUES **************************** * DO 30 J=1,N T(1,J)=T(2,J) 30 CONTINUE * ************************** FIND HEAT TRANSFER ************************* * ND = (T(2,1)-T(2,2))/(Y(2)-Y(1)) IF( X.GT.0.0) THEN SND=SND+(ND+NDM)*DX/2.0 ENDIF NDM=ND * ************************** WRITE THE RESULTS ************************** * WRITE(1,4000) ND,DE WRITE(6,4200) X,ND WRITE(2,4100) X,ND,DE WRITE(1,4001) DO 40 J=1,N WRITE(1,4002)J,Y(J),T(2,J) 40 CONTINUE * ************************** CHANGE DX ********************************** * IF (DX.LT.DXMAX) THEN DX = 1.05*DX ELSE DX = DXMAX ENDIF * X = X + DX * ************************** CHECK FOR MAXIMUM X ************************ * IF(X.GT.XMAX) GO TO 102 GO TO 100 101 WRITE (1,5000) WRITE (6,5000) GO TO 103 102 WRITE(1,5001) WRITE(6,5001) WRITE(1,5500) SND/XMAX WRITE(6,5500) SND/XMAX 103 CONTINUE CLOSE(1) CLOSE(2) STOP * ************************** FORMAT STATEMENTS ************************** * 1000 FORMAT (10X,'BOUNDARY LAYER FLOW OVER AN ISOTHERMAL CYLINDER') 1001 FORMAT (10X,'-----------------------------------------',//) 2000 FORMAT (/,5X,'SOLUTION FOR X =',F12.4,/ & ,5X,'___________________________',/) 3000 FORMAT (5X,'****** NUMBER OF GRID POINTS INCREASED ******') 4000 FORMAT (5X,'ND/PD 1/2 =',E12.4,' DELT =',E12.4,/) 4001 FORMAT (2X,'J',7X,'Y',12X,'T') 4002 FORMAT (I4,2E12.4) 4100 FORMAT (F10.5,',',F10.5,',',F10.5) 4200 FORMAT (5X,' X =',F10.4,' ND/PD 1/2 =',F10.4) 5000 FORMAT (5X,'N EXCEEDS 500') 5001 FORMAT (/,5X,'********** XMAX REACHED **********') 5500 FORMAT (5X,' MEAN ND/PD 1/2 =',F10.4) 6000 FORMAT(////////) END * * ******************************************************************** * SUBROUTINE TRISOL(NN,A,B,C,D,H) * *********** THIS IS A TRI-DIAGONAL MATRIX SOLVER ******************* * --- --- * * THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM * DIMENSION A(500),B(500),C(500),D(500),H(500),W(500),Q(500),G(500) W(1)=A(1) G(1)=D(1)/W(1) DO 2 K=2,NN K1=K-1 Q(K1)=B(K1)/W(K1) W(K)=A(K)-C(K)*Q(K1) G(K)=(D(K)-C(K)*G(K1))/W(K) 2 CONTINUE H(NN)=G(NN) N1=NN-1 DO 3 K=1,N1 KK=NN-K H(KK)=G(KK)-Q(KK)*H(KK+1) 3 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE VELDV (X,U,DUX) * ****** THIS DETERMINES X-WISE VELOCITY AND VELOCITY GRADIENT ****** * U = 2.0*SIN(2.0*X) DUX = 4.0*COS(2.0*X) RETURN END C C ********************************************************************* C C PROGRAM "RECDUCT" C C THIS PROGRAM SOLVES FOR THE HEAT TRANSFER RATE IN FULLY- C DEVELOPED LAMINAR FLOW IN A RECTANGULAR DUCT. C C ********************************************************************* C REAL U(91,91),T(91,91),A(91),B(91),C(91), $ D(91),Q(91),H(91),QL(91),SUM1(91),SUM2(91),QTOP(91),QSID(91) C C ******************************************************************** C OPEN (1,FILE='RECDUCPR.DAT') OPEN (3,FILE='RECDUCPL.DAT') C C ********************************************************************* C INPUT C ********************************************************************* WRITE(6,6000) 6000 FORMAT(////////) WRITE(6,6001) 6001 FORMAT (10X,'FULLY-DEVELOPED FLOW IN A RECTANGULAR DUCT') WRITE(6,6002) 6002 FORMAT (10X,'------------------------------------------',//) WRITE(6,6003) 6003 FORMAT(/,' INPUT THE DUCT ASPECT RATIO',///) READ(5,*) AR WRITE(6,6004) 6004 FORMAT(/,' INPUT THE NUMBER OF NODAL POINTS', & /,' IN THE X- AND Y -DIRECTIONS ',///) READ(5,*) NX,NY C C ********************************************************************* WRITE(1,5247) 5247 FORMAT(/,5X,' FLOW IN A RECTANGULAR DUCT') WRITE(1,5237) AR 5237 FORMAT(/,' ASPECT RATIO= ', F12.3,//) C ********************************************************************* C C NX=31 C NY=31 DX=AR/(2.0*(NX-1)) DY=0.5/(NY-1) REU=1.0 C C ********************************************************************* C DX2=DX*DX DY2=DY*DY DX2I=1.0/DX2 DY2I=1.0/DY2 NX1=NX-1 NY1=NY-1 C C ********************************************************************* WRITE(3,7237) AR 7237 FORMAT(F12.4) C ********************************************************************* C C ********************************************************************* C INITIAL GUESSED VALUES C ********************************************************************* C DO 100 I=1,NX X=(I-1)*DX/(0.5*(NX-1)) DO 100 J=1,NY Y=(J-1)*2.0*DY/(AR*(NY-1)) U(I,J)=(1.0-X*X)*(1.0-Y*Y) 100 CONTINUE C C ********************************************************************* C SOLVE FOR VELOCITY DISTRIBUTION C ********************************************************************* C INTR=0 UCENX=1.0 1000 CONTINUE C C ********************************************************************* C INTR=INTR+1 C C C ********************************************************************* C FIND VELOCITY ON Y-DIRECTION LINES C ********************************************************************* C DO 5700 I=2,NX1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 5600 J=2,NY1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(U(I+1,J)+U(I-1,J))-1.0 5600 CONTINUE A(NY)=1.0 B(NY)=0.0 C(NY)=0.0 D(NY)=0.0 CALL TRISOL(NY,A,B,C,D,H) DO 5800 J=1,NY U(I,J)=U(I,J)+REU*(H(J)-U(I,J)) 5800 CONTINUE 5700 CONTINUE DO 5900 J=1,NY U(1,J)=U(2,J) U(NX,J)=0.0 5900 CONTINUE C C ********************************************************************* C FIND VELOCITY ON X-DIRECTION LINES C ********************************************************************* C DIFFX=0.0 DO 6110 J=2,NY1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 7001 I=2,NX1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(U(I,J+1)+U(I,J-1))-1.0 7001 CONTINUE A(NX)=1.0 B(NX)=0.0 C(NX)=0.0 D(NX)=0.0 CALL TRISOL(NX,A,B,C,D,H) DO 6211 I=1,NX U(I,J)=U(I,J)+REU*(H(I)-U(I,J)) 6211 CONTINUE 6110 CONTINUE DO 6250 I=1,NX U(I,1)=U(I,2) U(I,NY)=0.0 6250 CONTINUE U(1,1)=(2.0*U(2,1)*DX2I+2.0*U(1,2)*DY2I+1.0)/(2*DX2I+2.0*DY2I) C C ********************************************************************* C WRITE(6,2500) INTR 2500 FORMAT(4X,'VELOCITY ITERATION NUMBER = ',I5) C C ********************************************************************* C CHECK FOR ENDING C ********************************************************************* C IF(INTR.LT.1000) GO TO 1000 DIFFX=ABS(UCEN-UCENX) IF(DIFFX.LT.0.000001) GO TO 4107 IF(INTR.GT.3000) GO TO 4107 UCENX=UCEN GO TO 1000 4107 CONTINUE C C ********************************************************************* C UCEN=U(1,1) DO 1777 I=1,NX DO 1777 J=1,NY U(I,J)=U(I,J)/UCEN WRITE(1,1778) I,J,U(I,J) 1778 FORMAT(' I = ',I4,' J = ',I4,' U(I,J) = ',F12.3) X=(I-1)*DX/(0.5*(NX-1)) Y=(J-1)*2.0*DY/(AR*(NY-1)) WRITE(3,1779) X,Y,U(I,J) 1779 FORMAT(F12.5,',',F12.5,',',F12.5) 1777 CONTINUE C C ********************************************************************* C SOLVE FOR TEMPERATURE DISTRIBUTION C ********************************************************************* C DO 1788 I=1,NX DO 1788 J=1,NY T(I,J)=U(I,J) 1788 CONTINUE INTR=0 TCENX=1.0 1010 CONTINUE C C ********************************************************************* C INTR=INTR+1 C C C ********************************************************************* C FIND TEMPERATURE ON Y-DIRECTION LINES C ********************************************************************* C DO 5710 I=2,NX1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 5610 J=2,NY1 A(J)=-2.0*DX2I-2.0*DY2I+U(I,J) B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(T(I+1,J)+T(I-1,J)) 5610 CONTINUE A(NY)=1.0 B(NY)=0.0 C(NY)=0.0 D(NY)=0.0 CALL TRISOL(NY,A,B,C,D,H) DO 5810 J=1,NY T(I,J)=T(I,J)+REU*(H(J)-T(I,J)) 5810 CONTINUE 5710 CONTINUE DO 5910 J=1,NY T(1,J)=T(2,J) T(NX,J)=0.0 5910 CONTINUE C C ********************************************************************* C FIND TEMPERATURE ON X-DIRECTION LINES C ********************************************************************* C DIFFX=0.0 DO 6120 J=2,NY1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 7011 I=2,NX1 A(I)=-2.0*DX2I-2.0*DY2I+U(I,J) B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(T(I,J+1)+T(I,J-1)) 7011 CONTINUE A(NX)=1.0 B(NX)=0.0 C(NX)=0.0 D(NX)=0.0 CALL TRISOL(NX,A,B,C,D,H) DO 6221 I=1,NX T(I,J)=T(I,J)+REU*(H(I)-T(I,J)) 6221 CONTINUE 6120 CONTINUE DO 6260 I=1,NX T(I,1)=T(I,2) T(I,NY)=0.0 6260 CONTINUE T(1,1)=(2.0*T(2,1)*DX2I+2.0*T(1,2)*DY2I)/(2*DX2I+2.0*DY2I-1.0) C C ********************************************************************* C WRITE(6,2501) INTR 2501 FORMAT(4X,'TEMPERATURE ITERATION NUMBER = ',I5) C C ********************************************************************* C CHECK FOR ENDING C ********************************************************************* C IF(INTR.LT.1000) GO TO 1010 DIFFX=ABS(TCEN-TCENX) IF(DIFFX.LT.0.000001) GO TO 4117 TCENX=TCEN IF(INTR.GT.3000) GO TO 4117 GO TO 1010 4117 CONTINUE TCEN=T(1,1) DO 7777 I=1,NX DO 7777 J=1,NY T(I,J)=T(I,J)/TCEN 7777 CONTINUE DO 1889 I=1,NX SUM1(I)=0.0 SUM2(I)=0.0 DO 1889 J=2,NY SUM1(I)=SUM1(I)+0.5*(U(I,J)*T(I,J)+U(I,J-1)*T(I,J-1))*DY SUM2(I)=SUM2(I)+0.5*(U(I,J)+U(I,J-1))*DY 1889 CONTINUE SU1=0.0 SU2=0.0 DO 1890 I=2,NX SU1=SU1+0.5*(SUM1(I)+SUM1(I-1))*DX SU2=SU2+0.5*(SUM2(I)+SUM2(I-1))*DX 1890 CONTINUE TMEAN=SU1/SU2 DO 2889 I=1,NX QTOP(I)=(T(I,NY)+T(I,NY1))/DY 2889 CONTINUE DO 3889 J=1,NY QSID(J)=(T(NX,J)+T(NX1,J))/DX 3889 CONTINUE QSUM=0.0 DO 4889 I=2,NX QSUM=QSUM+0.5*(QTOP(I)+QTOP(I-1))*DX 4889 CONTINUE DO 5889 J=2,NY QSUM=QSUM+0.5*(QSID(J)+QSID(J-1))*DY 5889 CONTINUE XNU=(QSUM/(0.5+AR/2.0))/TMEAN C C ********************************************************************* C DO 1787 I=1,NX DO 1787 J=1,NY WRITE(1,1789) I,J,T(I,J) 1789 FORMAT(' I = ',I4,' J = ',I4,' T(I,J) = ',F12.3) X=(I-1)*DX/(0.5*(NX-1)) Y=(J-1)*2.0*DY/(AR*(NY-1)) WRITE(3,1779) X,Y,T(I,J) 1787 CONTINUE WRITE(1,1989) TMEAN WRITE(6,1989) TMEAN 1989 FORMAT(' THETA MEAN = ',F12.5) WRITE(1,6989) XNU WRITE(6,6989) XNU 6989 FORMAT(' MEAN NU = ',F12.5) C C ********************************************************************* C STOP END C C ********************************************************************* C END OF MAIN PROGRAM C ********************************************************************* C C C ********************************************************************* C SUBROUTINE TRISOL(N,A,B,C,D,H) C C ********************************************************************* C TRI-DIAGONAL MATRIX SOLVER C ********************************************************************* C REAL A(91),B(91),C(91),D(91),H(91),W(91),R(91),G(91) W(1)=A(1) G(1)=D(1)/W(1) DO 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 CONTINUE H(N)=G(N) N1=N-1 DO 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 CONTINUE RETURN END ************************************************************************ ***** ***** ***** SIMPLATE ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES SIMILARITY SOLUTION RESULTS FOR ***** ***** ***** ***** LAMINAR BOUNDARY LAYER FLOW OVER AN ISOTHERMAL FLAT PLATE ***** ***** ***** ************************************************************************ * DIMENSION T(1500),F(1500),G(1500),H(1500) & ,HGUESS(1500),FPIN(1500),P(1500) * *********************************************************************** * OPEN(UNIT=1,FILE='SIMPLAPR.DAT') OPEN(UNIT=2,FILE='SIMPLAPL.DAT') * *********************************************************************** * WRITE(6,4250) WRITE(6,4290) WRITE(6,4260) READ(5,*) PR WRITE(6,4250) DETA=0.01 N=1500 WRITE(1,4290) WRITE(1,4270) PR * *********************************************************************** * * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * INTR=1 HGUESS(INTR)=1.0 100 F(1)=0.0 G(1)=0.0 H(1)=HGUESS(INTR) DO 1000 I=2,N F(I)=F(I-1)+G(I-1)*DETA G(I)=G(I-1)+H(I-1)*DETA H(I)=H(I-1)-F(I-1)*H(I-1)*DETA/2.0 1000 CONTINUE FPIN(INTR)=G(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=HGUESS(1)-0.000001 GO TO 100 ELSE IF(INTR.GT.200) GO TO 200 IF(ABS(1.0-FPIN(INTR)).LT.0.0000005) GO TO 300 INTR=INTR+1 DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+DGUESS GO TO 100 ENDIF 200 WRITE(6,2000) GO TO 7777 300 WRITE(6,3000) INTR 7777 CONTINUE DO 5000 I=1,N P(I)=H(I)**PR 5000 CONTINUE T(1)=0.0 DO 6000 I=2,N T(I)=T(I-1)+DETA*(P(I)+P(I-1))/2.0 6000 CONTINUE WRITE(1,4500) DO 7000 I=1,N T(I)=T(I)/T(N) ETA=(I-1)*DETA WRITE(1,4000) ETA,F(I),G(I),H(I),T(I) WRITE(2,5500) ETA,G(I),T(I) 7000 CONTINUE TGRAD=(T(2)-T(1))/DETA WRITE(1,4280) TGRAD STOP CLOSE(1) CLOSE(2) 2000 FORMAT(' FAILURE TO CONVERGE') 3000 FORMAT(' CONVERGENCE IN ',I5,' ITERATIONS') 4000 FORMAT(5F10.6) 4250 FORMAT(////) 4260 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 4270 FORMAT(//,' PRANDTL NUMBER = ',F12.3,//) 4280 FORMAT(//,' THETA GRADIENT AT WALL = ',F12.5) 4290 FORMAT(/,' SIMILARITY SOLUTION FOR FLOW OVER ISOTHERMAL PLATE',//, $ ' Output written to files SIMPLAPR.DAT and SIMPLAPL.DAT', //) 4500 FORMAT(' ETA F dF/DETA D2F/DETA2 THETA') 5500 FORMAT(F10.6,',',F10.6,',',F10.6) END ************************************************************************ ***** ***** ***** SIMPLCOM ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES FIRST ORDER SERIES SIMILARITY TYPE ***** ***** ***** ***** SOLUTION RESULTS FOR COMBINED CONVECTIVE ***** ***** ***** ***** LAMINAR BOUNDARY LAYER FLOW OVER AN ISOTHERMAL FLAT PLATE ***** ***** ***** ************************************************************************ * *********************************************************************** * DIMENSION TO(1500),FO(1500),GO(1500),HO(1500) & ,T1(1500),F1(1500),G1(1500),H1(1500),HT(1500) & ,HGUESS(1500),FPIN(1500),P(1500) REAL KF0,KG0,KH0,KF1,KG1,KH1,KF2,KG2,KH2,KF3,KG3,KH3 & ,KT0,KT1,KT2,KT3 COMMON DETA,PR * *********************************************************************** * OPEN(UNIT=1,FILE='SIMPLCPR.DAT') OPEN(UNIT=2,FILE='SIMPLCPL.DAT') * *********************************************************************** * WRITE(6,4250) WRITE(6,4290) WRITE(6,4260) READ(5,*) PR WRITE(6,4360) READ(5,*) DIR IF (DIR.EQ.2) THEN DIRM=-1 ELSE DIRM=1 ENDIF WRITE(6,4250) DETA=0.01 N=1500 WRITE(1,4290) WRITE(1,4270) PR * *********************************************************************** * * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * * FIRST ORDER ( FORCED CONVECTION ) SOLUTION * *********************************************************************** * INTR=1 HGUESS(INTR)=1.0 100 FO(1)=0.0 GO(1)=0.0 HO(1)=HGUESS(INTR) DO 1000 I=2,N F0=FO(I-1) G0=GO(I-1) H0=HO(I-1) FC=F0 GC=G0 HC=H0 CALL RUTKUN(FC,GC,HC,KF0,KG0,KH0) FC=F0+0.5*KF0 GC=G0+0.5*KG0 HC=H0+0.5*KH0 CALL RUTKUN(FC,GC,HC,KF1,KG1,KH1) FC=F0+0.5*KF1 GC=G0+0.5*KG1 HC=H0+0.5*KH1 CALL RUTKUN(FC,GC,HC,KF2,KG2,KH2) FC=F0+KF2 GC=G0+KG2 HC=H0+KH2 CALL RUTKUN(FC,GC,HC,KF3,KG3,KH3) FO(I)=FO(I-1)+(KF0+2.0*KF1+2.0*KF2+KF3)/6.0 GO(I)=GO(I-1)+(KG0+2.0*KG1+2.0*KG2+KG3)/6.0 HO(I)=HO(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 1000 CONTINUE FPIN(INTR)=GO(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=0.9*HGUESS(1) GO TO 100 ELSE IF(INTR.GT.200) GO TO 200 IF(ABS(1.0-FPIN(INTR)).LT.0.0000005) GO TO 300 INTR=INTR+1 DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+0.2*DGUESS GO TO 100 ENDIF 200 WRITE(6,2000) GO TO 7777 300 WRITE(6,3000) INTR 7777 CONTINUE DO 5000 I=1,N P(I)=HO(I)**PR 5000 CONTINUE TO(1)=0.0 DO 6000 I=2,N TO(I)=TO(I-1)+DETA*(P(I)+P(I-1))/2.0 6000 CONTINUE WRITE(1,4500) DO 7000 I=1,N TO(I)=TO(I)/TO(N) ETA=(I-1)*DETA WRITE(1,4000) ETA,FO(I),GO(I),HO(I),TO(I) WRITE(2,5500) ETA,GO(I),TO(I) 7000 CONTINUE TGRAD=(TO(1)-TO(2))/DETA * *********************************************************************** * * SECOND ORDER VELOCITY SOLUTION * *********************************************************************** * INTR=1 HGUESS(INTR)=0.1 110 F1(1)=0.0 G1(1)=0.0 H1(1)=HGUESS(INTR) DO 1010 I=2,N F0=F1(I-1) G0=G1(I-1) H0=H1(I-1) TL=DIRM*(1.0-(TO(I)+TO(I-1))/2.0) FL=(FO(I)+FO(I-1))/2.0 HL=(HO(I)+HO(I-1))/2.0 FC=F0 GC=G0 HC=H0 CALL RUTK1N(FC,GC,HC,KF0,KG0,KH0,FL,HL,TL) FC=F0+0.5*KF0 GC=G0+0.5*KG0 HC=H0+0.5*KH0 CALL RUTK1N(FC,GC,HC,KF1,KG1,KH1,FL,HL,TL) FC=F0+0.5*KF1 GC=G0+0.5*KG1 HC=H0+0.5*KH1 CALL RUTK1N(FC,GC,HC,KF2,KG2,KH2,FL,HL,TL) FC=F0+KF2 GC=G0+KG2 HC=H0+KH2 CALL RUTK1N(FC,GC,HC,KF3,KG3,KH3,FL,HL,TL) F1(I)=F1(I-1)+(KF0+2.0*KF1+2.0*KF2+KF3)/6.0 G1(I)=G1(I-1)+(KG0+2.0*KG1+2.0*KG2+KG3)/6.0 H1(I)=H1(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 1010 CONTINUE FPIN(INTR)=G1(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=0.9*HGUESS(1) GO TO 110 ELSE IF(INTR.GT.200) GO TO 210 IF(ABS(FPIN(INTR)).LT.0.0000005) GO TO 310 INTR=INTR+1 DGUESS=+(FPIN(INTR-1))*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+0.2*DGUESS GO TO 110 ENDIF 210 WRITE(6,2000) GO TO 7787 310 WRITE(6,3010) INTR 7787 CONTINUE * *********************************************************************** * * SECOND ORDER TEMPERATURE SOLUTION * *********************************************************************** * INTR=1 HGUESS(INTR)=0.1 130 T1(1)=0.0 HT(1)=HGUESS(INTR) DO 1030 I=2,N T0=T1(I-1) H0=HT(I-1) TL=-(TO(I)-TO(I-1))/DETA FL=(FO(I)+FO(I-1))/2.0 FU=(F1(I)+F1(I-1))/2.0 TC=T0 HC=H0 CALL RUTKTN(HC,KT0,KH0,FL,FU,TL) HC=H0+0.5*KH0 CALL RUTKTN(HC,KT1,KH1,FL,FU,TL) HC=H0+0.5*KH1 CALL RUTKTN(HC,KT2,KH2,FL,FU,TL) HC=H0+KH2 CALL RUTKTN(HC,KT3,KH3,FL,FU,TL) T1(I)=T1(I-1)+(KT0+2.0*KT1+2.0*KT2+KT3)/6.0 HT(I)=HT(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 1030 CONTINUE FPIN(INTR)=T1(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=0.9*HGUESS(1) GO TO 130 ELSE IF(INTR.GT.200) GO TO 430 IF(ABS(FPIN(INTR)).LT.0.0000005) GO TO 410 INTR=INTR+1 DGUESS=+(FPIN(INTR-1))*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+0.2*DGUESS GO TO 130 ENDIF 430 WRITE(6,2000) GO TO 7797 410 WRITE(6,3020) INTR 7797 CONTINUE * *********************************************************************** * WRITE(1,4510) DO 7030 I=1,N ETA=(I-1)*DETA WRITE(1,4000) ETA,F1(I),G1(I),HT(I),T1(I) WRITE(2,5500) ETA,G1(I),T1(I) 7030 CONTINUE WRITE(1,4281) HO(1),H1(1) WRITE(6,4281) HO(1),H1(1) WRITE(1,4282) TGRAD,HT(1) WRITE(6,4282) TGRAD,HT(1) * *********************************************************************** * CLOSE(1) CLOSE(2) STOP 2000 FORMAT(' ******** FAILURE TO CONVERGE *********') 3000 FORMAT(' 1 ST. ORDER VELOCITY FUNC. CONVERGENCE IN ',I5,' ITERS') 3010 FORMAT(' 2 ND. ORDER VELOCITY FUNC. CONVERGENCE IN ',I5,' ITERS') 3020 FORMAT(' 2 ND. ORDER TEMP. FUNC. CONVERGENCE IN ',I5,' ITERS') 4000 FORMAT(5F10.6) 4250 FORMAT(////) 4260 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 4360 FORMAT(' INPUT ASSISTING FLOW (1) OR OPPOSING FLOW (2) ') 4270 FORMAT(//,' PRANDTL NUMBER = ',F12.3,//) 4280 FORMAT(//,' THETA GRADIENT AT WALL = ',F9.3) 4281 FORMAT(//,' FIRST ORDER VEL. GRADIENT AT WALL = ',F9.3, & /,' SECOND ORDER VEL. GRADIENT AT WALL = ',F9.3) 4282 FORMAT(//,' FIRST ORDER THETA GRADIENT AT WALL = ',F9.3, & /,' SECOND ORDER THETA GRADIENT AT WALL = ',F9.3) 4290 FORMAT(/,' FIRST ORDER COMBINED CONVECTION OVER ISOTHERMAL PLATE', $ /,' __________________________________________________',//) 4500 FORMAT(' ETA F dF/DETA D2F/DETA2 THETA') 4510 FORMAT(' ETA F dF/DETA DTHETA/DETA THETA') 5500 FORMAT(F10.6,',',F10.6,',',F10.6) END * *********************************************************************** * * 0 ORDER RUTTA-KUNGE SOLUTION ROUTINE * *********************************************************************** * SUBROUTINE RUTKUN(F,G,H,KF,KG,KH) REAL KF,KG,KH COMMON DETA,PR KF=DETA*G KG=DETA*H KH=-0.5*F*H*DETA RETURN END * *********************************************************************** * * 1 ST ORDER RUTTA-KUNGE SOLUTION ROUTINE * *********************************************************************** * SUBROUTINE RUTK1N(F,G,H,KF,KG,KH,FL,HL,TL) REAL KF,KG,KH COMMON DETA,PR KF=DETA*G KG=DETA*H KH=-0.5*FL*H*DETA-0.5*F*HL*DETA-TL*DETA RETURN END * *********************************************************************** * * 1 ST ORDER TEMPERATURE RUTTA-KUNGE SOLUTION ROUTINE * *********************************************************************** * SUBROUTINE RUTKTN(H,KT,KH,FL,FU,TL) REAL KF,KT,KH COMMON DETA,PR KT=DETA*H KH=-0.5*PR*(FU*TL+FL*H)*DETA RETURN END ************************************************************************ ***** ***** ***** SIMPLNAT ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES SIMILARITY SOLUTION RESULTS FOR ***** ***** ***** ***** LAMINAR FREE CONVECTIVE BOUNDARY LAYER FLOW OVER A ***** ***** ***** ***** VERTICAL ISOTHERMAL FLAT PLATE. ***** ***** ***** ************************************************************************ * DIMENSION T(5000),F(5000),G(5000),H(5000),A(5000) & ,HGUESS(3),AGUESS(3),FPIN(3),TPIN(3) REAL KF0,KG0,KH0,KF1,KG1,KH1,KF2,KG2,KH2,KF3,KG3,KH3 & ,KT0,KA0,KT1,KA1,KT2,KA2,KT3,KA3 COMMON DETA,PR * *********************************************************************** * OPEN(UNIT=1,FILE='SIMPLNPR.DAT') OPEN(UNIT=2,FILE='SIMPLNPL.DAT') * *********************************************************************** * WRITE(6,4250) WRITE(6,4290) WRITE(6,4260) READ(5,*) PR WRITE(6,4250) DETA=0.0015 N=5000 WRITE(1,4290) WRITE(1,4270) PR * *********************************************************************** * * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * INTR=1 HGUESS(1)=0.45 AGUESS(1)=-1.1 IF(PR.LT.15) THEN HGUESS(1)=0.6 AGUESS(1)=-0.9 ENDIF IF(PR.LT.10) THEN HGUESS(1)=0.7 AGUESS(1)=-0.7 ENDIF IF(PR.LT.5) THEN HGUESS(1)=0.8 AGUESS(1)=-0.5 ENDIF IF(PR.LT.1) THEN HGUESS(1)=0.9 AGUESS(1)=-0.4 ENDIF 100 CONTINUE DO 1001 J=1,3 F(1)=0.0 G(1)=0.0 T(1)=1.0 H(1)=HGUESS(J) A(1)=AGUESS(J) DO 1000 I=2,N F0=F(I-1) G0=G(I-1) H0=H(I-1) A0=A(I-1) T0=T(I-1) FC=F0 GC=G0 HC=H0 AC=A0 TC=T0 CALL RUTKUN(FC,GC,HC,TC,KF0,KG0,KH0) CALL RUTTUN(TC,AC,FC,KT0,KA0) FC=F0+0.5*KF0 GC=G0+0.5*KG0 HC=H0+0.5*KH0 TC=T0+0.5*KT0 AC=A0+0.5*KA0 CALL RUTKUN(FC,GC,HC,TC,KF1,KG1,KH1) CALL RUTTUN(TC,AC,FC,KT1,KA1) FC=F0+0.5*KF1 GC=G0+0.5*KG1 HC=H0+0.5*KH1 TC=T0+0.5*KT1 AC=A0+0.5*KA1 CALL RUTKUN(FC,GC,HC,TC,KF2,KG2,KH2) CALL RUTTUN(TC,AC,FC,KT2,KA2) FC=F0+KF2 GC=G0+KG2 HC=H0+KH2 TC=T0+KT2 AC=A0+KA2 CALL RUTKUN(FC,GC,HC,TC,KF3,KG3,KH3) CALL RUTTUN(TC,AC,FC,KT3,KA3) F(I)=F(I-1)+(KF0+2.0*KF1+2.0*KF2+KF3)/6.0 G(I)=G(I-1)+(KG0+2.0*KG1+2.0*KG2+KG3)/6.0 H(I)=H(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 T(I)=T(I-1)+(KT0+2.0*KT1+2.0*KT2+KT3)/6.0 A(I)=A(I-1)+(KA0+2.0*KA1+2.0*KA2+KA3)/6.0 1000 CONTINUE FPIN(J)=G(N) TPIN(J)=T(N) IF(ABS(FPIN(1)).LT.0.0000005) GO TO 300 IF (J.EQ.1) THEN HGUESS(2)=HGUESS(1)+0.001 AGUESS(2)=AGUESS(1) ENDIF IF (J.EQ.2) THEN HGUESS(3)=HGUESS(1) AGUESS(3)=AGUESS(1)+0.001 ENDIF 1001 CONTINUE IF(INTR.GT.100) GO TO 200 DEFVF=(FPIN(2)-FPIN(1))/0.001 DEFVT=(FPIN(3)-FPIN(1))/0.001 DETVF=(TPIN(2)-TPIN(1))/0.001 DETVT=(TPIN(3)-TPIN(1))/0.001 DH=(TPIN(1)/DETVT - FPIN(1)/DEFVT)/(DETVF/DETVT-DEFVF/DEFVT) DA=(TPIN(1)/DETVF - FPIN(1)/DEFVF)/(DETVT/DETVF-DEFVT/DEFVF) write(6,9999) INTR,FPIN(1),TPIN(1) 9999 format(2x' Iter. No. = ',I4,' F prime infin. = ', F10.5, & ' T infin. = ',F10.5) HGUESS(1)=HGUESS(1)-0.5*DH AGUESS(1)=AGUESS(1)-0.5*DA INTR=INTR+1 GO TO 100 200 WRITE(6,2000) GO TO 7777 300 WRITE(6,3000) INTR 7777 CONTINUE WRITE(1,4500) DO 7000 I=1,N,10 ETA=(I-1)*DETA WRITE(1,4000) ETA,F(I),G(I),H(I),T(I) WRITE(2,5500) ETA,G(I),T(I) 7000 CONTINUE WRITE(1,4280) A(1) WRITE(6,4280) A(1) WRITE(1,4281) H(1) WRITE(6,4281) H(1) CLOSE(1) CLOSE(2) STOP 2000 FORMAT(' ******** FAILURE TO CONVERGE *********') 3000 FORMAT(' VELOCITY FUNCTION CONVERGENCE IN ',I5,'ITERATIONS') 4000 FORMAT(5F10.6) 4250 FORMAT(////) 4260 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 4270 FORMAT(//,' PRANDTL NUMBER = ',F12.3,//) 4280 FORMAT(//,' THETA GRADIENT AT WALL = ',F9.3) 4281 FORMAT(//,' F PRIME GRADIENT AT WALL = ',F9.3) 4290 FORMAT(/,'SIMILARITY SOLUTION FOR FREE CONVECTIVE FLOW',/ $ ' OVER A VERTICAL ISOTHERMAL PLATE',/, $ ' ____________________________________________',//) 4500 FORMAT(' ETA F dF/DETA D2F/DETA2 THETA') 5500 FORMAT(F10.6,',',F10.6,',',F10.6) END * *********************************************************************** * * RUTTA-KUNGE SOLUTION ROUTINE FOR VELOCITY * *********************************************************************** * SUBROUTINE RUTKUN(F,G,H,T,KF,KG,KH) REAL KF,KG,KH COMMON DETA,PR KF=DETA*G KG=DETA*H KH=(-0.75*F*H+G*G/2.0-T)*DETA RETURN END * *********************************************************************** * * RUTTA-KUNGE SOLUTION ROUTINE FOR TEMPERATURE * *********************************************************************** * SUBROUTINE RUTTUN(T,A,F,KT,KA) REAL KT,KA COMMON DETA,PR KT=DETA*A KA=-0.75*PR*F*A*DETA RETURN END ************************************************************************ ***** ***** ***** SIMVART ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES SIMILARITY SOLUTION RESULTS FOR ***** ***** ***** ***** LAMINAR BOUNDARY LAYER FLOW OVER A FLAT PLATE WHOSE ***** ***** ***** ***** SURFACE TEMPERATURE IS ( T1 + C x**n ) ***** ***** ***** ************************************************************************ * *********************************************************************** * DIMENSION T(1500),F(1500),G(1500),H(1500) & ,HGUESS(1500),FPIN(1500),P(1500) REAL KF0,KG0,KH0,KF1,KG1,KH1,KF2,KG2,KH2,KF3,KG3,KH3 & ,KT0,KT1,KT2,KT3 COMMON DETA,XN,PR * *********************************************************************** * OPEN(UNIT=1,FILE='SIMVARPR.DAT') OPEN(UNIT=2,FILE='SIMVARPL.DAT') * *********************************************************************** * WRITE(6,4250) WRITE(6,4290) WRITE(6,4260) READ(5,*) PR WRITE(6,4250) WRITE(6,4265) READ(5,*) XN WRITE(6,4250) DETA=0.01 N=1500 WRITE(1,4290) WRITE(1,4270) PR * *********************************************************************** * * N=NUMBER OF GRID POINTS * PR=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * INTR=1 HGUESS(INTR)=1.0 100 F(1)=0.0 G(1)=0.0 H(1)=HGUESS(INTR) DO 1000 I=2,N F0=F(I-1) G0=G(I-1) H0=H(I-1) FC=F0 GC=G0 HC=H0 CALL RUTKUN(FC,GC,HC,KF0,KG0,KH0) FC=F0+0.5*KF0 GC=G0+0.5*KG0 HC=H0+0.5*KH0 CALL RUTKUN(FC,GC,HC,KF1,KG1,KH1) FC=F0+0.5*KF1 GC=G0+0.5*KG1 HC=H0+0.5*KH1 CALL RUTKUN(FC,GC,HC,KF2,KG2,KH2) FC=F0+KF2 GC=G0+KG2 HC=H0+KH2 CALL RUTKUN(FC,GC,HC,KF3,KG3,KH3) F(I)=F(I-1)+(KF0+2.0*KF1+2.0*KF2+KF3)/6.0 G(I)=G(I-1)+(KG0+2.0*KG1+2.0*KG2+KG3)/6.0 H(I)=H(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 1000 CONTINUE FPIN(INTR)=G(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=HGUESS(1)-0.0001 GO TO 100 ELSE IF(INTR.GT.200) GO TO 200 IF(ABS(1.0-FPIN(INTR)).LT.0.0000005) GO TO 300 INTR=INTR+1 DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+DGUESS GO TO 100 ENDIF 200 WRITE(6,2000) GO TO 7777 300 WRITE(6,3000) INTR 7777 CONTINUE HU=H(1) INTR=1 HGUESS(INTR)=G(1) 101 T(1)=0.0 H(1)=HGUESS(INTR) DO 1001 I=2,N FC=(F(I)+F(I-1))/2.0 GC=(G(I)+G(I-1))/2.0 T0=T(I-1) H0=H(I-1) TC=T0 HC=H0 CALL RUTTUN(TC,HC,KT0,KH0,FC,GC) TC=T0+0.5*KT0 HC=H0+0.5*KH0 CALL RUTTUN(TC,HC,KT1,KH1,FC,GC) TC=T0+0.5*KT1 HC=H0+0.5*KH1 CALL RUTTUN(TC,HC,KT2,KH2,FC,GC) TC=T0+KT2 HC=H0+KH2 CALL RUTTUN(TC,HC,KT3,KH3,FC,GC) T(I)=T(I-1)+(KT0+2.0*KT1+2.0*KT2+KT3)/6.0 H(I)=H(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0 1001 CONTINUE FPIN(INTR)=T(N) IF (INTR.EQ.1) THEN INTR=2 HGUESS(INTR)=HGUESS(1)-0.0001 GO TO 101 ELSE IF(INTR.GT.200) GO TO 201 IF(ABS(1.0-FPIN(INTR)).LT.0.00001) GO TO 301 INTR=INTR+1 DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)- & HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1)) HGUESS(INTR)=HGUESS(INTR-1)+DGUESS GO TO 101 ENDIF 201 WRITE(6,2000) GO TO 7778 301 WRITE(6,3001) INTR 7778 CONTINUE WRITE(1,4500) DO 7000 I=1,N ETA=(I-1)*DETA WRITE(1,4000) ETA,F(I),G(I),H(I),T(I) WRITE(2,5500) ETA,G(I),T(I) 7000 CONTINUE TGRAD=(T(2)-T(1))/DETA WRITE(1,4280) TGRAD WRITE(6,4280) TGRAD WRITE(1,4281) HU WRITE(6,4281) HU STOP CLOSE(1) CLOSE(2) 2000 FORMAT(' ******** FAILURE TO CONVERGE *********') 3000 FORMAT(' VELOCITY FUNCTION CONVERGENCE IN ',I5,' ITERATIONS') 3001 FORMAT(' TEMPERATURE FUNCTION CONVERGENCE IN ',I5,' ITERATIONS') 4000 FORMAT(5F10.6) 4250 FORMAT(////) 4260 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 4265 FORMAT(' INPUT THE VALUE OF n THEN ') 4270 FORMAT(//,' PRANDTL NUMBER = ',F12.3,//) 4280 FORMAT(//,' THETA GRADIENT AT WALL = ',F9.3) 4281 FORMAT(//,' F PRIME GRADIENT AT WALL = ',F9.3) 4290 FORMAT(/,' SIMILARITY SOLUTION FOR FLOW OVER ISOTHERMAL PLATE',/, $ ' __________________________________________________',//) 4500 FORMAT(' ETA F dF/DETA D2F/DETA2 THETA') 5500 FORMAT(F10.6,',',F10.6,',',F10.6) END * *********************************************************************** * * RUTTA-KUNGE SOLUTION ROUTINE FOR VELOCITY * *********************************************************************** * SUBROUTINE RUTKUN(F,G,H,KF,KG,KH) REAL KF,KG,KH COMMON DETA,XN,PR KF=DETA*G KG=DETA*H KH=-0.5*F*H*DETA RETURN END * *********************************************************************** * * RUTTA-KUNGE SOLUTION ROUTINE FOR TEMPERATURE * *********************************************************************** * SUBROUTINE RUTTUN(T,H,KT,KH,F,G) REAL KT,KH COMMON DETA,XN,PR KT=DETA*H KH=-0.5*PR*F*H*DETA+XN*G*PR*(T-1.0)*DETA RETURN END ************************************************************************ ***** ***** ***** TURBINRK ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES SOLVES THE MOMENTUM INTEGRAL EQUATION ***** ***** ***** ***** TOGETHER WITH AUXILARY EQUATIONS FOR THE SHAPE FACTOR AND ***** ***** ***** ***** WALL SHEAR STRESS FOR TURBULENT BOUNDARY LAYER FLOW OVER ***** ***** ***** ***** AN ISOTHERMAL SURFACE THE HEAT TRANSFER RATE BEING OBTAINED ***** ***** ***** ***** USING THE ANALOGY SOLUTION. THE RUNGE-KUTTA SOLUTION ***** ***** ***** ***** PROCEDURE IS ADOPTED. ***** ***** ***** ************************************************************************ * real NUX,NUL common DD2X,DHX,RE,H,D2,U,UUDX,CF2 * *********************************************************************** * open(UNIT=1,FILE='TURBINOU.DAT') * *********************************************************************** * write(6,1000) write(6,1100) write(6,1200) read(5,*) Re write(6,1300) read(5,*) A,B,C,D write(6,1400) read(5,*) PR write(6,1500) write(1,2000) write(1,2100) Re write(1,2200) A,B,C,D write(1,2300) Pr write(1,3500) * *********************************************************************** * * Re=REYNOLDS NUMBER * A,B,C,D=VELOCITY DISTRIBUTION COEFFICINTS * Pr=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * DX=0.002 X=DX U=A+B*X+C*X*X+D*X*X*X H=1.4 D2=( 0.0412*(X**0.7886))/((RE*U)**0.2114) CF2=(0.123*(10**(-0.678*H)))/((RE*U*D2)**0.268) * *********************************************************************** * 100 X=X+DX write(6,1600) X H0=H D20=D2 U=A+B*X+C*X*X+D*X*X*X UUDX=(B+2*C*X+3*D*X*X)/U CF2=(0.123*(10**(-0.678*H)))/((RE*U*D2)**0.268) if(CF2.lt.0.0001) then write(6,4000) write(1,4000) goto 200 endif call RUTKUN D2=D20+0.5*DD2X*DX H=H0+0.5*DHX*DX DD2X0=DD2X DHX0=DHX call RUTKUN D2=D20+0.5*DD2X*DX H=H0+0.5*DHX*DX DD2X1=DD2X DHX1=DHX call RUTKUN D2=D20+DD2X*DX H=H0+DHX*DX DD2X2=DD2X DHX2=DHX call RUTKUN DD2X3=DD2X DHX3=DHX H=H0+(DHX0+2.0*DHX1+2.0*DHX2+DHX3)*DX/6.0 D2=D20+(DD2X0+2.0*DD2X1+2.0*DD2X2+DD2X3)*DX/6.0 NUL=CF2*RE*U NUX=NUL*X REX=RE*X*U write(1,3000) X,REX,CF2,NUX IF(X.LT.1.0) GOTO 100 200 close(1) stop 1000 FORMAT(/,' TURBULENT BOUNDARY LAYER SOLUTION FOR FLOW OVER ',/, &' AN ISOTHERMAL SURFACE',//) 1100 FORMAT(////) 1200 FORMAT(' INPUT THE VALUE OF THE REYNOLDS NUMBER THEN ') 1300 FORMAT(' INPUT THE VALUES OF VEL. COEFFS. A,B,C,D THEN ') 1400 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 1500 FORMAT(////) 1600 FORMAT(' X = ',F12.6) 2000 FORMAT(/,' TURBULENT BOUNDARY LAYER SOLUTION FOR FLOW OVER ',/, &' AN ISOTHERMAL SURFACE',//) 2100 FORMAT(//,' REYNOLDS NUMBER= ',F12.1,/) 2200 FORMAT(//,' A= ',F9.4,' B= ',F9.4,' C= ',F9.4,' D= ',F9.4,/) 2300 FORMAT(//,' PRANDTL NUMBER= ',F12.3,//) 3000 FORMAT(F8.5,F14.2,F10.6,F12.4) 3500 FORMAT(' X Rex Cf/2 Nux',/) 4000 FORMAT(//,' ****** SEPARATION IMMINENT ******') end * *********************************************************************** * * RUTTA-KUNGE SOLUTION ROUTINE * *********************************************************************** * SUBROUTINE RUTKUN common DD2X,DHX,RE,H,D2,U,UUDX,CF2 DD2X=CF2-(H+2)*D2*UUDX UDRE=(U*D2*RE)**(1/6) DHX=exp(5*(H-1.4))*(-UDRE*D2*UUDX-0.0135*(H-1.4))/(UDRE*D2) return end ************************************************************************ ***** ***** ***** TURBOUND ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR TURBULENT BOUNDARY LAYER FLOW ***** ***** ***** ***** USING A FINITE DIFFERENCE METHOD. THE BOUNDARY LAYER IS ***** ***** ***** ***** LAMINAR UP TO THE TRANSITION REYNOLDS NUMBER. ***** ***** ***** ************************************************************************ * dimension U(2,250),V(2,250),T(2,250),A(250),B(250),C(250),D(250), * H(250),Y(250),E(250) common AH,BH,CH,DH,AV,BV,CV,DV,XMAX real NX,L(250),NV C C ******************************************************************** C open (1,file='TURBONPR.DAT') open (2,file='TURBONPL.DAT') C C ************************* INPUT ************************************* C write(6,6000) write(6,1000) write(6,1001) write(6,6001) read(5,*) XMAX,XTRAN,PR write(6,6003) read(5,*) AV,BV,CV,DV write(6,6002) read(5,*) AH,BH,CH,DH write (1,1002) XMAX,XTRAN,PR,AH,BH,CH,DH,AV,BV,CV,DV * C ********************************************************************* * X = 0.0 PRT = 0.9 N = 150 * *************** DEFINE NODAL POINTS IN Y-DIRECTION ****************** * Y(N)=1.0*sqrt(XMAX) Y(1) = 0.0 N1 = N-1 N2 = N-2 EXPD=1.05 DY1=((EXPD-1.0)/(EXPD**N-1))*(Y(N)-Y(1)) DY=DY1 do 1 J = 2,N Y(J) = Y(J-1)+DY 1 DY=EXPD*DY * *********************** ASSIGN INITIAL VALUES *********************** * DU = Y(3) call VELDP(X,U(2,N),DPX) call TEMP (X,T(1,1)) U(1,1) = 0.0 V(1,1) = 0.0 * do 3 J = 2, N U(1,J) = 1.0 V(1,J) = 0.0 3 T(1,J) = 0.0 DX = XMAX/050000.0 DXMIN=DX DXMAX = 10.0*DX * *************************** SOLUTION BEGINS ************************** * 100 DX = 1.1*DX if(DX.GT.DXMAX) DX=DXMAX X = X+DX * ************************** FIND EDDY VISCOSITY ************************ * YSC = SQRT(U(2,2)/Y(2)-Y(2)*DPX/2) do 10 J = 2, N1 if(X.LE.XTRAN) then E(J)=0.0 go to 10 endif YD = Y(J)/DU if(YD.GT.0.6) GO TO 11 if(YD.LT.0.1) GO TO 12 YS = Y(J)*YSC L(J) = 0.41*(1.0-EXP(-YS/26.0))*YD-1.54*(YD-0.1)**2 & +2.76*(YD-0.1)**3-1.88*(YD-0.1)**4 L(J) = L(J)*DU go to 13 11 L(J) = 0.089*DU go to 13 12 YS = Y(J)*YSC L(J) = 0.41*(1.0-EXP(-YS/26.0))*Y(J) 13 continue DYP=Y(J+1)-Y(J) DYM=Y(J)-Y(J-1) DY = (Y(J+1)-Y(J-1)) E(J) = (L(J)**2)*ABS((DYM/(DYP*DY))*(U(1,J+1)-U(1,J))+ & (DYP/(DYM*DY))*(U(1,J)-U(1,J-1))) 10 continue E(1)=0.0 E(N)=0.0 * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ******************* * U(2,1) = 0.0 call VELDP (X,U(2,N),DPX) A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=U(2,1) do 20 J = 2, N1 DYP = Y(J+1)-Y(J) DYM = Y(J)-Y(J-1) DY = (Y(J+1)-Y(J-1)) A(J) = (U(1,J)/DX)+V(1,J)*(DYP/(DYM*DY)-DYM/(DYP*DY))+ & ((2.0+E(J+1)+E(J))/DYP+(2.0+E(J)+E(J-1))/DYM)/DY B(J) = V(1,J)*(DYM/(DYP*DY))-((2.0+E(J+1)+E(J))/DYP)/DY C(J) = -V(2,J)*(DYP/(DYM*DY))-((2.0+E(J)+E(J-1))/DYM)/DY 20 D(J) = U(1,J)**2/DX-DPX A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=U(2,N) call TRISOL (A,B,C,D,H,N) do 21 J = 1, N 21 U(2,J) = H(J) * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * call TEMP (X,T(2,1)) T(2,N) = 0.000001 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,1) do 22 J = 2, N1 DYP = Y(J+1)-Y(J) DYM = Y(J)-Y(J-1) DY = Y(J+1)-Y(J-1) A(J) = (U(1,J)/DX)+V(1,J)*(DYP/(DYM*DY)-DYM/(DYP*DY))+ & ((2.0/PR+E(J+1)/PRT+E(J)/PRT)/DYP+ & (2.0/PR+E(J)/PRT+E(J-1)/PRT)/DYM)/DY B(J)=V(1,J)*(DYM/(DYP*DY))-((2.0/PR+E(J+1)/PRT+E(J)/PRT)/DYP)/DY C(J)=-V(2,J)*(DYP/(DYM*DY))-((2.0/PR+E(J)/PRT+E(J-1)/PRT)/DYM)/DY 22 D(J) = U(1,J)*T(1,J)/DX A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=T(2,N) call TRISOL (A,B,C,D,H,N) DO 23 J = 1, N 23 T(2,J) = H(J) * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** LAM01330 * V(2,1) = 0.0 do 24 J = 2,N 24 V(2,J) = V(2,J-1)-((Y(J)-Y(J-1))/(2.0*DX))*(U(2,J) & -U(1,J)+U(2,J-1)-U(1,J-1)) * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UU = 0.99*U(2,N) TT = 0.01*T(2,1) do 30 J = 1, N1 if(U(2,J).GE.UU) GO TO 31 30 continue 31 DU = Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) do 32 J =1, N1 if(T(2,J).le.TT) go to 33 32 continue 33 DT = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * if(DU.LE.Y(N-5).and.DT.LE.Y(N-5)) go to 34 NMIN = N+1 NMAX = N+10 if(NMAX.gt.250) go to 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) do 35 J = NMIN, NMAX Y(J) =Y(J-1)+DY U(2,J) = U(2,N) T(2,J) = T(2,N) 35 V(2,J) = V(2,J-1)+DVY*DY write (1,3000) N = NMAX N1=N-1 N2=N-2 write (1,2000)X write (1,4001) do 50 J = 1, N 50 write (1,4002) J,Y(J),U(2,J),V(2,J),T(2,J),E(J) 34 continue * ********************** TRANSFER VALUES ********************************* * do 40 J = 1, N U(1,J) = U(2,J) V(1,J) = V(2,J) T(1,J) = T(2,J) 40 continue * ******************** WRITE THE RESULTS ********************************* * DTY = (T(2,1)-T(2,2))/(Y(2)-Y(1)) NX = DTY*X/T(2,1) write (6,4000)X,DTY,DU,DT write (1,6500)X,DTY,NX write (2,6550)X,NX * ******************* CHECK FOR MAXIMUM X ******************************** * if(X.GT.XMAX) go to 102 go to 100 101 write (6,5000) go to 103 102 write (6,5001) 103 continue write (1,2000) X do 650 J = 1, N 650 write (1,4002) J,Y(J),U(2,J),V(2,J),T(2,J),E(J) close(1) close(2) stop * ************************** FORMAT STATEMENTS ************************** * 1000 format(10X,'TURBULENT BOUNDARY LAYER FLOW') 1001 format(10X,'-----------------------------') 1002 format(5X,' MAXIMUM REX =', F10.1, $ /,5X,' TRANSITION REX =',F10.1, $ /,5X,' PRANDTL NO. =',F10.4, $ /,5X, 'TEMPERATURE COEFFS. =',4F10.4, $ /,5X, 'VELOCITY COEFFS. =',4F10.4,//) 2000 format(5X,'SOLUTION FOR X = ', E12.4,/) 3000 format(5X,'NUMBER OF GRID POINTS INCREASED',/) 4000 format(5X,'X=',E12.4,' dT/dY=',E12.4,' DELU',E12.4,' DELT=',E12.4) 4001 format(3X,'J',7X,'Y',12X,'U',12X,'V',12X,'T',12X,'E') 4002 format(I5,5E12.4) 5000 format(15X,'N EXCEEDS 250') 5001 format(/,15X,'XMAX REACHED') 6000 format(////////) 6001 format(/,' INPUT THE VALUES OF THE MAXIMUM REX TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN, THE TRANSITION REX AND', & /,' OF THE PRANDTL NUMBER - (XMAX,XTRAN,PR)',//) 6002 format(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE WALL TEMPERATURE VARIATION - & (AH,BH,CH,DH) ',//) 6003 format(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE FREESTEAM VELOCITY VARIATION - & (AV,BV,CV,DV) ',//) 6500 format (' X = ',F14.4,' dT/dY = ',F14.7,' Nux = ',F14.4) 6550 format (F14.4,',',F14.4) end C C ********************************************************************* C SUBROUTINE TRISOL(A,B,C,D,H,N) C C *************** TRI-DIAGONAL MATRIX SOLVER ************************** C C *** THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM ******** C dimension A(250),B(250),C(250),D(250),H(250),W(250),R(250),G(250) W(1)=A(1) G(1)=D(1)/W(1) do 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 continue H(N)=G(N) N1=N-1 do 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 continue return end * ******************************************************************** * SUBROUTINE TEMP(X,TW) common AH,BH,CH,DH,AV,BV,CV,DV,XMAX * *********** THIS DETERMINES THE WALL TEMPERATURE ******************* * XR=X/XMAX TW = AH+BH*XR+CH*XR*XR+DH*XR*XR*XR return end * ************************************************************************ * SUBROUTINE VELDP (X,UF,DPX) common AH,BH,CH,DH,AV,BV,CV,DV,XMAX * ****** THIS DETERMINES FREESTREAM VELOCITY AND PRESSURE GRADIENT ****** * XR=X/XMAX UF = AV+BV*XR+CV*XR*XR+DV*XR*XR*XR DPX = -UF*(BV+2.0*CV*X+3.0*DV*X*X)/XMAX return end ************************************************************************ ***** ***** ***** TURBOUNQ ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR TURBULENT BOUNDARY LAYER FLOW ***** ***** ***** ***** USING A FINITE DIFFERENCE METHOD. THE BOUNDARY LAYER IS ***** ***** ***** ***** LAMINAR UP TO THE TRANSITION REYNOLDS NUMBER. THERE IS A ***** ***** ***** ***** SPECIFIED HEAT FLUX AT THE SURFACE. ***** ***** ***** ************************************************************************ * dimension U(2,250),V(2,250),T(2,250),A(250),B(250),C(250),D(250), * H(250),Y(250),E(250) common AQ,BQ,CQ,DQ,AV,BV,CV,DV,XMAX real NX,L(250),NV C C ******************************************************************** C open (1,file='TURBOQPR.DAT') open (2,file='TURBOQPL.DAT') C C ************************* INPUT ************************************* C write(6,6000) write(6,1000) write(6,1001) write(6,6001) read(5,*) XMAX,XTRAN,PR write(6,6003) read(5,*) AV,BV,CV,DV write(6,6002) read(5,*) AQ,BQ,CQ,DQ write (1,1002) XMAX,XTRAN,PR,AQ,BQ,CQ,DQ,AV,BV,CV,DV * C ********************************************************************* * X = 0.0 PRT = 0.9 N = 150 * *************** DEFINE NODAL POINTS IN Y-DIRECTION ****************** * Y(N)=1.0*sqrt(XMAX) Y(1) = 0.0 N1 = N-1 N2 = N-2 EXPD=1.05 DY1=((EXPD-1.0)/(EXPD**N-1))*(Y(N)-Y(1)) DY=DY1 do 1 J = 2,N Y(J) = Y(J-1)+DY 1 DY=EXPD*DY * *********************** ASSIGN INITIAL VALUES *********************** * DU = Y(3) call VELDP(X,U(2,N),DPX) call FLUX (X,QW) U(1,1) = 0.0 T(1,1) = Y(2)*QW V(1,1) = 0.0 * do 3 J = 2, N U(1,J) = 1.0 V(1,J) = 0.0 3 T(1,J) = 0.0 DX = XMAX/050000.0 DXMIN=DX DXMAX = 10.0*DX * *************************** SOLUTION BEGINS ************************** * 100 DX = 1.1*DX if(DX.GT.DXMAX) DX=DXMAX X = X+DX * ************************** FIND EDDY VISCOSITY ************************ * YSC = SQRT(U(2,2)/Y(2)-Y(2)*DPX/2) do 10 J = 2, N1 if(X.LE.XTRAN) then E(J)=0.0 go to 10 endif YD = Y(J)/DU if(YD.GT.0.6) GO TO 11 if(YD.LT.0.1) GO TO 12 YS = Y(J)*YSC L(J) = 0.41*(1.0-EXP(-YS/26.0))*YD-1.54*(YD-0.1)**2 & +2.76*(YD-0.1)**3-1.88*(YD-0.1)**4 L(J) = L(J)*DU go to 13 11 L(J) = 0.089*DU go to 13 12 YS = Y(J)*YSC L(J) = 0.41*(1.0-EXP(-YS/26.0))*Y(J) 13 continue DYP=Y(J+1)-Y(J) DYM=Y(J)-Y(J-1) DY = (Y(J+1)-Y(J-1)) E(J) = (L(J)**2)*ABS((DYM/(DYP*DY))*(U(1,J+1)-U(1,J))+ & (DYP/(DYM*DY))*(U(1,J)-U(1,J-1))) 10 continue E(1)=0.0 E(N)=0.0 * ***************** SOLVE MOMENTUM EQUATION TO GET "U" ******************* * U(2,1) = 0.0 call VELDP(X,U(2,N),DPX) A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=U(2,1) do 20 J = 2, N1 DYP = Y(J+1)-Y(J) DYM = Y(J)-Y(J-1) DY = (Y(J+1)-Y(J-1)) A(J) = (U(1,J)/DX)+V(1,J)*(DYP/(DYM*DY)-DYM/(DYP*DY))+ & ((2.0+E(J+1)+E(J))/DYP+(2.0+E(J)+E(J-1))/DYM)/DY B(J) = V(1,J)*(DYM/(DYP*DY))-((2.0+E(J+1)+E(J))/DYP)/DY C(J) = -V(2,J)*(DYP/(DYM*DY))-((2.0+E(J)+E(J-1))/DYM)/DY 20 D(J) = U(1,J)**2/DX-DPX A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=U(2,N) call TRISOL (A,B,C,D,H,N) do 21 J = 1, N 21 U(2,J) = H(J) * ****************** SOLVE ENERGY EQUATION TO GET "T" ******************** * call FLUX (X,QW) T(2,N) = 0.000001 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=QW*Y(2) do 22 J = 2, N1 DYP = Y(J+1)-Y(J) DYM = Y(J)-Y(J-1) DY = Y(J+1)-Y(J-1) A(J) = (U(1,J)/DX)+V(1,J)*(DYP/(DYM*DY)-DYM/(DYP*DY))+ & ((2.0/PR+E(J+1)/PRT+E(J)/PRT)/DYP+ & (2.0/PR+E(J)/PRT+E(J-1)/PRT)/DYM)/DY B(J)=V(1,J)*(DYM/(DYP*DY))-((2.0/PR+E(J+1)/PRT+E(J)/PRT)/DYP)/DY C(J)=-V(2,J)*(DYP/(DYM*DY))-((2.0/PR+E(J)/PRT+E(J-1)/PRT)/DYM)/DY 22 D(J) = U(1,J)*T(1,J)/DX A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=T(2,N) call TRISOL (A,B,C,D,H,N) DO 23 J = 1, N 23 T(2,J) = H(J) * ************* SOLVE THE CONTINUITY EQUATION TO GET "V" *************** LAM01330 * V(2,1) = 0.0 do 24 J = 2,N 24 V(2,J) = V(2,J-1)-((Y(J)-Y(J-1))/(2.0*DX))*(U(2,J) & -U(1,J)+U(2,J-1)-U(1,J-1)) * ************* FIND THE BOUNDARY LAYER THICKNESSES ******************** * UU = 0.99*U(2,N) TT = 0.01*T(2,1) do 30 J = 1, N1 if(U(2,J).GE.UU) GO TO 31 30 continue 31 DU = Y(J)-(Y(J)-Y(J-1))*(U(2,J)-UU)/(U(2,J)-U(2,J-1)) do 32 J =1, N1 if(T(2,J).le.TT) go to 33 32 continue 33 DT = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J)) * ************* CHECK POSITION OF OUTER GRID POINTS ******************** * if(DU.LE.Y(N-5).and.DT.LE.Y(N-5)) go to 34 NMIN = N+1 NMAX = N+10 if(NMAX.gt.250) go to 101 DVY = (V(2,N)-V(2,N-1))/(Y(N)-Y(N-1)) do 35 J = NMIN, NMAX Y(J) =Y(J-1)+DY U(2,J) = U(2,N) T(2,J) = T(2,N) 35 V(2,J) = V(2,J-1)+DVY*DY write (1,3000) N = NMAX N1=N-1 N2=N-2 write (1,2000)X write (1,4001) do 50 J = 1, N 50 write (1,4002) J,Y(J),U(2,J),V(2,J),T(2,J),E(J) 34 continue * ********************** TRANSFER VALUES ********************************* * do 40 J = 1, N U(1,J) = U(2,J) V(1,J) = V(2,J) T(1,J) = T(2,J) 40 continue * ******************** WRITE THE RESULTS ********************************* * write (6,4000)X,T(2,1),DU,DT NX=X/T(2,1) write (1,6500)X,T(2,1),NX write (2,6550)X,T(2,1),NX * ******************* CHECK FOR MAXIMUM X ******************************** * if(X.GT.XMAX) go to 102 go to 100 101 write (6,5000) go to 103 102 write (6,5001) 103 continue write (1,2000) X do 650 J = 1, N 650 write (1,4002) J,Y(J),U(2,J),V(2,J),T(2,J),E(J) close(1) close(2) stop * ************************** FORMAT STATEMENTS ************************** * 1000 format(10X,'TURBULENT BOUNDARY LAYER FLOW') 1001 format(10X,'-----------------------------') 1002 format(5X,' MAXIMUM REX =', F10.1, $ /,5X,' TRANSITION REX =',F10.1, $ /,5X,' PRANDTL NO. =',F10.4, $ /,5X, 'HEAT FLUX COEFFS. =',4F10.4, $ /,5X, 'VELOCITY COEFFS. =',4F10.4,//) 2000 format(5X,'SOLUTION FOR X = ', E12.4,/) 3000 format(5X,'NUMBER OF GRID POINTS INCREASED',/) 4000 format(5X,'X=',E12.4,' T w=',E12.4,' DELU',E12.4,' DELT=',E12.4) 4001 format(3X,'J',7X,'Y',12X,'U',12X,'V',12X,'T',12X,'E') 4002 format(I5,5E12.4) 5000 format(15X,'N EXCEEDS 250') 5001 format(/,15X,'XMAX REACHED') 6000 format(////////) 6001 format(/,' INPUT THE VALUES OF THE MAXIMUM REX TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN, THE TRANSITION REX AND', & /,' OF THE PRANDTL NUMBER - (XMAX,XTRAN,PR)',//) 6002 format(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE WALL HEAT FLUX VARIATION - & (AQ,BQ,CQ,DQ) ',//) 6003 format(/,' INPUT THE VALUES OF THE COEFFICIENTS OF THE 3RD ORDER', & /,' POLYNOMIAL THAT DESCRIBES THE FREESTEAM VELOCITY VARIATION - & (AV,BV,CV,DV) ',//) 6500 format (' X = ',F14.4,' T w = ',F14.7,' Nux = ',F14.4) 6550 format (F14.4,',',F14.4,',',F14.4) end C C ********************************************************************* C SUBROUTINE TRISOL(A,B,C,D,H,N) C C *************** TRI-DIAGONAL MATRIX SOLVER ************************** C C *** THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM ******** C dimension A(250),B(250),C(250),D(250),H(250),W(250),R(250),G(250) W(1)=A(1) G(1)=D(1)/W(1) do 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 continue H(N)=G(N) N1=N-1 do 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 continue return end * ******************************************************************** * SUBROUTINE FLUX(X,QW) common AQ,BQ,CQ,DQ,AV,BV,CV,DV,XMAX * *********** THIS DETERMINES THE WALL TEMPERATURE ******************* * XR=X/XMAX QW = AQ+BQ*XR+CQ*XR*XR+DQ*XR*XR*XR return end * ************************************************************************ * SUBROUTINE VELDP (X,UF,DPX) common AQ,BQ,CQ,DQ,AV,BV,CV,DV,XMAX * ****** THIS DETERMINES FREESTREAM VELOCITY AND PRESSURE GRADIENT ****** * XR=X/XMAX UF = AV+BV*XR+CV*XR*XR+DV*XR*XR*XR DPX = -UF*(BV+2.0*CV*XR+3.0*DV*XR*XR)/XMAX return end ************************************************************************ ***** ***** ***** TURDUCD ***** ***** _______ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM SOLVES FOR TURBULENT THERMALLY ***** ***** ***** ***** DEVELOPING FLOW IN A CIRCULAR PIPE WITH A UNIFORM ***** ***** ***** ***** WALL TEMPERATURE. ***** ***** ***** ************************************************************************ * dimension U(250),T(2,250),A(250),B(250),C(250),D(250), * H(250),R(250),E(250),ER(250),PRD(250) real ND,NDA C C ******************************************************************** C open (1,file='TURDUCPR.DAT') open (2,file='TURDUCPL.DAT') C C ************************* INPUT ************************************* C write(6,6000) write(6,1000) write(6,1001) write(6,6001) read(5,*) ZMAX,PR,RE write (1,1002) ZMAX,PR,RE write(6,6002) read(5,*) ITQ if(ITQ.eq.1) then write (1,1003) else write (1,1004) endif * C ********************************************************************* * Z = 0.0 PRT = 0.9 N = 100 * *************** DEFINE NODAL POINTS IN Y-DIRECTION ****************** * R(N)=0.5 R(1) = 0.0 N1 = N-1 N2 = N-2 EXPD=1.05 DR1=((EXPD-1.0)/(EXPD**(N-1)-1.0))*(R(N)-R(1)) DR=DR1 do 1 K = 2,N J=N+1-K R(J) = R(J+1)-DR 1 DR=EXPD*DR * *********************** ASSIGN INITIAL VALUES *********************** * do 3 J = 1, N RR=1.0-2.0*R(J) U(J) = (60.0/49.0)*(RR**(1.0/7.0)) 3 T(1,J) = 0.0 T(1,N)=1.0 DZ = ZMAX/050000.0 DZMIN=DZ DZMAX = 10.0*DZ * ************************** FIND EDDY VISCOSITY ************************ * CF2 = 0.03955/RE**0.25 CFS=sqrt(CF2) PRT=0.9 PRRT=PR/PRT do 10 J = 1, N YD=0.5-R(J) YPLUS=YD*RE*CFS if(YPLUS.LT.5.0) E(J)=0.0 if(YPLUS.GE.5.0.and.YPLUS.LE.30.0) E(J)=YPLUS/5.0-1.0 if(YPLUS.GT.30.0) E(J)=0.8*(1.0-YD/0.5)*YPLUS**(6.0/7.0)-1.0 ER(J) = R(J)*(1.0+PRRT*E(J)) 10 continue * *************************** SOLUTION BEGINS ************************** * 100 DZ = 1.05*DZ if(DZ.GT.DZMAX) DZ=DZMAX Z = Z+DZ * ***************** SOLVE ENERGY EQUATION TO GET "T" ******************* * A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 do 20 J = 2, N1 DRP = R(J+1)-R(J) DRM = R(J)-R(J-1) DR = R(J+1)-R(J-1) A(J) = U(J)/DZ + & ((ER(J+1)+ER(J))/DRP+(ER(J)+ER(J-1))/DRM)/(R(J)*DR) B(J) = - ((ER(J+1)+ER(J))/DRP)/(R(J)*DR) C(J) = -((ER(J)+ER(J-1))/DRM)/(R(J)*DR) 20 D(J) = U(J)*T(1,J)/DZ if(ITQ.eq.1) then A(N)=1.0 B(N)=0.0 C(N)=0.0 D(N)=1.0 else A(N)=1.0 B(N)=0.0 C(N)=-1.0 D(N)=R(N)-R(N1) endif call TRISOL (A,B,C,D,H,N) do 21 J = 1, N 21 T(2,J) = H(J) * ********************** TRANSFER VALUES ********************************* * do 40 J = 1, N T(1,J) = T(2,J) 40 continue * ******************** WRITE THE RESULTS ********************************* * DTR = (T(2,N)-T(2,N1))/(R(N)-R(N1)) do 45 J = 1, N PRD(J)=U(J)*T(2,J)*R(J) 45 continue SUT=0.0 do 46 J = 2, N SUT=SUT+(PRD(J)+PRD(J-1))*(R(J)-R(J-1))/2.0 46 continue TAVG=8.0*SUT if(ITQ.eq.1) then ND = DTR NDA=ND/(1.0-TAVG) else ND = 1.0/T(2,N) NDA=1.0/(T(2,N)-TAVG) endif if(ITQ.eq.1) then write (6,4000)Z,DTR write (1,6500)Z,DTR,ND,NDA write (2,6550)Z,ND,NDA else write (6,4010)Z,T(2,N) write (1,6510)Z,T(2,N),ND,NDA write (2,6560)Z,T(2,N),ND,NDA endif * ******************* CHECK FOR MAXIMUM Z ******************************** * if(Z.GT.ZMAX) go to 102 go to 100 102 write (6,5001) 103 continue write (1,2000) Z write (1,4001) do 650 J = 1, N 650 write (1,4002) J,R(J),U(J),T(2,J),E(J) close(1) close(2) stop * ************************** FORMAT STATEMENTS ************************** * 1000 format(10X,' THERMALLY DEVELOPING TURBULENT PIPE FLOW') 1001 format(10X,' ----------------------------------------') 1002 format(5X,' MAXIMUM Z =', F10.4, $ /,5X,' PRANDTL NO. =',F10.4, $ /,5X,' REYNOLDS NO. =',F10.1,//) 1003 format(5X,' WALL TEMPERATURE UNIFORM',//) 1004 format(5X,' WALL HEAT FLUX UNIFORM',//) 2000 format(5X,'SOLUTION FOR Z = ', F12.8,/) 4000 format(5X,'Z=',F10.8,' dT/dR=',F14.7) 4001 format(3X,'J',7X,'R',12X,'U',12X,'T',12X,'E') 4002 format(I5,4E12.4) 4010 format(5X,'Z=',F10.8,' T WALL=',F14.7) 5001 format(/,15X,'ZMAX REACHED') 6000 format(////////) 6001 format(/,' INPUT THE VALUES OF THE MAXIMUM Z TO WHICH THE', & /,' SOLUTION IS TO BE UNDERTAKEN, OF THE PRANDTL NUMBER', & /,' AND OF THE REYNOLDS NUMBER - (ZMAX,PR,RE)',//) 6002 format(/,' IS WALL (A) TEMPERATURE OR', & /,' (B) HEAT FLUX UNIFORM - (A)=1,(B)=2',//) 6500 format (' Z = ',F10.8,' dT/dR = ',F10.4,' NuD = ',F10.4, &' NuDa = ',F10.4) 6510 format (' Z = ',F10.8,' T WALL = ',F10.4,' NuDi = ',F10.4, &' NuDa = ',F10.4) 6550 format (F10.8,',',F10.4,',',F10.4) 6560 format (F10.8,',',F10.4,',',F10.4,',',F10.4) end C C ********************************************************************* C SUBROUTINE TRISOL(A,B,C,D,H,N) C C *************** TRI-DIAGONAL MATRIX SOLVER ************************** C C *** THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM ******** C dimension A(250),B(250),C(250),D(250),H(250),W(250),R(250),G(250) W(1)=A(1) G(1)=D(1)/W(1) do 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 continue H(N)=G(N) N1=N-1 do 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 continue return end ************************************************************************ ***** ***** ***** TURINDEV ***** ***** ________ ***** ***** ***** ************************************************************************ ***** ***** ***** THIS PROGRAM GIVES SOLVES THE MOMENTUM INTEGRAL EQUATION ***** ***** ***** ***** TOGETHER WITH AUXILARY EQUATIONS FOR THE SHAPE FACTOR AND ***** ***** ***** ***** WALL SHEAR STRESS FOR TURBULENT DEVELOPING DUCT FLOW ***** ***** ***** ***** THE HEAT TRANSFER RATE BEING GIVEN BY THE ANALOGY SOLUTION. ***** ***** ***** ***** ************************************************************************ * REAL NUW,NUWA,NP * *********************************************************************** * OPEN(UNIT=1,FILE='TURINDOU.DAT') * *********************************************************************** * WRITE(6,1000) WRITE(6,1100) WRITE(6,1200) READ(5,*) Re WRITE(6,1400) READ(5,*) PR WRITE(6,1500) WRITE(1,2000) WRITE(1,2100) Re WRITE(1,2300) Pr WRITE(1,3500) * *********************************************************************** * * Re=REYNOLDS NUMBER * Pr=PRANDTL NUMBER * DETA=ETA STEP SIZE * *********************************************************************** * U=1.0 UUDX=0.0 DX=0.01 X=DX H=1.4 D2=( 0.0412*(X**0.7886))/((RE*U)**0.2114) CF2=(0.123*(10**(-0.678*H)))/((RE*U*D2)**0.268) DD2X=CF2 DHX=0.0 UUDX=(2.0/(1.0-2.0*H*D2))*(H*DD2X+D2*DHX) * *********************************************************************** * 100 X=X+DX WRITE(6,1600) X U=U+UUDX*DX DD2X=CF2-(H+2)*D2*UUDX UDRE=(U*D2*RE)**(1/6) DHX=exp(5*(H-1.4))*(-UDRE*D2*UUDX-0.0135*(H-1.4))/(UDRE*D2) D2=D2+DD2X*DX H=H+DHX*DX CF2=(0.123*(10**(-0.678*H)))/((RE*U*D2)**0.268) UUDX=(2.0/(1.0-2.0*H*D2))*(H*DD2X+D2*DHX) DB=(H*(H+1.0)/(H-1.0))*D2 NUW=CF2*RE*U*(PR**0.4) NP=2.0/(H-1.0) TRAT=1.0-((2.0*NP)/((NP+1.0)*(NP+2.0)))*U*DB NUWA=NUW/TRAT WRITE(1,3000) X,NUW,DB,NUWA,U IF(DB.LT.0.5) GOTO 100 STOP CLOSE(1) 1000 FORMAT(/,' DEVELOPING TURBULENT FLOW IN A DUCT',//) 1100 FORMAT(////) 1200 FORMAT(' INPUT THE VALUE OF THE REYNOLDS NUMBER THEN ') 1400 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ') 1500 FORMAT(////) 1600 FORMAT(' X = ',F12.6) 2000 FORMAT(/,' DEVELOPING TURBULENT FLOW IN A DUCT',//) 2100 FORMAT(' REYNOLDS NUMBER= ',F12.1,/) 2300 FORMAT(' PRANDTL NUMBER= ',F12.3,//) 3000 FORMAT(F8.5,F12.4,F10.5,F12.4,F10.5) 3500 FORMAT(' X Nuw B.L.Thick. Nuw mean U',/) END C C ********************************************************************* C C PROGRAM "MIXSQCYL" C C THIS PROGRAM SOLVES THE MIXED CONVECTIVE FLOW OVER A SQUARE CYLINDER. C STEADY SYMMETRICAL FLOW IS ASSUMED. C THE FULL NAVIER-STOKES AND ENERGY EQUATIONS C ARE USED IN OBTAINING THE SOLUTION. THE PROGRAM IS INTENDED TO C BE USED TO STUDY LOW REYNOLDS NUMBER FLOWS. C C ********************************************************************* C real LU,LD,LN,S(2,91,65),V(2,91,65),T(2,91,65),A(91),B(91),C(91), $ D(91),Q(91),H(91),QL(91) C C ******************************************************************** C open (1,file='MIXSQCPR.DAT') open (2,file='MIXSQCPL.DAT') C C ********************************************************************* C INPUT C ********************************************************************* C C RE=REYNOLDS NUMBER C PR=PRANDTL NUMBER C GR=GRASHOF NUMBER C IDIR=BUOYANCY FORCE DIRECTION PARAMETER ( +1 ASSISTING C FLOW, -1 OPPOSING FLOW) C LU=DIMENSIONLESS UPSTREAM DISTANCE C LD=DIMENSIONLESS DOWNSTREAM DISTANCE C LN=DIMENSIONLESS CROSS-STREAM DISTANCE C C ********************************************************************* C write (1,1002) write (1,1003) C C ********************************************************************* C write(6,1400) write(6,1040) read(5,*) RE if (RE.gt.100.0) then write(6,6996) 6996 format(' ******** RE100 *********') stop endif write(6,1050) read(5,*) PR write(6,1102) read(5,*) GR write(6,1500) read(5,*) IDIR LU=1.0 LD=2.0 LN=1.5 GB=GR/(RE*RE) C C ********************************************************************* C write(1,5237) RE,PR,GR,IDIR 5237 format(//,' REYNOLDS NUMBER= ', F10.1,/,' PRANDTL NUMBER= ',F12.3, $/,' GRASHOF NUMBER= ', F12.1,/,' DIRECTION PARAMETER= ',I4,//) C ********************************************************************* C REX=0.2 RES=0.7 DX=0.050 NT=15 DY=0.5/(NT-1) C C ********************************************************************* C DX2=DX*DX DY2=DY*DY DX2I=1.0/DX2 DY2I=1.0/DY2 RX2=1.0/(RE*DX2) RY2=1.0/(RE*DY2) PX2=RX2/PR PY2=RY2/PR DXI=1.0/(2.0*DX) DYI=1.0/(2.0*DY) NU=INT(LU/DX)+1 NP=INT(1.0/DX)+NU ND=INT(LD/DX)+NP NN=INT(LN/DY)+1 NN1=NN-1 ND1=ND-1 NU1=NU-1 NP1=NP-1 NT1=NT-1 NUP=NU+1 NTP=NT+1 NPP=NP+1 NTN=NN-NT+1 NDN=ND-NP+1 IW=NP+(ND-NP)/2 ITM=NU+NP/2 C C ********************************************************************* write(2,7237) RE,GR,PR,NU,NP,ND,NT,NN,DX,DY 7237 format(F9.2,','F10.1,',',F9.4,',',5(I4,','),F10.6,',',F10.6) C C ********************************************************************* C INITIAL VALUES C ********************************************************************* C do 100 I=1,ND do 100 J=1,NN S(1,I,J)=0.0 if(J.eq.NN) S(1,I,J)=LN if(I.eq.1) S(1,I,J)=(J-1)*LN/(NN-1) if(I.eq.ND) S(1,I,J)=(J-1)*LN/(NN-1) if(I.eq.NU.and.J.LE.NT) S(1,I,J)=0.0 if(I.eq.NP.and.J.LE.NT) S(1,I,J)=0.0 V(1,I,J)=0.0 T(1,I,J)=0.0 100 continue do 176 I=NU,NP S(1,I,NT)=0.0 176 continue do 200 I=NU,NP T(1,I,NT)=1.0 200 continue do 253 J=1,NT T(1,NU,J)=1.0 T(1,NP,J)=1.0 253 continue DO 6254 I=1,ND DO 6254 J=1,NN T(2,I,J)=T(1,I,J) 6254 CONTINUE C C ********************************************************************* C SOLVE INVISCID FLOW EQUATIONS TO GET INITIAL STREAM FUNCTION VALUES C ********************************************************************* C INTR=0 3100 continue C C ********************************************************************* C INTR=INTR+1 C C C ********************************************************************* C FIND STREAM FUNCTION ON Y-DIRECTION LINES C ********************************************************************* C do 5507 J=1,NN S(2,1,J)=S(1,1,J) 5507 continue do 5700 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 do 5600 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 5600 continue A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN call TRISOL(NN,A,B,C,D,H) do 5800 J=1,NN S(2,I,J)=S(1,I,J)+RES*(H(J)-S(1,I,J)) 5800 continue 5700 continue do 6791 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 do 6691 K=NTP,NN1 J=K-NT+1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,K)+S(1,I-1,K)) 6691 continue A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=LN call TRISOL(NTN,A,B,C,D,H) do 6891 K=1,NTN J=NT+K-1 S(2,I,J)=S(1,I,J)+RES*(H(K)-S(1,I,J)) 6891 continue 6791 continue do 6781 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 do 9600 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 9600 continue A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=S(1,I,NN) call TRISOL(NN,A,B,C,D,H) do 9800 J=1,NN S(2,I,J)=S(1,I,J)+RES*(H(J)-S(1,I,J)) 9800 continue 6781 continue do 5900 J=1,NN S(2,ND,J)=S(1,ND,J) 5900 continue C C ********************************************************************* C do 6207 I=1,ND do 6207 J=1,NN S(1,I,J)=S(2,I,J) 6207 continue C C ********************************************************************* C FIND STREAM FUNCTION ON X-DIRECTION LINES C ********************************************************************* C do 6510 I=1,ND S(2,I,1)=S(1,I,1) 6510 continue do 6110 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) do 7001 I=2,NU1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 7001 continue A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=S(1,NU,J) call TRISOL(NU,A,B,C,D,H) DIFFS=0.0 do 6211 I=1,NU DIFF=ABS(H(I)-S(1,I,J)) if(DIFF.gt.DIFFS) DIFFS=DIFF S(2,I,J)=S(1,I,J)+RES*(H(I)-S(1,I,J)) 6211 continue 6110 continue do 5910 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,NP,J) do 8182 K=NPP,ND1 I=K-NP+1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,K,J+1)+S(1,K,J-1)) 8182 continue A(NDN)=1.0 B(NDN)=0.0 C(NDN)=0.0 D(NDN)=S(1,ND,J) call TRISOL(NDN,A,B,C,D,H) DIFFX=0.0 do 6011 K=1,NDN I=NP+K-1 DIFF=ABS(H(K)-S(1,I,J)) if(DIFF.gt.DIFFS) DIFFS=DIFF S(2,I,J)=S(1,I,J)+RES*(H(K)-S(1,I,J)) 6011 continue 5910 continue do 8197 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) do 9091 I=2,ND1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 9091 continue A(ND)=1.0 B(ND)=0.0 C(ND)=0.0 D(ND)=S(1,ND,J) call TRISOL(ND,A,B,C,D,H) do 6097 I=1,ND DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFS) DIFFS=DIFF S(2,I,J)=S(1,I,J)+RES*(H(I)-S(1,I,J)) 6097 continue 8197 continue do 6250 I=1,ND do 6250 J=1,NN S(1,I,J)=S(2,I,J) 6250 continue C C ********************************************************************* C BOTTOM OF INVISCID FLOW STREAM FUNCTION LOOP C ********************************************************************* C write(6,5584) INTR 5584 format(' INVISCID FLOW LOOP, ITERATION NUMBER ',I6) C C ********************************************************************* C if(INTR.lt.300) GO TO 3100 if(DIFFS.lt.0.0001) GO TO 4107 if(INTR.gt.1500) GO TO 4107 go to 3100 4107 continue C C ********************************************************************* C INTR=0 C C ********************************************************************* C TOP OF VORTICITY - STREAM FUNCTION - TEMPERATURE LOOP C ********************************************************************* C 1000 continue C C ********************************************************************* C INTR=INTR+1 C C ********************************************************************* C CALCULATION OF WALL VORTICITY VALUES C ********************************************************************* C do 300 I=1,ND IF(I.GE.NU.AND.I.LE.NP) THEN V(2,I,NT)=V(1,I,NT)+REX*(-2.0*(S(1,I,NT+1)- $ S(1,I,NT))/DY2-V(1,I,NT)) ENDIF IF(I.LT.NU) V(2,I,1)=0.0 IF(I.GT.NP) V(2,I,1)=0.0 V(2,I,NN)=0.0 300 CONTINUE DO 353 J=1,NN V(2,1,J)=0.0 IF(J.LE.NT) THEN V(2,NU,J)=V(1,NU,J)+REX*(-2.0*(S(1,NU-1,J)- $ S(1,NU,1))/DY2-V(1,NU,J)) V(2,NP,J)=V(1,NP,J)+REX*(-2.0*(S(1,NP+1,J)- $ S(1,NP,1))/DY2-V(1,NP,J)) ENDIF 353 CONTINUE V(2,NU,NT)=V(1,NU,NT)+REX*(-2.0*(S(1,NU-1,NT)- $ S(1,NU,NT))/DX2-2.0*(S(1,NU,NT+1)- $ S(1,NU,NT))/DY2-V(1,NU,NT)) V(2,NP,NT)=V(1,NP,NT)+REX*(-2.0*(S(1,NP+1,NT)- $ S(1,NP,NT))/DX2-2.0*(S(1,NP,NT+1)- $ S(1,NP,NT))/DY2-V(1,NP,NT)) C C ********************************************************************* C FIND VORTICITY ON Y-DIRECTION LINES C ********************************************************************* C DO 600 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,1) DO 500 J=2,NN1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(V(1,I+1,J)-V(1,I-1,J)) $+RX2*(V(1,I+1,J)+V(1,I-1,J)) $-GB*IDIR*DYI*(T(1,I,J+1)-T(1,I,J-1)) 500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 700 J=1,NN V(2,I,J)=V(1,I,J)+REX*(H(J)-V(1,I,J)) 700 CONTINUE 600 CONTINUE DO 691 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,NT) DO 591 K=NTP,NN1 J=K-NT+1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-RY2 C(J)=+DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-RY2 D(J)=-DYI*(S(1,I,K+1)-S(1,I,K-1))*DXI*(V(1,I+1,K)-V(1,I-1,K)) $+RX2*(V(1,I+1,K)+V(1,I-1,K)) $-GB*IDIR*DYI*(T(1,I,K+1)-T(1,I,K-1)) 591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=0.0 CALL TRISOL(NTN,A,B,C,D,H) DO 791 K=1,NTN J=NT+K-1 V(2,I,J)=V(1,I,J)+REX*(H(K)-V(1,I,J)) 791 CONTINUE 691 CONTINUE DO 681 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,I,1) DO 581 J=2,NN1 A(J)=2.0*RX2+2.0*RY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-RY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(V(1,I+1,J)-V(1,I-1,J)) $+RX2*(V(1,I+1,J)+V(1,I-1,J)) $-GB*IDIR*DYI*(T(1,I,J+1)-T(1,I,J-1)) 581 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 781 J=1,NN V(2,I,J)=V(1,I,J)+REX*(H(J)-V(1,I,J)) 781 CONTINUE 681 CONTINUE DO 800 J=1,NN V(2,ND,J)=V(2,ND-1,J) 800 CONTINUE DO 2100 I=1,ND DO 2100 J=1,NN V(1,I,J)=V(2,I,J) 2100 CONTINUE C C ********************************************************************* C DO 2410 I=1,ND V(2,I,1)=V(1,I,1) 2410 CONTINUE C C ********************************************************************* C FIND VORTICITY ON X-DIRECTION LINES C ********************************************************************* C DO 1101 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 991 I=2,NU1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(V(1,I,J+1)-V(1,I,J-1)) $+RY2*(V(1,I,J+1)+V(1,I,J-1)) $-GB*IDIR*DYI*(T(1,I,J+1)-T(1,I,J-1)) 991 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=V(2,NU,J) CALL TRISOL(NU,A,B,C,D,H) DO 1191 I=1,NU V(2,I,J)=V(1,I,J)+REX*(H(I)-V(1,I,J)) 1191 CONTINUE 1101 CONTINUE DO 1192 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=V(2,NP,J) DO 1082 K=NPP,ND1 I=K-NP+1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-RX2 C(I)=-DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-RX2 D(I)=+DXI*(S(1,K+1,J)-S(1,K-1,J))*DYI*(V(1,K,J+1)-V(1,K,J-1)) $+RY2*(V(1,K,J+1)+V(1,K,J-1)) $-GB*IDIR*DYI*(T(1,K,J+1)-T(1,K,J-1)) 1082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DO 1374 K=1,NDN I=NP+K-1 V(2,I,J)=V(1,I,J)+REX*(H(K)-V(1,I,J)) 1374 CONTINUE 1192 CONTINUE DO 1010 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 900 I=2,ND1 A(I)=2.0*RX2+2.0*RY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-RX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(V(1,I,J+1)-V(1,I,J-1)) $+RY2*(V(1,I,J+1)+V(1,I,J-1)) $-GB*IDIR*DYI*(T(1,I,J+1)-T(1,I,J-1)) 900 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 1100 I=1,ND V(2,I,J)=V(1,I,J)+REX*(H(I)-V(1,I,J)) 1100 CONTINUE 1010 CONTINUE DO 1200 I=1,ND V(2,I,NN)=0.0 1200 CONTINUE DO 2150 I=1,ND DO 2150 J=1,NN V(1,I,J)=V(2,I,J) 2150 CONTINUE C C ********************************************************************* C FIND STREAM FUNCTION ON Y-DIRECTION LINES C ********************************************************************* C DO 3400 J=1,NN S(2,1,J)=S(1,1,J) 3400 CONTINUE DO 3600 I=2,NU1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 3500 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,J)-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 3500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN CALL TRISOL(NN,A,B,C,D,H) DO 3700 J=1,NN S(2,I,J)=S(1,I,J)+REX*(H(J)-S(1,I,J)) 3700 CONTINUE 3600 CONTINUE DO 4691 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 4591 K=NTP,NN1 J=K-NT+1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,K)-DX2I*(S(1,I+1,K)+S(1,I-1,K)) 4591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=LN CALL TRISOL(NTN,A,B,C,D,H) DO 4791 K=1,NTN J=NT+K-1 S(2,I,J)=S(1,I,J)+REX*(H(K)-S(1,I,J)) 4791 CONTINUE 4691 CONTINUE DO 4681 I=NPP,ND1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 7500 J=2,NN1 A(J)=-2.0*DX2I-2.0*DY2I B(J)=DY2I C(J)=DY2I D(J)=-V(1,I,J)-DX2I*(S(1,I+1,J)+S(1,I-1,J)) 7500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=LN CALL TRISOL(NN,A,B,C,D,H) DO 7700 J=1,NN S(2,I,J)=S(1,I,J)+REX*(H(J)-S(1,I,J)) 7700 CONTINUE 4681 CONTINUE DO 3800 J=1,NN S(2,ND,J)=S(2,ND-1,J) 3800 CONTINUE C C ********************************************************************* C DO 4100 I=1,ND DO 4100 J=1,NN S(1,I,J)=S(2,I,J) 4100 CONTINUE C C ********************************************************************* C FIND STREAM FUNCTION ON X-DIRECTION LINES C ********************************************************************* C DO 4410 I=1,ND S(2,I,1)=0.0 4410 CONTINUE DO 4010 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 4900 I=2,NU1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,I,J)-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 4900 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=S(1,NU,J) CALL TRISOL(NU,A,B,C,D,H) DIFFX=0.0 DO 4111 I=1,NU DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(I)-S(1,I,J)) 4111 CONTINUE 4010 CONTINUE DO 8010 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,NP,J) DO 6082 K=NPP,ND1 I=K-NP+1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,K,J)-DY2I*(S(1,K,J+1)+S(1,K,J-1)) 6082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DO 8111 K=1,NDN I=NP+K-1 DIFF=ABS(H(K)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(K)-S(1,I,J)) 8111 CONTINUE 8010 CONTINUE DO 6091 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=S(1,1,J) DO 6991 I=2,ND1 A(I)=-2.0*DX2I-2.0*DY2I B(I)=DX2I C(I)=DX2I D(I)=-V(1,I,J)-DY2I*(S(1,I,J+1)+S(1,I,J-1)) 6991 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 8191 I=1,ND DIFF=ABS(H(I)-S(1,I,J)) IF(DIFF.GT.DIFFX) DIFFX=DIFF S(2,I,J)=S(1,I,J)+REX*(H(I)-S(1,I,J)) 8191 CONTINUE 6091 CONTINUE DO 4200 I=1,ND S(2,I,NN)=LN 4200 CONTINUE DO 4150 I=1,ND DO 4150 J=1,NN S(1,I,J)=S(2,I,J) 4150 CONTINUE C C ********************************************************************* C TOP OF TEMPERATURE LOOP C ********************************************************************* C C C ********************************************************************* C FIND TEMPERATURE ON Y-DIRECTION LINES C ********************************************************************* C DO 6400 J=1,NN T(2,1,J)=0.0 6400 CONTINUE DO 6600 I=2,NU1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 6500 J=2,NN1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(T(1,I+1,J)-T(1,I-1,J)) $+PX2*(T(1,I+1,J)+T(1,I-1,J)) 6500 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 6700 J=1,NN T(2,I,J)=T(1,I,J)+REX*(H(J)-T(1,I,J)) 6700 CONTINUE 6600 CONTINUE DO 6696 I=NU,NP A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,I,NT) DO 6591 K=NTP,NN1 J=K-NT+1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-PY2 C(J)=+DXI*(S(1,I+1,K)-S(1,I-1,K))*DYI-PY2 D(J)=-DYI*(S(1,I,K+1)-S(1,I,K-1))*DXI*(T(1,I+1,K)-T(1,I-1,K)) $+PX2*(T(1,I+1,K)+T(1,I-1,K)) 6591 CONTINUE A(NTN)=1.0 B(NTN)=0.0 C(NTN)=0.0 D(NTN)=0.0 CALL TRISOL(NTN,A,B,C,D,H) DO 6799 K=1,NTN J=NT+K-1 T(2,I,J)=T(1,I,J)+REX*(H(K)-T(1,I,J)) 6799 CONTINUE 6696 CONTINUE DO 6681 I=NPP,ND1 A(1)=1.0 B(1)=-1.0 C(1)=0.0 D(1)=0.0 DO 6581 J=2,NN1 A(J)=2.0*PX2+2.0*PY2 B(J)=-DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 C(J)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI-PY2 D(J)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI*(T(1,I+1,J)-T(1,I-1,J)) $+PX2*(T(1,I+1,J)+T(1,I-1,J)) 6581 CONTINUE A(NN)=1.0 B(NN)=0.0 C(NN)=0.0 D(NN)=0.0 CALL TRISOL(NN,A,B,C,D,H) DO 6783 J=1,NN T(2,I,J)=T(1,I,J)+REX*(H(J)-T(1,I,J)) 6783 CONTINUE 6681 CONTINUE DO 6800 J=1,NN T(2,ND,J)=T(2,ND-1,J) 6800 CONTINUE C C ********************************************************************* C DO 6109 I=1,ND DO 6109 J=1,NN T(1,I,J)=T(2,I,J) 6109 CONTINUE C C ********************************************************************* C FIND TEMPERATURE ON X-DIRECTION LINES C ********************************************************************* C DO 6410 I=1,ND T(2,I,1)=T(1,I,1) 6410 CONTINUE DO 7101 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 7991 I=2,NU1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(T(1,I,J+1)-T(1,I,J-1)) $+PY2*(T(1,I,J+1)+T(1,I,J-1)) 7991 CONTINUE A(NU)=1.0 B(NU)=0.0 C(NU)=0.0 D(NU)=T(2,NU,J) CALL TRISOL(NU,A,B,C,D,H) DO 8127 I=1,NU T(2,I,J)=T(1,I,J)+REX*(H(I)-T(1,I,J)) 8127 CONTINUE 7101 CONTINUE DO 7192 J=2,NT A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=T(2,NP,J) DO 7082 K=NPP,ND1 I=K-NP+1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-PX2 C(I)=-DYI*(S(1,K,J+1)-S(1,K,J-1))*DXI-PX2 D(I)=+DXI*(S(1,K+1,J)-S(1,K-1,J))*DYI*(T(1,K,J+1)-T(1,K,J-1)) $+PY2*(T(1,K,J+1)+T(1,K,J-1)) 7082 CONTINUE A(NDN)=1.0 B(NDN)=0.0 C(NDN)=-1.0 D(NDN)=0.0 CALL TRISOL(NDN,A,B,C,D,H) DO 7374 K=1,NDN I=NP+K-1 T(2,I,J)=T(1,I,J)+REX*(H(K)-T(1,I,J)) 7374 CONTINUE 7192 CONTINUE DO 7010 J=NTP,NN1 A(1)=1.0 B(1)=0.0 C(1)=0.0 D(1)=0.0 DO 6900 I=2,ND1 A(I)=2.0*PX2+2.0*PY2 B(I)=DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 C(I)=-DYI*(S(1,I,J+1)-S(1,I,J-1))*DXI-PX2 D(I)=+DXI*(S(1,I+1,J)-S(1,I-1,J))*DYI*(T(1,I,J+1)-T(1,I,J-1)) $+PY2*(T(1,I,J+1)+T(1,I,J-1)) 6900 CONTINUE A(ND)=1.0 B(ND)=0.0 C(ND)=-1.0 D(ND)=0.0 CALL TRISOL(ND,A,B,C,D,H) DO 7100 I=1,ND T(2,I,J)=T(1,I,J)+REX*(H(I)-T(1,I,J)) 7100 CONTINUE 7010 CONTINUE DO 7200 I=1,ND T(2,I,NN)=0.0 T(2,I,1)=T(2,I,2) 7200 CONTINUE DO 6159 I=1,ND DO 6159 J=1,NN T(1,I,J)=T(2,I,J) 6159 CONTINUE C C ********************************************************************* C BOTTOM OF VORTICITY, STREAM FUNCTION, TEMPERATURE LOOP C ********************************************************************* C WRITE(6,5500) INTR,DIFFX 5500 format(' ITERATION NUMBER ',I5,' MAX. DIFF. IN STR. FUNC. ',F8.5) C C ********************************************************************* C IF(INTR.LT.300) GO TO 1000 IF(DIFFX.LT.0.00005) GO TO 2000 IF(INTR.GT.2500) GO TO 2000 GO TO 1000 2000 CONTINUE C C ********************************************************************* C WRITE OUT STREAM FUNCTION, VORTICITY AND TEMPERATURE VALUES C ********************************************************************* C WRITE(1,5931) 5931 format(/,' X Y S V $ T',/) DO 6149 I=1,ND DO 6149 J=1,NN X=(I-1)*(LU+1.0+LD)/(ND-1) Y=(J-1)*LN/(NN-1) WRITE(1,5287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) 5287 format(5F12.5) IF(I.LE.NU) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF IF(I.GE.NP) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF IF(I.GT.NU.AND.I.LT.NP) THEN IF(J.GE.NT) THEN WRITE(2,7287) X,Y,S(2,I,J),V(2,I,J),T(2,I,J) ENDIF ENDIF 7287 format(4(F12.5,','),F12.5) 6149 CONTINUE C C ********************************************************************* C WRITE OUT SURFACE VORTICITY AND HEAT TRANSFER RATE VALUES C ********************************************************************* C WRITE(1,5991) 5991 FORMAT(/,' POINT DISTANCE VORTICITY HEAT FLUX',/) IL=0 DO 5103 J=1,NT IL=IL+1 QL(IL)=(1.0-T(2,NU1,J))/DX XREL=0.0 WRITE(1,5200) IL,XREL,V(2,NU,J),(1.0-T(2,NU1,J))/DX 5103 CONTINUE DO 5100 I=NU,NP IL=IL+1 QL(IL)=(1.0-T(2,I,NTP))/DY XREL=REAL(I-NU)/REAL(NP-NU) WRITE(1,5200) IL,XREL,V(2,I,NT),(1.0-T(2,I,NTP))/DY 5200 format(I4,3F12.5) 5100 CONTINUE DO 5106 J=1,NT IL=IL+1 K=NT-J+1 QL(IL)=(1.0-T(2,NPP,K))/DX XREL=1.0 WRITE(1,5200) IL,XREL,V(2,NP,K),(1.0-T(2,NPP,K))/DX 5106 CONTINUE C C ********************************************************************* C DETERMINE AND WRITE OUT MEAN HEAT TRANSFER RATE C ********************************************************************* C QM=QL(1)*DY/2.0 DO 4928 J=2,NT1 QM=QM+QL(J)*DY 4928 CONTINUE QM=QM+QL(NT)*DY/2.0 QM=QM+QL(NTP)*DX/2.0 DO 4927 I=NUP,NP1 K=I-NU+1+NT QM=QM+QL(K)*DX 4927 CONTINUE QM=QM+QL(NP-NU+1+NT)*DX/2.0 QM=QM+QL(NP-NU+2+NT)*DY/2.0 DO 4926 J=2,NT1 K=NT+(NP-NU+1)+J QM=QM+QL(K)*DY 4926 CONTINUE QM=QM+QL(2*NT+NP-NU+1)*DY/2.0 write(1,5283) QM/2.0 5283 format(/,' MEAN NUSSELT NUMBER = ',F12.5) GR=IDIR*GR write(6,5247) RE,PR,GR 5247 format(/,' REYNOLDS NUMBER= ', F12.1,/,' PRANDTL NUMBER=',F12.3, & /,' GRASHOF NUMBER= ', F12.1,/) write(6,5283) QM/2.0 close(1) close(2) stop * ************************** FORMAT STATEMENTS ************************** * 1002 format (10X,'MIXED CONVECTIVE FLOW OVER A SQUARE CYLINDER') 1003 format (10X,'--------------------------------------------',//) 1040 format(/////,5X,'Input the Value of the Reynolds Number ',/) 1050 format(//,5X,'Input the Value of the Prandtl Number ',/) 1102 format(//,5X,'Input the Value of the Grashof Number ',/) 1400 format(////////,5X,'MIXED CONVECTION OVER A SQUARE CYLINDER', & /,5X,'---------------------------------------',///, & /,5X,' Output is written to a file named MIXSQCPR.DAT',///) 1500 format(//,5X,'Assisting Flow (+1) or Opposing Flow (-1) ',/) END C C ********************************************************************* C END OF MAIN PROGRAM C ********************************************************************* C C C ********************************************************************* C SUBROUTINE TRISOL(N,A,B,C,D,H) C C ********************************************************************* C TRI-DIAGONAL MATRIX SOLVER C ********************************************************************* C REAL A(91),B(91),C(91),D(91),H(91),W(91),R(91),G(91) W(1)=A(1) G(1)=D(1)/W(1) DO 100 I=2,N I1=I-1 R(I1)=B(I1)/W(I1) W(I)=A(I)-C(I)*R(I1) G(I)=(D(I)-C(I)*G(I1))/W(I) 100 CONTINUE H(N)=G(N) N1=N-1 DO 200 I=1,N1 II=N-I H(II)=G(II)-R(II)*H(II+1) 200 CONTINUE RETURN END