-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBL_MAIN.f90
149 lines (109 loc) · 5.59 KB
/
BL_MAIN.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
!****************************************************************************
! * FILE = MULTIPLE FILES *
! * AUTHOR = CHANDAN KUMAR ([email protected]) *
! * INSTITUTE = INDIAN INSTITUTE OF TECHNOLOGY (IIT), KHARAGPUR *
!****************************************************************************
! * THIS PROGRAM IS DISTRIBUTED IN A HOPE THAT IT WILL BE USEFUL. *
!****************************************************************************
PROGRAM COMPRESSIBLE_BL
USE GLOBAL
IMPLICIT NONE
DOUBLE PRECISION :: C1, C2, LD, P1, Q1, C11, C22, NAD
INTEGER :: MA1, NA
TTOT = 311.0D0 ! TOTAL TEMPERATURE (GIVEN)
TINF = TTOT/(1.0D0 + 0.5D0*(GAMMAN -1.0D0)*M**2) ! STATIC TEMPERATURE
SOS =SQRT(GAMMAN * RCONST * TINF) ! SPEED OF SOUND
UINF = SOS*M ! FREESTREAM VELOCITY
IF (TINF .GT. 110.4D0) THEN
FMU =(1.458D0*TINF**1.5)/(TINF+110.4D0)*1.0D-05 ! FREESTREAM VISCOSITY
ELSE
FMU = (0.693873D0*1.0D-06)*TINF
END IF
DO I = 1, N
T(I) = 1.0D0 ! NONDIMENSIONAL TEMPERATURE
MU(I) = 1.0D0 ! NONDIMENSIONAL VISCOSITY
DMU(I) = 0.0D0 ! DERIVATIVE OF NONDIMENISONAL VISCOSITY
END DO
!====================================================================================
!====================================================================================
DO NA =1, 40
!-----------SOLVING BL MOMENTUM EQUATION---------------------------------------------
! BC
Y(1,1) = 0.0D0 ! U AT SURFACE
Y(3,1) = 0.0D0 ! G AT SUARFCE
P = 0.01D0 ! TWO INITIAL GUESSES FOR U' AT SURFACE TO USE SECANT ITERATION
Q = 0.02D0
Y(2,1) =P
CALL RK41
C1 = Y(1,N)-1.0D0 ! Y(1,N) IS U IN FREESTREAM THAT MUST SATISFY 1.
! HENCE, FOR PERFECT INITIAL GUESS C1 OR C2 MUST BE ZERO.
Y(2,1) = Q
CALL RK41
C2 = Y(1,N)-1.0D0
DO MA = 1, 10 ! SECANT ITERATION STARTS WITH INITIAL GUESS
R = P - (C1 *(P-Q))/(C1- C2)
C2 = C1
Q = P
Y(2,1) = R
CALL RK41
C1 = Y(1,N)-1.0D0
P = R
IF (DABS(C1) .LT. 1.0D-10) THEN
EXIT
END IF
END DO
!-----------SOLVING BL ENERGY EQUATION---------------------------------------------
! BC
Y(5,1) = 0.0D0 ! THETA' (DERIVATIVE OF SCALED TEMPERATURE) AT SURFACE
P1 = 0.4D0 ! THETA (SCALED TEMPERATURE) AT SUFACE AS INTIAL GUESS
Q1 = 0.5D0
Y(4,1) = P1
CALL RK42
C11 = Y(4,N) ! Y(4,N) SCALED TEMPRATURE IN FREESTREAM MUST SATISFY ZERO
Y(4,1) = Q1
CALL RK42
C22 = Y(4,N)
DO MA1 = 1, 10 ! SECANT ITERATION STARTS WITH INITIAL GUESS
R = P1-(C11*(P1-Q1))/(C11-C22)
C22 = C11
Q1 = P1
Y(4,1) = R
CALL RK42
C11 = Y(4,N)
P1 = R
IF (DABS(C11) .LT. 1.0D-10) THEN
EXIT
END IF
END DO
DO J =1, N ! UPDATE VISCOITY USING UPDATED TEMPERATURE
T(J) = 1+0.5D0*(GAMMAN-1.0D0)*M*M*Y(4,J)
TINST = T(J)*TINF
IF (TINST .GT. 110.4D0) THEN
WMU =(1.458D0*TINST**1.5)/(TINST+110.4D0)*1.0D-05
ELSE
WMU = (0.693873D0*1.0D-06)*TINST
END IF
MU(J) = WMU/FMU
END DO
DO J = 1, N ! UPDATE VISCOITY DERIVATIVE
IF ( J .EQ. 1) THEN
DMU(J)= (-3.0D0*MU(J) + 4.0D0*MU(J+1) - MU(J+2))/(2.0D0*H)
ELSE IF (J .EQ. N) THEN
DMU(J) = (3.0D0*MU(J)-4.0D0*MU(J-1)+MU(J-2))/(2.0D0*H)
ELSE
DMU(J) = (MU(J+1)-MU(J-1))/(2.0D0*H)
END IF
END DO
PRINT*, 'ITERATION NO.' , NA
END DO
!================================================================================
!================================================================================
OPEN(2,FILE="DATA.DAT")
WRITE(2,*),'NON-DIMENSIONAL U,', 'NON-DIMENSIONAL T'
PRINT*, 'NON-DIMENSIONAL U,', 'NON-DIMENSIONAL T'
DO J =1, N
T(J) = 1+0.5D0*(GAMMAN-1.0D0)*M*M*Y(4,J)
PRINT*, Y(1,J), T(J)
WRITE(2,*),Y(1,J), T(J)
END DO
END PROGRAM COMPRESSIBLE_BL