-
Notifications
You must be signed in to change notification settings - Fork 0
/
DEWP.f90
247 lines (201 loc) · 7.17 KB
/
DEWP.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
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
SUBROUTINE BULBP ()
IMPLICIT NONE
END SUBROUTINE
PROGRAM BULBP
IMPLICIT NONE
REAL :: X1 = 0.0, X2 = 0.0, DELTAX1 = 0.01, Y1, Y2, ERROR = 0.01
REAL :: T, PHI1, PHI2, PSAT1, PSAT2, GAMMA1, GAMMA2, PA, PA2
REAL :: GAMMA12, GAMMA22
!ANTOINE EQUATION CONSTANTS
REAL :: A1 = 73.304, B1 = -7122.3, C1 = -7.1424, D1 = 2.8853E-6
REAL :: A2 = 107.02, B2 = -10237.0, C2 = -11.695, D2 = 6.8003E-18
INTEGER :: E1 = 2, E2 = 6
!VAN LAAR CONSTANTS
REAL :: A = 489.38, B = 1173.28
REAL :: R = 8.314
!PENG ROBINSON VARIABLES
REAL :: ACUB1, BCUB1, TC1 = 514.29, PC1 = 6137000, AC1, BZ1, W1 = 0.649
REAL :: ACUB2, BCUB2, TC2 = 577.2, PC2 = 3.9013E+06, AC2, BZ2, W2 = 0.5890
!CUBIC RESOLUTION VARIABLES
REAL :: CUB11, CUB21, CUB31, Z11, Z21, Z31, Z2I1, Z3I1, Z1
REAL :: CUB12, CUB22, CUB32, Z12, Z22, Z32, Z2I2, Z3I2, Z2
!FIXED TEMPERATURE IN K
T = 343
!PRINT FORMAT
PRINT*, "PROCESSING..."
OPEN(10, FILE = "OUTPUT.TXT")
WRITE(10,*) "CONSTANT TEMPERATURE:", T, "K"
WRITE(10,*) "P(PA) X1 X2 Y1 Y2"
!DETERMINATION OF VAPOR PRESSURE
PSAT1 = EXP(A1 + (B1/T) + C1*LOG(T) + D1*T**E1)
PSAT2 = EXP(A2 + (B2/T) + C2*LOG(T) + D2*T**E2)
DO WHILE (Y1 < 1)
Y2 = 1 - Y1
!INITIALIZE WITH FUGACITY COEFFICIENT EQUAL TO 1
PHI1 = 1
PHI2 = 1
!INITIALIZE WITH ACTIVITY COEFFICIENT EQUAL TO 1
GAMMA1 = 1
GAMMA2 = 1
!DETERMINATION OF P
PA = 1/((Y1*PHI1)/(GAMMA1*PSAT1)+(Y2*PHI2)/(GAMMA2*PSAT2))
!DETERINATION OF CONCENTRATIONS IN LIQUID-PHASE
X1 = (Y1*PHI1*PA)/(GAMMA1*PSAT1)
X2 = (Y2*PHI2*PA)/(GAMMA2*PSAT2)
!DETERMINATION OF ACTIVITY COEFFICIENT
GAMMA1 = EXP((A*((B*X2)/(A*X1 + B*X2))**2)/(R*T))
GAMMA2 = EXP((B*((A*X1)/(A*X1 + B*X2))**2)/(R*T))
DO WHILE (SQRT(PA-PA2)**2 > ERROR)
PA = PA2
DO WHILE (SQRT(GAMMA1-GAMMA12)**2 > ERROR .OR. SQRT(GAMMA2-GAMMA22)**2 > ERROR)
GAMMA1 = GAMMA12
GAMMA2 = GAMMA22
!DETERMINATION OF FUGACITY COEFFICIENT
!TRANSFORMATION OF VALUES TO CUBIC FORMAT
CALL PRECUBIC (T, TC1, PC1, W1, PA, R, CUB11, CUB21, CUB31, ACUB1, BCUB1, AC1, BZ1)
CALL PRECUBIC (T, TC2, PC2, W2, PA, R, CUB12, CUB22, CUB32, ACUB2, BCUB2, AC2, BZ2)
!DETERMINATION OF CUBIC ROOTS
CALL SOLVECUBIC (CUB11, CUB21, CUB31, Z11, Z21, Z31, Z2I1, Z3I1)
CALL SOLVECUBIC (CUB12, CUB22, CUB32, Z12, Z22, Z32, Z2I2, Z3I2)
!CHOOSING CORRECT ROOTS
CALL CHOOSEZ (Z11, Z21, Z31, Z2I1, Z3I1, ACUB1, BCUB1, Z1)
CALL CHOOSEZ (Z12, Z22, Z32, Z2I2, Z3I2, ACUB2, BCUB2, Z2)
!DETERINATION OF CONCENTRATIONS IN LIQUID-PHASE
X1 = (Y1*PHI1*PA)/(GAMMA1*PSAT1)
X2 = (Y2*PHI2*PA)/(GAMMA2*PSAT2)
!NORMALIZATION OF X VALUES
X1 = X1/(X1+X2)
X2 = X2/(X1+X2)
GAMMA12 = EXP((A*((B*X2)/(A*X1 + B*X2))**2)/(R*T))
GAMMA22 = EXP((B*((A*X1)/(A*X1 + B*X2))**2)/(R*T))
END DO
!DETERMINATION OF NEW PRESSURE VALUE
PA2 = (X1*GAMMA1*PSAT1)/PHI1 + (X2*GAMMA2*PSAT2)/PHI2
!DETERMINATION OF CONCENTRATIONS IN VAPOR-PHASE
Y1 = X1*GAMMA1*PSAT1/(PHI1*PA)
Y2 = X2*GAMMA2*PSAT2/(PHI2*PA)
END DO
!PRINT FORMAT
WRITE(10,10) PA, X1, X2, Y1, Y2
10 FORMAT (F12.4, F12.4,F12.4,F12.4,F12.4)
!NEXT CONCENTRATION IN LIQUID-PHASE
X1 = X1 + DELTAX1
END DO
END PROGRAM BULBP
SUBROUTINE PRECUBIC (T, TC, PC, W, P, R, CUB1, CUB2, CUB3, ACUB, BCUB, AC, BZ)
IMPLICIT NONE
REAL, INTENT(IN) :: T, TC, PC, W, P, R
REAL, INTENT(OUT) :: CUB1, CUB2, CUB3, ACUB, BCUB, AC, BZ
REAL :: TR, AT, K
IF (W .LE. 0.49) THEN
K = 0.3746 + 1.54226*W - 0.26992*W**2
ELSE
K = 0.379642 + 1.48503*W -0.164423*W**2 + 0.016666*W**3
ENDIF
TR = T/TC
AC = 0.45724*R*R*TC*TC/PC
BZ = 0.07780*R*TC/PC
AT = AC*(1 + K*(1 - SQRT(TR)))**2
ACUB = AT*P/((R*T)*(R*T))
BCUB = BZ*P/(R*T)
CUB1 = BCUB - 1
CUB2 = ACUB - 2*BCUB - 3*BCUB*BCUB
CUB3 = BCUB*BCUB*BCUB + BCUB*BCUB - ACUB*BCUB
RETURN
END SUBROUTINE PRECUBIC
SUBROUTINE SOLVECUBIC (C1 , C2 , C3, X1, X2, X3, X2I, X3I)
REAL, INTENT(IN) :: C1 , C2 , C3
REAL :: Q, R, THETA, A, B
REAL, INTENT(OUT) :: X1, X2, X3, X2I, X3I
Q = (C1*C1 - 3*C2)/9
R = (2*C1*C1*C1 - 9*C1*C2 + 27*C3)/54
IF (R**2 .LT. Q**3) THEN
THETA = ACOS(R/SQRT(Q**3))
X1 = -2*SQRT(Q)*COS(THETA/3)-C1/3
X2 = -2*SQRT(Q)*COS((THETA+2*3.1415)/3)-C1/3
X3 = -2*SQRT(Q)*COS((THETA-2*3.1415)/3)-C1/3
X2I = 0
X3I = 0
ELSE
A = -(R + SQRT(R*R - Q*Q*Q))**(1/3)
IF (A .EQ. 0) THEN
B = Q/A
ELSE
B = 0
ENDIF
X1 = (A + B) - C1/3
X2 = -(1/2)*(A + B) - C1/3
X2I = SQRT(3.0)*(A - B)/2
X3 = -(1/2)*(A + B) - C1/3
X2I = -SQRT(3.0)*(A - B)/2
ENDIF
RETURN
END SUBROUTINE SOLVECUBIC
SUBROUTINE CALCPHI (A1, A2, X1, X2, B1, B2, Z1, Z2, ACUB1, BCUB1, ACUB2, BCUB2, PHI1, PHI2)
REAL, INTENT(IN) :: A1, A2, X1, X2, B1, B2, Z1, Z2, ACUB1, BCUB1, ACUB2, BCUB2
REAL, INTENT(OUT) :: PHI1, PHI2
REAL :: A, B, PHI11, PHI12, PHI21, PHI22
A = A1*X1**2 + 2*X1*X2*SQRT(A1*A2) + A2*X2**2
B = X1*B1 + X2*B2
PHI11 = -LOG(Z1-BCUB1) + (Z1-1)*B1/B - (ACUB1/(2*SQRT(2.0)*BCUB1))*((1/A)*(2*(X1*A1 + X2*SQRT(A1*A2))) - (B1/B))
PHI12 = LOG((Z1 + (SQRT(2.0) + 1.0)*BCUB1)/(Z1-(SQRT(2.0) - 1.0)*BCUB1))
PHI1 = EXP(PHI11*PHI12)
PHI21 = -LOG(Z2-BCUB2) + (Z2-1)*B2/B - (ACUB2/(2*SQRT(2.0)*BCUB2))*((1/A)*(2*(X1*SQRT(A1*A2) + X2*A2)) - (B2/B))
PHI22 = LOG((Z2 + (SQRT(2.0) + 1.0)*BCUB2)/(Z2-(SQRT(2.0) - 1.0)*BCUB2))
PHI2 = EXP(PHI21*PHI22)
RETURN
END SUBROUTINE CALCPHI
SUBROUTINE CHOOSEZ (Z1, Z2, Z3, Z2I, Z3I, ACUB, BCUB, Z)
REAL, INTENT(IN) :: Z1, Z2, Z3, Z2I, Z3I, ACUB, BCUB
REAL, INTENT(OUT) :: Z
REAL :: ZMIN, ZMAX, FUGMIN, FUGMAX, FUGRATIO
IF (Z2I .NE. 0 .AND. Z3I .NE. 0) THEN
Z = Z1
ELSE IF (Z2I .NE. 0 ) THEN
ZMIN = MIN(Z1, Z3)
ZMAX = MAX(Z1, Z3)
IF (ZMIN < 0) THEN
Z = ZMAX
ELSE
FUGMIN = ((ZMIN - 1)/(ZMIN-BCUB))*((ZMIN+(SQRT(2.0)+1)*BCUB)/(ZMIN-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGMAX = ((ZMAX - 1)/(ZMAX-BCUB))*((ZMAX+(SQRT(2.0)+1)*BCUB)/(ZMAX-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGRATIO = FUGMAX/FUGMIN
IF(FUGRATIO>1) THEN
Z = ZMIN
ELSE
Z = ZMAX
ENDIF
ENDIF
ELSE IF (Z3I .NE. 0) THEN
ZMIN = MIN(Z1, Z2)
ZMAX = MAX(Z1, Z2)
IF (ZMIN < 0) THEN
Z = ZMAX
ELSE
FUGMIN = ((ZMIN - 1)/(ZMIN-BCUB))*((ZMIN+(SQRT(2.0)+1)*BCUB)/(ZMIN-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGMAX = ((ZMAX - 1)/(ZMAX-BCUB))*((ZMAX+(SQRT(2.0)+1)*BCUB)/(ZMAX-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGRATIO = FUGMAX/FUGMIN
IF(FUGRATIO>1) THEN
Z = ZMIN
ELSE
Z = ZMAX
ENDIF
ENDIF
ELSE
ZMIN = MIN(Z1, Z2, Z3)
ZMAX = MAX(Z1, Z2, Z3)
IF (ZMIN < 0) THEN
Z = ZMAX
ELSE
FUGMIN = ((ZMIN - 1)/(ZMIN-BCUB))*((ZMIN+(SQRT(2.0)+1)*BCUB)/(ZMIN-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGMAX = ((ZMAX - 1)/(ZMAX-BCUB))*((ZMAX+(SQRT(2.0)+1)*BCUB)/(ZMAX-(SQRT(2.0)+1)*BCUB))**(ACUB/(2*SQRT(2.0)*BCUB))
FUGRATIO = FUGMAX/FUGMIN
IF(FUGRATIO>1) THEN
Z = ZMIN
ELSE
Z = ZMAX
ENDIF
ENDIF
ENDIF
RETURN
END SUBROUTINE CHOOSEZ