GOAT/tpdss.m at main · smithh2710/GOAT · GitHub

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
%% Tangent Plane Distance analysis by using successive substitution
% -------------------------------------------------------------------------
% The Definition of Variables.
% comp_overall : Overall composition
% press        : Pressure
% temp         : Temperature
% pressc       : Critical pressure
% tempc        : Critical temperature
% acentric     : Acentric factor
% BIP          : Binary interaction parameters
% eos_type     : 'PR' (default) or 'SRK'
% twophase     : true if phase split occurs, false otherwise
% -------------------------------------------------------------------------

function twophase = tpdss(comp_overall, press, temp, pressc, tempc, acentric, BIP, eos_type)

% Default to PR if eos_type not specified
if nargin < 8 || isempty(eos_type)
    eos_type = 'PR';
end

% The initial estimate of equilibrium constant
K = wilsoneq(press, temp, pressc, tempc, acentric);

% Vapor-like phase
comp_vap = compvapest(K, comp_overall);
vaporphase = checkphasesplit(comp_vap, comp_overall, press, temp, pressc, tempc, acentric, BIP, eos_type);

% Liquid-like phase
comp_liq = compliqest(K, comp_overall);
liquidphase = checkphasesplit(comp_liq, comp_overall, press, temp, pressc, tempc, acentric, BIP, eos_type);

if (vaporphase == true) || (liquidphase == true)
    twophase = true;
else
    twophase = false;
end

end

%% Check if phase split occurs
function phasesplit = checkphasesplit(comp_check, comp_overall, press, temp, pressc, tempc, acentric, BIP, eos_type)

[fugcoef_overall, ~] = fugacitycoef_multicomp(comp_overall, press, temp, pressc, tempc, acentric, BIP, eos_type);

ncomp = size(comp_overall, 1);
x = comp_check;
tol = 1e-6;
eps = 1;
maxiter = 20;

for iter = 1:maxiter

    xsum = sum(x);
    x = x/xsum;

    if max(abs(x - comp_overall)) < tol
        phasesplit = false;
        return;
    end

    if eps < tol
        break;
    end

    [fugcoef_x, ~] = fugacitycoef_multicomp(x, press, temp, pressc, tempc, acentric, BIP, eos_type);
    eps = calcrsd(x, fugcoef_x, comp_overall, fugcoef_overall);
    x = updatex(comp_overall, fugcoef_overall, fugcoef_x);

end

if iter >= maxiter
    fprintf('The iteration in checkphasesplit() did not converge. One phase is assumed.\n');
end

if xsum > 1
    phasesplit = true;
else
    phasesplit = false;
end

end

%% Estimate vapor composition from K-values
function comp_vap = compvapest(K, comp_overall)

ncomp = size(comp_overall,1);

comp_vap = zeros(ncomp,1);
for i = 1:ncomp
    comp_vap(i) = K(i)*comp_overall(i);
end

end

%% Estimate liquid composition from K-values
function comp_liq = compliqest(K, comp_overall)

ncomp = size(comp_overall,1);

comp_liq = zeros(ncomp,1);
for i = 1:ncomp
    comp_liq(i) = comp_overall(i)/K(i);
end

end

%% Calculate residual
function rsd = calcrsd(x, fugcoef_x, comp_overall, fugcoef_overall)

ncomp = size(x,1);

rsd = 0;
for i = 1:ncomp
    rsd = rsd + ( log(x(i)/comp_overall(i)*fugcoef_x(i)/fugcoef_overall(i)) )^2;
end

end

%% Update trial composition
function x = updatex(comp_overall, fugcoef_overall, fugcoef_x)

ncomp = size(comp_overall,1);

x = zeros(ncomp,1);
for i = 1:ncomp
    x(i) = comp_overall(i)*fugcoef_overall(i)/fugcoef_x(i);
end

end