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