2007_elhamsadeghi
جمعه 30 بهمن 1388, 16:16 عصر
با سلام به همگی
من یک دانشجوی کارشناسی ارشد مهندسی شیمی هستم.پروژهای که محاسبات فلش برای یک سیستم دو جزئی انجام می دهد را با VB6 نوشتم که همه Subroutine ها در یک ماژول تعریف کردم. نمی دونم چرا در برنامه ای که در فرم اصلی برنامه نوشتم به حلقه while ایراد میگیره؟
نا آخر روز شنبه هم بیشتر فرصت ندارم!کسی میتونه کمکم کنه؟
خطایی که میده:Wend whitout While
در ضمن من چون قبلا برنامه رو با متلب نوشتم،مدل متلب نوشتن رو در برنامه VB پاک نکردم و بصورت توضیحات آوردم.
برنامه اصلی :
Public z1, z2, m1, m2, err, VF, Psat1, Psat2, Pdew, Pbub, Pfac, xdew1, xdew2, x1_old, ybub1, ybub2, y1_old, gam1, gam2, gam_dew1, gam_dew2, gam_bub1, gam_bub2, phi1, phi2, phi_dew1, phi_dew2, phi_bub1, phi_bub2, K1, K2, tol As Double
Public T, V1, V2, Rou1, Rou2 As Double
Public P, iter As Integer
Private Sub Flash_Bottom_Click()
' 1:Water 2:Vinylacetate
'----------------------------------
'----------------------------------
' Input Data
T = 80 + 273.15 'K
P = 1 'bar
V1 = 3 'Lit
V2 = 1.1 'Lit
'----------------------------------
'----------------------------------
' Calculating Mol Fraction of the feed (zi)
'------------------------------------------
' Water :1
V1 = V1 / 1000 'm^3
Rou1 = 971.8 'Kg/m^3
'Mw1 = 18.01528; %Kg/Kmol
m1 = V1 * Rou1 'Kg
' Vinyacetate :2
V2 = V2 / 1000 'm^3
Rou2 = 934# 'kg/m^3
'Mw2 = 86.09044; % Kg/Kmol
m2 = V2 * Rou2 'Kg
z1 = m1 / (m1 + m2)
z2 = 1 - z1
'----------------------------------
'----------------------------------
' Checking the State of mixture
'----------------------------------
Call DewP(z1, T, Pdew, xdew1, xdew2)
'[Pdew,xdew1,xdew2] = DewP(z1,T); % Dew Point Pressure an xi
Call Bublp(z1, T, Pbub, ybub1, ybub2)
'[Pbub,ybub1,ybub2] = BublP(z1,T); % Bubble Point Pressure an yi
If P > Pbub Then
Text1(0).Text = "iquid Phase"
Else
If P < Pdew Then
Text1(0).Text = "Vapor Phase"
Else
'----------------------------------
'----------------------------------
'Performing VLE Flash Calculation
'----------------------------------
' Calculating properties at Dew and bubble Point
Call Gamma_Wat(xdew1, T, gam_dew1)
'gam_dew1 = Gamma_Wat(xdew1,T);
Call Gamma_Vin(xdew2, T, gam_dew2)
'gam_dew2 = Gamma_Vin(xdew2,T);
Call Gamma_Wat(ybub1, T, gam_bub1)
'gam_bub1 = Gamma_Wat(ybub1,T);
Call Gamma_Vin(ybub2, T, gam_bub2)
'gam_bub2 = Gamma_Vin(ybub2,T);
Call CalcPhi(T, Pdew, xdew1, phi_dew1, phi_dew2)
'[phi_dew1,phi_dew2] = CalcPhi(T,Pdew,xdew1);
Call CalcPhi(T, Pbub, ybub1, phi_bub1, phi_bub2)
'[phi_bub1,phi_bub2] = CalcPhi(T,Pbub,ybub1);
'Initial Guess by Intrpolation
'initializing by means of Lever Rule
Pfac = (P - Pdew) / (Pbub - Pdew)
gam1 = Pfac * (gam_bub1 - gam_dew1) + gam_dew1
gam2 = Pfac * (gam_bub2 - gam_dew2) + gam_dew2
phi1 = Pfac * (phi_bub1 - phi_dew1) + phi_dew1
phi2 = Pfac * (phi_bub2 - phi_dew2) + phi_dew2
VF = (Pbub - P) / (Pbub - Pdew) 'Initial Guess for V/F
'--------------------------------------------
Call CalcPsat(T, Psat1, Psat2)
'[Psat1,Psat2] = CalcPsat(T);
err = 1#
x1_old = xdew1 ' Storing values of x1 ,y1 for control of method
y1_old = ybub1
iter = 0 ' Number of iteration
tol = 0.00000001 ' Tolerance
While err > tol
iter = iter + 1
K1 = gam1 * Psat1 / phi1 / P
K2 = gam2 * Psat2 / phi2 / P
K1 = 1 / (K1 - 1) ' Just for simplifying the form of function
K2 = 1 / (K2 - 1)
' VF = fzero(@(VF) z1/(k1+VF) + z2/(k2+VF) , VF); %Find the zeros of F(VF)
VF = -(K1 * z2 + K2 * z1) / (z1 + z2)
' Correct the value of VF
'-------------------------
If (VF < 0) Then
VF = 0.0000001
Else
If (VF > 1) Then
VF = 0.9999999
'Calculate new value of xi , yi
'------------------------------
x1 = z1 / (1 + VF * (K1 - 1))
y1 = K1 * x1
' x2 = z2/(1+VF*(K2-1))
' y2 = K2*x2
err = Abs(x1_old - x1) + Abs(y1_old - y1)
'err = err + abs(x2_old-x2) + abs(y2_old-y2);
'Updating values of Gamma_i , Phi_i
'----------------------------------
Call Gamma_Wat(x1, T, gam1)
'gam1 = Gamma_Wat(x1,T);
Call Gamma_Vin(1 - x1, T, gam2)
'gam2 = Gamma_Vin(1-x1,T);
Call CalcPhi(T, P, y1, phi1, phi2)
'[Phi1,Phi2] = CalcPhi(T,P,y1);
x1_old = x1
y1_old = y1
Wend
End If
'----------------------------------
'----------------------------------
'----------------------------------
' Printing The Results
Text1(0).Text = "two phase"
Text1(1).Text = Str(iter)
'disp (iter)
Text1(2).Text = Str(x1)
'disp('Water (mol%) In Liquid Phase: ')
'disp (x1)
Text1(3).Text = Str(y1)
'disp('Water (mol%) In Vapor Phase: ')
'disp (y1)
Text1(4).Text = Str(1 - x1)
'disp('Vinyl-Acetate (mol%) In Liquid Phase: ')
'disp (1 - x1)
Text1(5).Text = Str(1 - y1)
'disp('Vinyl-Acetate (mol%) In Vapor Phase: ')
'disp (1 - y1)
Text1(6).Text = Str(VF)
End Sub
Private Sub Text1_Change(Index As Integer)
End Sub
ماژول :
Public i As Integer
Public z1, z2, y1, y2, x1, x2, T, ps1, ps2, Pb, Pn, gamw, gamv, g1_old, g2_old, phi1, phi2 As Double
Public tol, errp, errg As Double
Public Sub CalcPsat(T, ps1, ps2)
'function [ps1,ps2] = CalcPsat(T)
A1 = 11.9647
B1 = 3984.93
C1 = -39.734
ps1 = Exp(A1 - B1 / (C1 + T)) 'bar
A2 = 7.51868
B2 = 1452.058
C2 = -32.562
ps2 = 10 ^ (A2 - B2 / (C2 + T)) 'mmHg
ps2 = ps2 * (1.0133 / 760) 'bar
End Sub
Public Sub Gamma_Wat(x1, T, gamw)
'function gamw= Gamma_Wat(x1,T)
' Calculate the activity of Water by NRTL
' ---------------------------------------
' Data from Aspen Library
aij = 0#
aji = 0#
bij = 1364.6
bji = 415.7
cij = 0.2
t12 = aij + bij / T 'Water
t21 = aji + bji / T 'Vinylacetate
a = cij
x2 = 1 - x1
G12 = Exp(-a * t12)
G21 = Exp(-a * t21)
gamw = x2 * x2 * (t21 * (G21 / (x1 + x2 * G21)) ^ 2 + t12 * G12 / (x2 + x1 * G12) ^ 2)
' gam2 = x1*x1 *( t12* (G12/(x2+x1*G12))^2 + t21*G21/(x1+x2*G21)^2 );
gamw = Exp(gamw)
End Sub
Public Sub Gamma_Vin(x2, T, gamv)
'function gamv= Gamma_Vin(x2,T)
' Calculate the activity of VinylAcetate by NRTL
' -----------------------------------
' Data from Aspen Library
aij = 0#
aji = 0#
bij = 1364.6
bji = 415.7
cij = 0.2
t12 = aij + bij / T 'Water
t21 = aji + bji / T 'Vinylacetate
a = cij
x1 = 1 - x2
G12 = Exp(-a * t12)
G21 = Exp(-a * t21)
'gamw= x2*x2 *( t21* (G21/(x1+x2*G21))^2 + t12*G12/(x2+x1*G12)^2 )
gamv = x1 * x1 * (t12 * (G12 / (x2 + x1 * G12)) ^ 2 + t21 * G21 / (x1 + x2 * G21) ^ 2)
gamv = Exp(gamv)
End Sub
Public Sub CalcPhi(T, P, y1, phi1, phi2)
'function [phi1,phi2]= CalcPhi(T,P,y1)
' Calculate the Fugacity Coefficient by Virial Equation(Two terms)
' ----------------------------------------------------------------
Dim MW, OMEGA, PC, Pc12, TC, Tc12, tr, VC, Vc12, ZC, Zc12, R, B0, B1, B12, Bii, B11, B22, d12 As Double
Dim k12 As Integer
' 1:Water; 2:Vinyl
MW = [18.01528, 86.09044] ' MW
OMEGA = [0.344861, 0.351307] ' OMEGA
PC = [220.55, 39.58] ' PC(bar)
TC = [647.13, 519.13] ' TC(K)
' TC=[373.98,245.98] ' TC(C)
VC = [0.0559478, 0.27] ' VC(m^3/kmol)
ZC = [0.229, 0.248] ' ZC
R = 8.3145 ' KJ/KmolK
' -------------------------------
' -------------------------------
' Calculating Pure Component properties
' -------------------------------------
tr = T / TC
' Pitzer Correlations for the Second Virial Coefficient
B0 = 0.083 - 0.422 / tr ^ 1.6
B1 = 0.139 - 0.172 / tr ^ 4.2
Bii = (TC * R / PC) * (B0 + OMEGA * B1)
B11 = Bii(1)
B22 = Bii(2)
' -------------------------------
' -------------------------------
'Calculating Mixture properties
' -------------------------------
k12 = 0# ' Emperical interaction Parameter ( No Data available so: set to zero)
Tc12 = (TC(1) * TC(2)) ^ 0.5 * (1 - k12)
Omega12 = 0.5 * Sum(OMEGA)
Zc12 = Sum(ZC) / 2
Vc12 = (0.5 * Sum(VC ^ (1 / 3))) ^ 3
Pc12 = Zc12 * R * Tc12 / Vc12
tr = T / Tc12
' Pitzer Correlations for the Second Virial Coefficient
' -----------------------------------------------------
B0 = 0.083 - 0.422 / tr ^ 1.6
B1 = 0.139 - 0.172 / tr ^ 4.2
B12 = (Tc12 * R / Pc12) * (B0 + Omega12 * B1)
d12 = 2 * B12 - B11 - B22
' -------------------------------
' -------------------------------
y2 = 1 - y1
LnPhi1 = P / (R * T) * (B11 + y2 * y2 * d12)
LnPhi2 = P / (R * T) * (B22 + y1 * y1 * d12)
phi1 = Exp(LnPhi1)
phi2 = Exp(LnPhi2)
End Sub
Public Sub Bublp(z1, T, Pb, y1, y2)
' function [pb,y1,y2] = BublP(z1,T)
' Calculate the Bubble point Pressure and xi at T temperature
z2 = 1 - z1
' Calculation of Saturated Pressure of each component at T (K)
' ?????? [Psat1,Psat2] = CalcPsat(T); 'bar
Call CalcPsat(T, Psat1, Psat2)
' -----------------------------------
' -----------------------------------
'Initial Guess
' -----------------------------------
phi1 = 1
phi2 = 1
x1 = z1
x2 = z2
'gam1 = Gamma_Wat(x1, T)
Call Gamma_Wat(x1, T, gam1)
'gam2 = Gamma_Vin(x2, T)
Call Gamma_Vin(x2, T, gam2)
' -----------------------------------
' -----------------------------------
' Initializing Pressure
' ---------------------
Pn = x1 * gam1 * Psat1 / phi1 + x2 * gam2 * Psat2 / phi2
Pb = Pn ' Store Pressure to calculate error
tol = 0.00001 ' Tolerance
errp = 1
While (errp > tol)
'Update yi and Fugacity Coefficients
' -----------------------------------
y1 = x1 * gam1 * Psat1 / phi1 / Pn
y2 = x2 * gam2 * Psat2 / phi2 / Pn
'??? [phi1,phi2]= CalcPhi(T,Pn,y1);
Call CalcPhi(T, Pn, y1, phi1, phi2)
' Update Pressure
' -----------------------------------
Pn = x1 * gam1 * Psat1 / phi1 + x2 * gam2 * Psat2 / phi2
errp = Abs(Pn - Pb)
Pb = Pn
Wend
End Sub
Public Sub DewP(z1, T, Pd, x1, x2)
'function [pd,x1,x2] = DewP(z1,T)
' Calculate the Dew point Pressure and xi at T temperature
z2 = 1 - z1
'Calculation of Saturated Pressure of each component at T (K)
' -----------------------------------------------------------
'[Psat1,Psat2] = CalcPsat(T); %bar
Call CalcPsat(T, Psat1, Psat2)
' -----------------------------------
'Initial Guess
phi1 = 1
phi2 = 1
gam1 = 1
gam2 = 1
y1 = z1
y2 = z2
Pn = 1 / (y1 * phi1 / gam1 / Psat1 + y2 * phi2 / gam2 / Psat2)
x1 = y1 * phi1 * Pn / Psat1 / gam1
x2 = y2 * phi2 * Pn / Psat2 / gam2
'???? gam1 = Gamma_Wat(x1,T);
Call Gamma_Wat(x1, T, gam1)
'???? gam2 = Gamma_Vin(x2,T);
Call Gamma_Vin(x2, T, gam2)
g1_old = gam1 ' Store Gamma to calculate error
g2_old = gam2
Pd = Pn ' Store Pressure to calculate error
tol = 0.00000001
errp = 1 ' error for Calculating Pressure
errg = 1 ' error for Calculating Gamma
While (errp > tol)
'[phi1,phi2]= CalcPhi(T,Pn,y1);
Call CalcPhi(T, Pn, y1, phi1, phi2)
While (errg > tol)
' Calulate xi
'------------
x1 = y1 * phi1 * Pn / Psat1 / gam1
x2 = y2 * phi2 * Pn / Psat2 / gam2
'--- Normilizing xi ---
x12 = x1 + x2
x1 = x1 / x12
x2 = x2 / x12
' Update Gamma
'----------------------
'???? gam1 = Gamma_Wat(x1,T);
Call Gamma_Wat(x1, T, gam1)
'???? gam2 = Gamma_Vin(x2,T);
Call Gamma_Vin(x2, T, gam2)
' Calculate Error
'----------------
errg = Abs(gam1 - g1_old) + Abs(gam2 - g2_old)
g1_old = gam1
g2_old = gam2
Wend
' Update Pressure and Error
'--------------------------
Pn = 1# / (y1 * phi1 / gam1 / Psat1 + y2 * phi2 / gam2 / Psat2)
errp = Abs(Pn - Pd)
Pd = Pn
Wend
End Sub
من یک دانشجوی کارشناسی ارشد مهندسی شیمی هستم.پروژهای که محاسبات فلش برای یک سیستم دو جزئی انجام می دهد را با VB6 نوشتم که همه Subroutine ها در یک ماژول تعریف کردم. نمی دونم چرا در برنامه ای که در فرم اصلی برنامه نوشتم به حلقه while ایراد میگیره؟
نا آخر روز شنبه هم بیشتر فرصت ندارم!کسی میتونه کمکم کنه؟
خطایی که میده:Wend whitout While
در ضمن من چون قبلا برنامه رو با متلب نوشتم،مدل متلب نوشتن رو در برنامه VB پاک نکردم و بصورت توضیحات آوردم.
برنامه اصلی :
Public z1, z2, m1, m2, err, VF, Psat1, Psat2, Pdew, Pbub, Pfac, xdew1, xdew2, x1_old, ybub1, ybub2, y1_old, gam1, gam2, gam_dew1, gam_dew2, gam_bub1, gam_bub2, phi1, phi2, phi_dew1, phi_dew2, phi_bub1, phi_bub2, K1, K2, tol As Double
Public T, V1, V2, Rou1, Rou2 As Double
Public P, iter As Integer
Private Sub Flash_Bottom_Click()
' 1:Water 2:Vinylacetate
'----------------------------------
'----------------------------------
' Input Data
T = 80 + 273.15 'K
P = 1 'bar
V1 = 3 'Lit
V2 = 1.1 'Lit
'----------------------------------
'----------------------------------
' Calculating Mol Fraction of the feed (zi)
'------------------------------------------
' Water :1
V1 = V1 / 1000 'm^3
Rou1 = 971.8 'Kg/m^3
'Mw1 = 18.01528; %Kg/Kmol
m1 = V1 * Rou1 'Kg
' Vinyacetate :2
V2 = V2 / 1000 'm^3
Rou2 = 934# 'kg/m^3
'Mw2 = 86.09044; % Kg/Kmol
m2 = V2 * Rou2 'Kg
z1 = m1 / (m1 + m2)
z2 = 1 - z1
'----------------------------------
'----------------------------------
' Checking the State of mixture
'----------------------------------
Call DewP(z1, T, Pdew, xdew1, xdew2)
'[Pdew,xdew1,xdew2] = DewP(z1,T); % Dew Point Pressure an xi
Call Bublp(z1, T, Pbub, ybub1, ybub2)
'[Pbub,ybub1,ybub2] = BublP(z1,T); % Bubble Point Pressure an yi
If P > Pbub Then
Text1(0).Text = "iquid Phase"
Else
If P < Pdew Then
Text1(0).Text = "Vapor Phase"
Else
'----------------------------------
'----------------------------------
'Performing VLE Flash Calculation
'----------------------------------
' Calculating properties at Dew and bubble Point
Call Gamma_Wat(xdew1, T, gam_dew1)
'gam_dew1 = Gamma_Wat(xdew1,T);
Call Gamma_Vin(xdew2, T, gam_dew2)
'gam_dew2 = Gamma_Vin(xdew2,T);
Call Gamma_Wat(ybub1, T, gam_bub1)
'gam_bub1 = Gamma_Wat(ybub1,T);
Call Gamma_Vin(ybub2, T, gam_bub2)
'gam_bub2 = Gamma_Vin(ybub2,T);
Call CalcPhi(T, Pdew, xdew1, phi_dew1, phi_dew2)
'[phi_dew1,phi_dew2] = CalcPhi(T,Pdew,xdew1);
Call CalcPhi(T, Pbub, ybub1, phi_bub1, phi_bub2)
'[phi_bub1,phi_bub2] = CalcPhi(T,Pbub,ybub1);
'Initial Guess by Intrpolation
'initializing by means of Lever Rule
Pfac = (P - Pdew) / (Pbub - Pdew)
gam1 = Pfac * (gam_bub1 - gam_dew1) + gam_dew1
gam2 = Pfac * (gam_bub2 - gam_dew2) + gam_dew2
phi1 = Pfac * (phi_bub1 - phi_dew1) + phi_dew1
phi2 = Pfac * (phi_bub2 - phi_dew2) + phi_dew2
VF = (Pbub - P) / (Pbub - Pdew) 'Initial Guess for V/F
'--------------------------------------------
Call CalcPsat(T, Psat1, Psat2)
'[Psat1,Psat2] = CalcPsat(T);
err = 1#
x1_old = xdew1 ' Storing values of x1 ,y1 for control of method
y1_old = ybub1
iter = 0 ' Number of iteration
tol = 0.00000001 ' Tolerance
While err > tol
iter = iter + 1
K1 = gam1 * Psat1 / phi1 / P
K2 = gam2 * Psat2 / phi2 / P
K1 = 1 / (K1 - 1) ' Just for simplifying the form of function
K2 = 1 / (K2 - 1)
' VF = fzero(@(VF) z1/(k1+VF) + z2/(k2+VF) , VF); %Find the zeros of F(VF)
VF = -(K1 * z2 + K2 * z1) / (z1 + z2)
' Correct the value of VF
'-------------------------
If (VF < 0) Then
VF = 0.0000001
Else
If (VF > 1) Then
VF = 0.9999999
'Calculate new value of xi , yi
'------------------------------
x1 = z1 / (1 + VF * (K1 - 1))
y1 = K1 * x1
' x2 = z2/(1+VF*(K2-1))
' y2 = K2*x2
err = Abs(x1_old - x1) + Abs(y1_old - y1)
'err = err + abs(x2_old-x2) + abs(y2_old-y2);
'Updating values of Gamma_i , Phi_i
'----------------------------------
Call Gamma_Wat(x1, T, gam1)
'gam1 = Gamma_Wat(x1,T);
Call Gamma_Vin(1 - x1, T, gam2)
'gam2 = Gamma_Vin(1-x1,T);
Call CalcPhi(T, P, y1, phi1, phi2)
'[Phi1,Phi2] = CalcPhi(T,P,y1);
x1_old = x1
y1_old = y1
Wend
End If
'----------------------------------
'----------------------------------
'----------------------------------
' Printing The Results
Text1(0).Text = "two phase"
Text1(1).Text = Str(iter)
'disp (iter)
Text1(2).Text = Str(x1)
'disp('Water (mol%) In Liquid Phase: ')
'disp (x1)
Text1(3).Text = Str(y1)
'disp('Water (mol%) In Vapor Phase: ')
'disp (y1)
Text1(4).Text = Str(1 - x1)
'disp('Vinyl-Acetate (mol%) In Liquid Phase: ')
'disp (1 - x1)
Text1(5).Text = Str(1 - y1)
'disp('Vinyl-Acetate (mol%) In Vapor Phase: ')
'disp (1 - y1)
Text1(6).Text = Str(VF)
End Sub
Private Sub Text1_Change(Index As Integer)
End Sub
ماژول :
Public i As Integer
Public z1, z2, y1, y2, x1, x2, T, ps1, ps2, Pb, Pn, gamw, gamv, g1_old, g2_old, phi1, phi2 As Double
Public tol, errp, errg As Double
Public Sub CalcPsat(T, ps1, ps2)
'function [ps1,ps2] = CalcPsat(T)
A1 = 11.9647
B1 = 3984.93
C1 = -39.734
ps1 = Exp(A1 - B1 / (C1 + T)) 'bar
A2 = 7.51868
B2 = 1452.058
C2 = -32.562
ps2 = 10 ^ (A2 - B2 / (C2 + T)) 'mmHg
ps2 = ps2 * (1.0133 / 760) 'bar
End Sub
Public Sub Gamma_Wat(x1, T, gamw)
'function gamw= Gamma_Wat(x1,T)
' Calculate the activity of Water by NRTL
' ---------------------------------------
' Data from Aspen Library
aij = 0#
aji = 0#
bij = 1364.6
bji = 415.7
cij = 0.2
t12 = aij + bij / T 'Water
t21 = aji + bji / T 'Vinylacetate
a = cij
x2 = 1 - x1
G12 = Exp(-a * t12)
G21 = Exp(-a * t21)
gamw = x2 * x2 * (t21 * (G21 / (x1 + x2 * G21)) ^ 2 + t12 * G12 / (x2 + x1 * G12) ^ 2)
' gam2 = x1*x1 *( t12* (G12/(x2+x1*G12))^2 + t21*G21/(x1+x2*G21)^2 );
gamw = Exp(gamw)
End Sub
Public Sub Gamma_Vin(x2, T, gamv)
'function gamv= Gamma_Vin(x2,T)
' Calculate the activity of VinylAcetate by NRTL
' -----------------------------------
' Data from Aspen Library
aij = 0#
aji = 0#
bij = 1364.6
bji = 415.7
cij = 0.2
t12 = aij + bij / T 'Water
t21 = aji + bji / T 'Vinylacetate
a = cij
x1 = 1 - x2
G12 = Exp(-a * t12)
G21 = Exp(-a * t21)
'gamw= x2*x2 *( t21* (G21/(x1+x2*G21))^2 + t12*G12/(x2+x1*G12)^2 )
gamv = x1 * x1 * (t12 * (G12 / (x2 + x1 * G12)) ^ 2 + t21 * G21 / (x1 + x2 * G21) ^ 2)
gamv = Exp(gamv)
End Sub
Public Sub CalcPhi(T, P, y1, phi1, phi2)
'function [phi1,phi2]= CalcPhi(T,P,y1)
' Calculate the Fugacity Coefficient by Virial Equation(Two terms)
' ----------------------------------------------------------------
Dim MW, OMEGA, PC, Pc12, TC, Tc12, tr, VC, Vc12, ZC, Zc12, R, B0, B1, B12, Bii, B11, B22, d12 As Double
Dim k12 As Integer
' 1:Water; 2:Vinyl
MW = [18.01528, 86.09044] ' MW
OMEGA = [0.344861, 0.351307] ' OMEGA
PC = [220.55, 39.58] ' PC(bar)
TC = [647.13, 519.13] ' TC(K)
' TC=[373.98,245.98] ' TC(C)
VC = [0.0559478, 0.27] ' VC(m^3/kmol)
ZC = [0.229, 0.248] ' ZC
R = 8.3145 ' KJ/KmolK
' -------------------------------
' -------------------------------
' Calculating Pure Component properties
' -------------------------------------
tr = T / TC
' Pitzer Correlations for the Second Virial Coefficient
B0 = 0.083 - 0.422 / tr ^ 1.6
B1 = 0.139 - 0.172 / tr ^ 4.2
Bii = (TC * R / PC) * (B0 + OMEGA * B1)
B11 = Bii(1)
B22 = Bii(2)
' -------------------------------
' -------------------------------
'Calculating Mixture properties
' -------------------------------
k12 = 0# ' Emperical interaction Parameter ( No Data available so: set to zero)
Tc12 = (TC(1) * TC(2)) ^ 0.5 * (1 - k12)
Omega12 = 0.5 * Sum(OMEGA)
Zc12 = Sum(ZC) / 2
Vc12 = (0.5 * Sum(VC ^ (1 / 3))) ^ 3
Pc12 = Zc12 * R * Tc12 / Vc12
tr = T / Tc12
' Pitzer Correlations for the Second Virial Coefficient
' -----------------------------------------------------
B0 = 0.083 - 0.422 / tr ^ 1.6
B1 = 0.139 - 0.172 / tr ^ 4.2
B12 = (Tc12 * R / Pc12) * (B0 + Omega12 * B1)
d12 = 2 * B12 - B11 - B22
' -------------------------------
' -------------------------------
y2 = 1 - y1
LnPhi1 = P / (R * T) * (B11 + y2 * y2 * d12)
LnPhi2 = P / (R * T) * (B22 + y1 * y1 * d12)
phi1 = Exp(LnPhi1)
phi2 = Exp(LnPhi2)
End Sub
Public Sub Bublp(z1, T, Pb, y1, y2)
' function [pb,y1,y2] = BublP(z1,T)
' Calculate the Bubble point Pressure and xi at T temperature
z2 = 1 - z1
' Calculation of Saturated Pressure of each component at T (K)
' ?????? [Psat1,Psat2] = CalcPsat(T); 'bar
Call CalcPsat(T, Psat1, Psat2)
' -----------------------------------
' -----------------------------------
'Initial Guess
' -----------------------------------
phi1 = 1
phi2 = 1
x1 = z1
x2 = z2
'gam1 = Gamma_Wat(x1, T)
Call Gamma_Wat(x1, T, gam1)
'gam2 = Gamma_Vin(x2, T)
Call Gamma_Vin(x2, T, gam2)
' -----------------------------------
' -----------------------------------
' Initializing Pressure
' ---------------------
Pn = x1 * gam1 * Psat1 / phi1 + x2 * gam2 * Psat2 / phi2
Pb = Pn ' Store Pressure to calculate error
tol = 0.00001 ' Tolerance
errp = 1
While (errp > tol)
'Update yi and Fugacity Coefficients
' -----------------------------------
y1 = x1 * gam1 * Psat1 / phi1 / Pn
y2 = x2 * gam2 * Psat2 / phi2 / Pn
'??? [phi1,phi2]= CalcPhi(T,Pn,y1);
Call CalcPhi(T, Pn, y1, phi1, phi2)
' Update Pressure
' -----------------------------------
Pn = x1 * gam1 * Psat1 / phi1 + x2 * gam2 * Psat2 / phi2
errp = Abs(Pn - Pb)
Pb = Pn
Wend
End Sub
Public Sub DewP(z1, T, Pd, x1, x2)
'function [pd,x1,x2] = DewP(z1,T)
' Calculate the Dew point Pressure and xi at T temperature
z2 = 1 - z1
'Calculation of Saturated Pressure of each component at T (K)
' -----------------------------------------------------------
'[Psat1,Psat2] = CalcPsat(T); %bar
Call CalcPsat(T, Psat1, Psat2)
' -----------------------------------
'Initial Guess
phi1 = 1
phi2 = 1
gam1 = 1
gam2 = 1
y1 = z1
y2 = z2
Pn = 1 / (y1 * phi1 / gam1 / Psat1 + y2 * phi2 / gam2 / Psat2)
x1 = y1 * phi1 * Pn / Psat1 / gam1
x2 = y2 * phi2 * Pn / Psat2 / gam2
'???? gam1 = Gamma_Wat(x1,T);
Call Gamma_Wat(x1, T, gam1)
'???? gam2 = Gamma_Vin(x2,T);
Call Gamma_Vin(x2, T, gam2)
g1_old = gam1 ' Store Gamma to calculate error
g2_old = gam2
Pd = Pn ' Store Pressure to calculate error
tol = 0.00000001
errp = 1 ' error for Calculating Pressure
errg = 1 ' error for Calculating Gamma
While (errp > tol)
'[phi1,phi2]= CalcPhi(T,Pn,y1);
Call CalcPhi(T, Pn, y1, phi1, phi2)
While (errg > tol)
' Calulate xi
'------------
x1 = y1 * phi1 * Pn / Psat1 / gam1
x2 = y2 * phi2 * Pn / Psat2 / gam2
'--- Normilizing xi ---
x12 = x1 + x2
x1 = x1 / x12
x2 = x2 / x12
' Update Gamma
'----------------------
'???? gam1 = Gamma_Wat(x1,T);
Call Gamma_Wat(x1, T, gam1)
'???? gam2 = Gamma_Vin(x2,T);
Call Gamma_Vin(x2, T, gam2)
' Calculate Error
'----------------
errg = Abs(gam1 - g1_old) + Abs(gam2 - g2_old)
g1_old = gam1
g2_old = gam2
Wend
' Update Pressure and Error
'--------------------------
Pn = 1# / (y1 * phi1 / gam1 / Psat1 + y2 * phi2 / gam2 / Psat2)
errp = Abs(Pn - Pd)
Pd = Pn
Wend
End Sub