forked from navoj/ChemEngNumMeth
-
Notifications
You must be signed in to change notification settings - Fork 1
/
ex54.for
253 lines (253 loc) · 8.1 KB
/
ex54.for
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
C
C*********************** ABSTRACT ********************************
C
C THIS PROGRAM SOLVES THE FLOW DISTRIBUTION IN A DIVIDING FLOW
C MANIFOLD. USING THE DEVELOPMENT BY BAJURA AND JONES(7), THE FLOW
C IN A DIVIDING FLOW MANIFOLD IS DESCRIBED BY A SECOND ORDER NON-
C LINEAR DIFFERENTIAL WITH TWO CONSTANT VALUED BOUNDARY CONDITIONS.
C THIS BOUNDARY VALUE PROBLEM IS SOLVED USING FINITE DIFFEENCE
C APPROXIMATIONS WITH A SUCCESSIVE RELAXATION CONVERGENCE SCHEME.
C
C
C********************** NOMENCLATURE ******************************
C
C DX- DIMENSIONLESS SEPARATION BETWEEN GRID POINTS
C D1- MANIFOLD DIAMETER (M)
C D3- DIAMETER OF THE OPENING TO THE LATERAL PATH (M)
C ERLIM- CONVERGENCE CRITERIA (RELATIVE ERROR)
C FF- FRICTION FACTOR FOR THE MANIFOLD (DIMENSIONLESS)
C H- LATERAL RESISTANCE COEFICIENT IN VELOCITY HEADS (DIMENSIONLESS)
C ITER- ITERATION NUMBER
C L- LENGTH OF THE MANIFOLD (M)
C XM1- PARAMETER IN MANIFOLD FLOW EQUATION WHICH CONTAINS THE PRESSURE
C REOVERY EFFECT
C MAXITR- MAXIMUM NUMBER OF ITERATIONS
C PHI- PARAMETER IN THE MANIFOLD EQUATION WHICH CONTAINS THE FRICTION
C EFFECT
C SRF- SUCCESSIVE RELAXATION FACTOR
C Q(I)- DIMENSIONLESS MANIFOLD FLOW FOR GRID POINT I
C QLAT- DIMENSIONLESS MANIFOLD FLOW INTO A LATERAL PATH
C TEST- RELATIVE ERROR FOR A GRID POINT
C THETA- OVERALL MOMENTUM CORRETION FACTOR (DIMENSIONLESS)
C
C************************************************************************
C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION Q(30)
REAL L,M1
COMMON /ONE/XM1,PHI
EXTERNAL RECUR,BCL,BCR
C ****************
C********************** INPUT DATA ******************************
C ****************
C************** PHYSICAL PARAMETERS
C
H=1.
N=11
XN=FLOAT(N)-1.0
DX=1./FLOAT(N)
D3=.005
D1=.02
FF=.003
L=.5
THETA=1.1
PHI=XN*XN*(D3/D1)**4*FF*L/(H*2.*D1)
XM1 =XN*XN*(D3/D1)**4*THETA/H
C
C***************** NUMERICAL PARAMETERS
C
ERLIM=1.E-5
MAXITR=100
SRF=1.5
YG=.5
C
C***************** CALL SUBROUTINE SOR
C
XL=0.0
XR=1.0
DO 33 I=1,N
33 Q(I)=1.-DFLOAT(I-1)/DFLOAT(N-1)
CALL SOR(RECUR,BCL,BCR,N,SRF,ERLIM,MAXITR,XL,XR,Q,ITER)
C ***********************
C********************* PRINT OUT RESULTS *************************
C ***********************
WRITE(6,444)ITER
444 FORMAT( 5X,' AFTER ',I4,' ITERATIONS THE CONVERGED RESULTS ARE:')
WRITE(6,306)
306 FORMAT( //)
WRITE(6,305)
305 FORMAT( 5X,12HLATERAL PATH,5X,15HNORMALIZED FLOW)
C
C*********** CALCULATE NORMALIZED LATERAL FLOW
C
NM=N-1
DO 300 I=1,NM
QLAT=Q(I)-Q(1+I)
C PRINT OUT THE NORMALIZED LATERAL FLOWS
300 WRITE(6,301)I,QLAT
301 FORMAT( 9X,I3,12X,E11.4)
STOP
END
C
C***************************** ABSTRACT ******************************
C
C THIS SUBROUTINE APPLIES THE RECURSION RELATION FOR THE GOVERNING
C DIFFERENTIAL EQUATION. THIS SUBROUTINE IS USER SUPPLIED.
C
C*************************** NOMENCLATURE *****************************
C
C A- COEFICIENT IN THE SQUARE TERM IN THE RECURSION RELATION
C B- COEFICIENT IN THE LINEAR TERM IN THE RECURSION RELATION
C C- CONSTANT TERM IN THE RECURSION RELATION
C YC- THE CALCULATED VALUE OF Y(I) BASED UPON THE RECURSION RELATION
C YM- THE VALUE OF Y(I-1)
C YP- THE VALUE OF Y(I+1)
C DX- THE SEPARATION BETWEEN NODE POINTS
C
C*************************************************************************
C
SUBROUTINE RECUR(YM,YC,YP,DX)
IMPLICIT REAL*8(A-H,O-Z)
COMMON /ONE/XM1,PHI
A=PHI
B=(YP-YM)/(2.*DX)*(XM1-2./DX/DX)
C=(YP*YP-YM*YM)/(2.*DX*DX*DX)
C YC=(-B+DSQRT(B*B-4.*A*C))/2./A
YC=-C/(A*YC+B)
RETURN
END
C
C***************************** ABSTRACT ******************************
C
C THIS SUBROUTINE APPLIES THE RECURSION RELATION FOR THE LEFT
C BOUNDARY CONDITION. THIS SUBROUTINE IS USER SUPPLIED.
C
C*************************** NOMENCLATURE *****************************
C
C Y1- THE CALCULATED VALUE OF Y(1) BASED UPON THE RECURSION RELATION
C Y2- THE VALUE OF Y(2)
C Y3- THE VALUE OF Y(3)
C DX- THE SEPARATION BETWEEN NODE POINTS
C
C*************************************************************************
C
SUBROUTINE BCL(Y1,Y2,Y3,DX)
IMPLICIT REAL*8(A-H,O-Z)
Y1=1.0
RETURN
END
C
C***************************** ABSTRACT ******************************
C
C THIS SUBROUTINE APPLIES THE RECURSION RELATION FOR THE RIGHT
C BOUNDARY CONDITION. THIS SUBROUTINE IS USER SUPPLIED.
C
C*************************** NOMENCLATURE *****************************
C
C YN- THE CALCULATED VALUE OF Y(N) BASED UPON THE RECURSION RELATION
C YNN- THE VALUE OF Y(N-1)
C YNNN- THE VALUE OF Y(N-2)
C DX- THE SEPARATION BETWEEN NODE POINTS
C
C*************************************************************************
C
SUBROUTINE BCR(YN,YNN,YNNN,DX)
IMPLICIT REAL*8(A-H,O-Z)
YN=0.0
RETURN
END
C
C*********************** ABSTRACT ********************************
C
C THIS SOLVES ONE-DIMENSIONAL BVP'S USING THE SOR METHOD. THE USER
C MUST SUPPLY SUBROUTINES RECUR, BCL, AND BCR WHICH APPLY THE RECURSION
C RELATIONS FOR THE INTERIOR NODES, THE LEFT BOUNDARY, AND THE RIGHT
C BOUNDARY. ADDITIONAL INPUT PARAMETERS ARE SUPPLIED AS ARGUEMENTS
C TO THE CALL OF SUBROUTINE SOR.
C
C********************** NOMENCLATURE ******************************
C
C DX- SEPARATION BETWEEN GRID POINTS (M)
C ERLIM- CONVERGENCE CRITERIA (RELATIVE ERROR)
C ITER- ITERATION NUMBER
C N- NUMBER OF GRID POINTS
C MAXITR- MAXIMUM NUMBER OF ITERATIONS USED IN SEEKING A CONVERGED
C SOLUTION
C NUMERS- NUMBER OF POINTS THAT HAVE NOT MET CONVERGENCE CRITERIA FOR
C A PARTIULAR ITERATION
C SRF- SUCCESSIVE RELAXATION FACTOR
C TEST- RELATIVE ERROR FOR A GRID POINT
C Y(I)- THE DEPENDENT VARIABLE VALUE AT THE ITH NODE
C YC- THE DEPENDENT VARIABLE VALUE CALCULATED BY THE RECURSION RELATION
C YG- INITIAL GUESS FOR ALL THE GRID POINTS (DEG C)
C TOTAL- TOTAL ERROR
C
C
SUBROUTINE SOR(RECUR,BCL,BCR,N,SRF,ERLIM,MAXITR,XL,XR,Y,ITER)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION Y(1)
EXTERNAL RECUR,BCL,BCR
DX=(XR-XL)/FLOAT(N-1)
ITER=0
C *******************************
C******************* BEGIN ITERATIVE PROCEDURE *********************
C *******************************
100 ITER=ITER+1
NUMERS=0
TOTAL=0.0
C
C************* UPDATE BOUNDARY POINTS
C
Y1=Y(1)
Y2=Y(2)
Y3=Y(3)
C UPDATE LEFT BOUNDARY CONDITION
CALL BCL(Y1,Y2,Y3,DX)
Y(1)=Y1
YN=Y(N)
YNN=Y(N-1)
YNNN=Y(N-2)
C UPDATE THE RIGHT BOUNDARY CONDITION
CALL BCR(YN,YNN,YNNN,DX)
Y(N)=YN
NM=N-1
C
C************** UPDATE INTERIOR POINTS
C
DO 20 I=2,NM
YP=Y(I+1)
YM=Y(I-1)
YC=Y(I)
CALL RECUR(YM,YC,YP,DX)
TEST=ABS((YC-Y(I))/YC)
TOTAL=TOTAL+DABS(YC-Y(I))
Y(I)=Y(I)+SRF*(YC-Y(I))
20 IF(TEST.GT.ERLIM)NUMERS=NUMERS+1
IF(ITER.GT.MAXITR)GO TO 200
IF(NUMERS.NE.0)GO TO 100
C
C CONVERGED SOLUTION; RETURN TO MAIN PROGRAM
C
RETURN
200 CONTINUE
C ***************************************
C*************** PRINT OUT RESULTS IF NOT CONVERGED ******************
C ***************************************
202 FORMAT( ///)
WRITE(6,202)
WRITE(6,201)ITER,NUMERS,TOTAL
201 FORMAT( 7H AFTER ,I3,22H ITERATIONS THERE ARE ,I3,61H POINTS THAT
1HAVE NOT CONVERGED. THE TOTAL RELATIVE ERROR IS,E10.3)
C*********** PRINT OUT DEPENDENT VARIABLE VALUES
WRITE(6,202)
WRITE(6,203)
203 FORMAT( ' X Y')
WRITE(6,209)
209 FORMAT( /)
X=0.0
DO 210 I=1,N
WRITE(6,204)X,Y(I)
210 X=X+DX
204 FORMAT( 5X,E11.4,5X,E12.5)
RETURN
END