************************************************************************
***** *****
***** 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