Skip to content

Commit 84a5759

Browse files
simulkadesimulkade
committed
spent some claude code tokens on new features
1 parent 0ef4292 commit 84a5759

15 files changed

Lines changed: 1582 additions & 671 deletions

‎EOS/PR78EOS.m‎

Lines changed: 158 additions & 182 deletions
Original file line numberDiff line numberDiff line change
@@ -1,182 +1,158 @@
1-
function [liquid_z, vapor_z, fugacity, HR] = PR78EOS(mixture, thermo)
2-
% PR78EOS Peng-Robinson EOS with 1978 alpha-function correction.
3-
%
4-
% [liquid_z, vapor_z, fugacity, HR] = PR78EOS(mixture, thermo)
5-
%
6-
% Identical to PREOS except the alpha function is modified for components
7-
% with acentric factor > 0.491 (heavier hydrocarbons, polar molecules):
8-
% ω ≤ 0.491: m = 0.37646 + 1.54226·ω − 0.26992·ω² (same as PREOS)
9-
% ω > 0.491: m = 0.379642 + 1.48503·ω − 0.164423·ω² + 0.016666·ω³
10-
%
11-
% PARAMETERS:
12-
% mixture - Mixture object; relevant fields:
13-
% .temperature, .pressure, .mole_fraction, .components, .bip
14-
% thermo - ThermoModel object; relevant fields:
15-
% .mixingrule, .activity_model, .phase, .fugacity_switch
16-
%
17-
% RETURNS:
18-
% liquid_z - liquid compressibility factor
19-
% vapor_z - vapor compressibility factor
20-
% fugacity - [1 x N] fugacities [Pa] (zero if fugacity_switch == 0)
21-
% HR - residual molar enthalpy [J/mol]
22-
%
23-
% SEE ALSO: PREOS, SRKEOS, ThermoModel
24-
25-
%{
26-
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
27-
All rights reserved.
28-
29-
Redistribution and use in source and binary forms, with or
30-
without modification, are permitted provided that the following
31-
conditions are met:
32-
33-
* Redistributions of source code must retain the above copyright notice,
34-
this list of conditions and the following disclaimer.
35-
* Redistributions in binary form must reproduce the above
36-
copyright notice, this list of conditions and the following
37-
disclaimer in the documentation and/or other materials provided
38-
with the distribution.
39-
40-
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
41-
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
42-
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
43-
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
44-
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
45-
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
46-
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
47-
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
48-
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
49-
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
50-
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
51-
%}
52-
53-
54-
mixing_rule_num = thermo.mixingrule;
55-
activityfun = thermo.activity_model;
56-
phase1 = thermo.phase;
57-
fug_need = thermo.fugacity_switch;
58-
critical_pres = [mixture.components.Pc]; %[Pa]
59-
critical_temp = [mixture.components.Tc]; %[K]
60-
acentric_fact = [mixture.components.acentric_factor]; %[-]
61-
BIP = mixture.bip;
62-
x = mixture.mole_fraction;
63-
p = mixture.pressure; %[Pa]
64-
T = mixture.temperature;
65-
N = length(critical_temp);
66-
fugacity = zeros(1,N);
67-
R=8.314;
68-
s1 = 0.623225; %Huron Vidal
69-
70-
bi = 0.0777960739*R*critical_temp./critical_pres;
71-
aci=0.457235529*(R*critical_temp).^2 ./critical_pres;
72-
mi = (0.37646+(1.54226-0.26992*acentric_fact).*acentric_fact).*(acentric_fact<=0.491)+ ...
73-
(0.379642+1.48503*acentric_fact-0.164423*acentric_fact.^2+0.016666*acentric_fact.^3).*(acentric_fact>0.491);
74-
Tr = T./critical_temp;
75-
alfai = 1+mi.*(1-sqrt(Tr)); %//alfai=ai^0.5
76-
alfa = alfai.^2.0; %//alfa=ai
77-
ai = aci .* alfa;
78-
%Q is the parameter for MHV1 and MHV2 mixing rule
79-
%it depends on the EOS
80-
Q = (mixing_rule_num==3)*[-0.53 0]+(mixing_rule_num==4)*[-0.4347 -0.003654];
81-
[a, b] = mixing_rule(mixture, thermo, ai, bi, s1, Q);
82-
83-
A_coef=a*p/(R*T)^2;
84-
B_coef=b*p/(R*T);
85-
% poly_coef(4)=-B_coef*(A_coef-B_coef*(1+B_coef));
86-
% poly_coef(3)=A_coef-B_coef*(2.0+3*B_coef);
87-
% poly_coef(2)=-1+B_coef;
88-
% poly_coef(1)=1;
89-
poly_coef = [1 -1+B_coef A_coef-B_coef*(2.0+3*B_coef) ...
90-
-B_coef*(A_coef-B_coef*(1+B_coef))];
91-
92-
93-
z_root = roots(poly_coef);
94-
if (sum(imag(z_root)~=0)==0)
95-
liquid_z = min(z_root);
96-
vapor_z = max(z_root);
97-
else
98-
liquid_z = z_root(imag(z_root)==0);
99-
vapor_z = liquid_z;
100-
end
101-
% %{now we should calculate vapor and liquid
102-
% %compressibility factor by the following method :}
103-
% real_part = real(z_root);
104-
% imag_part = imag(z_root);
105-
% % findmaxmin(real_part,img_part,zv,zl);
106-
% max_no=real_part(3);
107-
% min_no=real_part(3);
108-
% for i=1:3
109-
% if ((real_part(i)>max_no) && (real_part(i)>0) && (imag_part(i)==0))
110-
% max_no=real_part(i);
111-
% end
112-
% if ((real_part(i)<min_no) && (real_part(i)>0) && (imag_part(i)==0))
113-
% min_no=real_part(i);
114-
% end
115-
% end
116-
% liquid_z = min_no;
117-
% vapor_z = max_no;
118-
if (phase1==1) %then //phase 1 is liquid
119-
zz=liquid_z;
120-
elseif (phase1==2) %then //phase 2 is vapor
121-
zz=vapor_z;
122-
end
123-
124-
if (fug_need==1)
125-
if (mixing_rule_num==1) %simple mixing rule
126-
part1=bi/b*(zz-1)-log(zz-b*p/(R*T));
127-
part2=x*(sqrt(ai'*ai).*(1-[BIP.EOScons]-[BIP.EOStdep]*T))';
128-
part3=A_coef/(2.828*B_coef)*(bi/b-2/a*part2) ...
129-
*log((zz+2.414*b*p/(R*T))/(zz-0.414*b*p/(R*T)));
130-
fugacity=exp(part1+part3);
131-
elseif (mixing_rule_num==2) %Huron Vidal mixing rule
132-
[~, gama] = activityfun(T, x, mixture.components, BIP);
133-
part1 = bi/b*(zz-1)-log(zz-b*p/(R*T));
134-
part3 = -1/(2*sqrt(2))*(ai./bi/R/T - ...
135-
log(gama)/0.623225)* ...
136-
log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
137-
fugacity = exp(part1+part3);
138-
elseif (mixing_rule_num==3) % MHV1 mixing rule
139-
[~, gama] = activityfun(T, x, mixture.components, BIP);
140-
q1 = -0.53; %Michelsen for PR
141-
logfi = bi/b*(zz-1) - log(zz-B_coef) - 1/(2*sqrt(2))*(ai./(bi*R*T) ...
142-
+ log(gama)/q1+log(b./bi)/q1+(bi/b-1)/q1)*log((zz+(1+sqrt(2))*B_coef)/ ...
143-
(zz+(1-sqrt(2))*B_coef));
144-
fugacity = exp(logfi);
145-
elseif (mixing_rule_num==4) % MHV2 mixing rule
146-
[~, gama] = activityfun(T, x, mixture.components, BIP);
147-
q1 = -0.4347;
148-
q2 = -0.003654;
149-
alphai = ai./(bi*R*T);
150-
alpha = a/(b*R*T);
151-
logfi = bi/b*(zz-1) - log(zz-B_coef) - ...
152-
1/(2*sqrt(2))*(q1*alphai+q2*(alpha^2+alphai.^2)+log(gama) ...
153-
+log(b./bi)+bi/b-1)/(q1+2*q2*alpha)* ...
154-
log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
155-
fugacity = exp(logfi);
156-
end
157-
end
158-
159-
160-
for i=1:N
161-
Abar(i)=0;
162-
for j=1:N
163-
Abar(i)=Abar(i)+sqrt(ai(j))/(R*T)*x(j);
164-
end
165-
Abar(i)=Abar(i)*sqrt(ai(i))/(R*T);
166-
end
167-
part1=0;
168-
for i=1:N
169-
part1=part1+x(i)*Abar(i)*mi(i)*(-1)/sqrt(critical_temp(i))/alfai(i);
170-
end
171-
part1=(part1*T^0.5/(A_coef/p)-1)*R*T*A_coef/(B_coef*2*sqrt(2))*log((zz+B_coef*(1+sqrt(2)))/(zz+B_coef*(1-sqrt(2))));
172-
HR=part1+R*T*(zz-1);
173-
% HR=0;
174-
175-
176-
177-
178-
179-
180-
181-
182-
1+
function [liquid_z, vapor_z, fugacity, HR, props] = PR78EOS(mixture, thermo)
2+
% PR78EOS Peng-Robinson EOS with 1978 alpha-function correction.
3+
%
4+
% [liquid_z, vapor_z, fugacity, HR] = PR78EOS(mixture, thermo)
5+
% [liquid_z, vapor_z, fugacity, HR, props] = PR78EOS(mixture, thermo)
6+
%
7+
% Identical to PREOS except the alpha function is modified for components
8+
% with acentric factor > 0.491 (heavier hydrocarbons, polar molecules):
9+
% ω ≤ 0.491: m = 0.37646 + 1.54226·ω − 0.26992·ω² (same as PREOS)
10+
% ω > 0.491: m = 0.379642 + 1.48503·ω − 0.164423·ω² + 0.016666·ω³
11+
%
12+
% PARAMETERS:
13+
% mixture - Mixture object; relevant fields:
14+
% .temperature, .pressure, .mole_fraction, .components, .bip
15+
% thermo - ThermoModel object; relevant fields:
16+
% .mixingrule, .activity_model, .phase, .fugacity_switch
17+
%
18+
% RETURNS:
19+
% liquid_z - liquid compressibility factor (smallest real root > B)
20+
% vapor_z - vapor compressibility factor (largest real root > B)
21+
% fugacity - [1 x N] fugacity coefficients (zero if fugacity_switch==0)
22+
% HR - residual molar enthalpy [J/mol]
23+
% props - struct with fields HR, SR, GR, VR, Cp_R, Cv_R (optional 5th output)
24+
%
25+
% SEE ALSO: PREOS, SRKEOS, select_z_roots, ThermoModel
26+
27+
%{
28+
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
29+
All rights reserved.
30+
31+
Redistribution and use in source and binary forms, with or
32+
without modification, are permitted provided that the following
33+
conditions are met:
34+
35+
* Redistributions of source code must retain the above copyright notice,
36+
this list of conditions and the following disclaimer.
37+
* Redistributions in binary form must reproduce the above
38+
copyright notice, this list of conditions and the following
39+
disclaimer in the documentation and/or other materials provided
40+
with the distribution.
41+
42+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
43+
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
44+
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
45+
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
46+
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
47+
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
48+
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
49+
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
50+
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
51+
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
52+
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
53+
%}
54+
55+
mixing_rule_num = thermo.mixingrule;
56+
activityfun = thermo.activity_model;
57+
phase1 = thermo.phase;
58+
fug_need = thermo.fugacity_switch;
59+
critical_pres = [mixture.components.Pc]; %[Pa]
60+
critical_temp = [mixture.components.Tc]; %[K]
61+
acentric_fact = [mixture.components.acentric_factor]; %[-]
62+
BIP = mixture.bip;
63+
x = mixture.mole_fraction;
64+
p = mixture.pressure; %[Pa]
65+
T = mixture.temperature;
66+
N = length(critical_temp);
67+
fugacity = zeros(1,N);
68+
R = 8.314;
69+
s1 = 0.623225; % Huron-Vidal constant for PR
70+
71+
bi = 0.0777960739*R*critical_temp./critical_pres;
72+
aci = 0.457235529*(R*critical_temp).^2 ./critical_pres;
73+
% PR78 alpha correction: different m for omega > 0.491
74+
mi = (0.37646+(1.54226-0.26992*acentric_fact).*acentric_fact).*(acentric_fact<=0.491)+ ...
75+
(0.379642+1.48503*acentric_fact-0.164423*acentric_fact.^2+0.016666*acentric_fact.^3).*(acentric_fact>0.491);
76+
Tr = T./critical_temp;
77+
alfai = 1+mi.*(1-sqrt(Tr));
78+
alfa = alfai.^2;
79+
ai = aci .* alfa;
80+
81+
Q = (mixing_rule_num==3)*[-0.53 0]+(mixing_rule_num==4)*[-0.4347 -0.003654];
82+
[a, b] = mixing_rule(mixture, thermo, ai, bi, s1, Q);
83+
84+
A_coef = a*p/(R*T)^2;
85+
B_coef = b*p/(R*T);
86+
87+
poly_coef = [1 -1+B_coef A_coef-B_coef*(2.0+3*B_coef) ...
88+
-B_coef*(A_coef-B_coef*(1+B_coef))];
89+
90+
z_root = roots(poly_coef);
91+
[liquid_z, vapor_z] = select_z_roots(z_root, B_coef);
92+
93+
if (phase1 == 1)
94+
zz = liquid_z;
95+
else
96+
zz = vapor_z;
97+
end
98+
99+
if (fug_need == 1)
100+
if (mixing_rule_num == 1)
101+
part1 = bi/b*(zz-1) - log(zz-b*p/(R*T));
102+
part2 = x*(sqrt(ai'*ai).*(1-[BIP.EOScons]-[BIP.EOStdep]*T))';
103+
part3 = A_coef/(2.828*B_coef)*(bi/b-2/a*part2) ...
104+
*log((zz+2.414*b*p/(R*T))/(zz-0.414*b*p/(R*T)));
105+
fugacity = exp(part1+part3);
106+
elseif (mixing_rule_num == 2)
107+
[~, gama] = activityfun(T, x, mixture.components, BIP);
108+
part1 = bi/b*(zz-1) - log(zz-b*p/(R*T));
109+
part3 = -1/(2*sqrt(2))*(ai./bi/R/T - ...
110+
log(gama)/0.623225)* ...
111+
log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
112+
fugacity = exp(part1+part3);
113+
elseif (mixing_rule_num == 3)
114+
[~, gama] = activityfun(T, x, mixture.components, BIP);
115+
q1 = -0.53;
116+
logfi = bi/b*(zz-1) - log(zz-B_coef) - 1/(2*sqrt(2))*(ai./(bi*R*T) ...
117+
+ log(gama)/q1+log(b./bi)/q1+(bi/b-1)/q1)*log((zz+(1+sqrt(2))*B_coef)/ ...
118+
(zz+(1-sqrt(2))*B_coef));
119+
fugacity = exp(logfi);
120+
elseif (mixing_rule_num == 4)
121+
[~, gama] = activityfun(T, x, mixture.components, BIP);
122+
q1 = -0.4347;
123+
q2 = -0.003654;
124+
alphai = ai./(bi*R*T);
125+
alpha = a/(b*R*T);
126+
logfi = bi/b*(zz-1) - log(zz-B_coef) - ...
127+
1/(2*sqrt(2))*(q1*alphai+q2*(alpha^2+alphai.^2)+log(gama) ...
128+
+log(b./bi)+bi/b-1)/(q1+2*q2*alpha)* ...
129+
log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
130+
fugacity = exp(logfi);
131+
end
132+
end
133+
134+
%% Residual properties (PR78 — same integral as PR, different mi)
135+
dadT = -(sum(x.*sqrt(ai))) .* sum(x.*sqrt(aci).*mi./sqrt(critical_temp)) ./ sqrt(T);
136+
137+
ln_term = log((zz+B_coef*(1+sqrt(2))) / (zz+B_coef*(1-sqrt(2))));
138+
HR = R*T*(zz-1) + (T*dadT - a) / (b*2*sqrt(2)) * ln_term;
139+
140+
if nargout >= 5
141+
SR = R*log(zz - B_coef) + dadT / (b*2*sqrt(2)) * ln_term;
142+
GR = HR - T*SR;
143+
VR = R*T*(zz-1) / p;
144+
145+
dsqrtaidT = -sqrt(aci) .* mi ./ (2*sqrt(critical_temp*T));
146+
d2sqrtaidT2 = sqrt(aci) .* mi ./ (4*sqrt(critical_temp) .* T^(3/2));
147+
d2adT2 = 2*(sum(x.*dsqrtaidT))^2 + 2*sum(x.*sqrt(ai))*sum(x.*d2sqrtaidT2);
148+
149+
Cv_R = -T * d2adT2 / (b*2*sqrt(2)) * ln_term;
150+
151+
V_mol = zz*R*T / p;
152+
den_pr = V_mol^2 + 2*b*V_mol - b^2;
153+
dPdT_V = R/(V_mol-b) - dadT/den_pr;
154+
dPdV_T = -R*T/(V_mol-b)^2 + 2*a*(V_mol+b)/den_pr^2;
155+
Cp_R = Cv_R - T*dPdT_V^2/dPdV_T - R;
156+
157+
props = struct('HR',HR, 'SR',SR, 'GR',GR, 'VR',VR, 'Cp_R',Cp_R, 'Cv_R',Cv_R);
158+
end

0 commit comments

Comments
 (0)