This is a listing of all the programs in the order given in the list.
NOTE: If you download this file you will have to separate the programs using a text editor.
************************************************************************
***** *****
***** 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