Hi, I keep getting y as a 3x3 matrix, but its supposed to come out as a 3x1 matrix. Could someone please find my mistake. I think Im missing a matrix decimal divide or multipl
Show older comments
% Take home problem 2
% Water (1), Acetone (2), Methanol (3)
R = 83.14; % bar cm^3 mol^-1 K^-1
P = 2; %bars
x = [0.35; 0.25; 0.40];
V = [18.07; 74.05; 40.73]; % cm^3/mol
a = [0 6062.5 1965.9; 1219.5 0 -677.8; 449.6 2441.4 0]; % J/mol
w = [0.345; 0.307; 0.564];
Tc =[647.1; 508.2; 512.6]; % K
Pc = [220.55; 47.01; 80.97]; % bar
Zc = [0.229; 0.233; 0.224];
Vc = [55.9; 209; 118]; % cm^3/mol
% Antoine’s equation constants (with these constants P is in kPa and T in Celsius
A = [16.3872; 14.3145; 16.5785];
B = [3885.70; 2756.22; 3638.27];
C = [230.170; 228.060; 239.500];
N = max(size(w)); % number of components
phi = ones(N);
Tsat = B./(A - log(100*P)) - C + 273.15
T = Tsat'*x
Psat = exp(A-B./(T-273.15+C))/100;
gamma = wilson(x, a, V, R*T);
% Let’s simply choose Pjsat as Psat(1)
s = 0; % initialize summation
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
T = B(1)/(A(1) - log(100*Pjsat)) - C(1) + 273.15; % temperature in Kelvins
Psat = exp(A-B./(T-273.15+C))/100;
y= x.*gamma.*Psat./(phi.*P);
[BB d] = virial2(T, Tc, Pc, Zc, Vc, w); %virial coefficients needed to calculate fugacity coefficients
phi = virialcoeff(BB, d, T, y, P); % fugacity coefficients
gamma = wilson(x, a, V, R*T);
s = 0;
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
T = B(1)/(A(1) - log(100*Pjsat)) - C(1) + 273.15;
er1 = 1; % initialize error to a value large enough that it can’t be satisfied the first iteration
while er1 > 1e-3
Psat = exp(A-B./(T-273.15+C))/100;
y=x.*gamma.*Psat./(phi.*P);
[BB d] = virial2(T, Tc, Pc, Zc, Vc, w); %virial coefficients needed to calculate fugacity coefficients
phi = virialcoeff(BB, d, T, y, P); % fugacity coefficients
gamma = wilson(x, a, V, R*T);
s = 0; % initialize summation
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
Tnew = B(1)/(A(1)-log(100*Pjsat))-C(1)+273.15; % T in K
er1 = abs(1-Tnew/T)
T = Tnew;
end
T
y
2 Comments
the cyclist
on 22 Jul 2022
It's difficult to debug this, without access to the subfunctions (e.g. wilson). I recommend setting breakpoints so that you can see what shape/size your other variables are.
Tony Mei
on 22 Jul 2022
Accepted Answer
% Take home problem 2
% Water (1), Acetone (2), Methanol (3)
R = 83.14; % bar cm^3 mol^-1 K^-1
P = 2; %bars
x = [0.35; 0.25; 0.40] ;
V = [18.07; 74.05; 40.73]; % cm^3/mol
a = [0 6062.5 1965.9; 1219.5 0 -677.8; 449.6 2441.4 0]; % J/mol
w = [0.345; 0.307; 0.564];
Tc =[647.1; 508.2; 512.6]; % K
Pc = [220.55; 47.01; 80.97]; % bar
Zc = [0.229; 0.233; 0.224];
Vc = [55.9; 209; 118]; % cm^3/mol
% Antoine’s equation constants (with these constants P is in kPa and T in Celsius
A = [16.3872; 14.3145; 16.5785];
B = [3885.70; 2756.22; 3638.27];
C = [230.170; 228.060; 239.500];
N = max(size(w)); % number of components
phi = ones(N);
Tsat = B./(A - log(100*P)) - C + 273.15;
T = Tsat'*x;
Psat = exp(A-B./(T-273.15+C))/100;
gamma = wilson(x, a, V, R*T);
% Let’s simply choose Pjsat as Psat(1)
s = 0; % initialize summation
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
T = B(1)/(A(1) - log(100*Pjsat)) - C(1) + 273.15; % temperature in Kelvins
Psat = exp(A-B./(T-273.15+C))/100;
y=x.*gamma.*Psat./(phi.*P);
[BB d] = virial2(T, Tc, Pc, Zc, Vc, w); %virial coefficients needed to calculate fugacity coefficients
phi = virialcoeff(BB, d, T, y, P) % fugacity coefficients
gamma = wilson(x, a, V, R*T);
s = 0;
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
T = B(1)/(A(1) - log(100*Pjsat)) - C(1) + 273.15;
er1 = 1; % initialize error to a value large enough that it can’t be satisfied the first iteration
while er1 > 1e-3
Psat = exp(A-B./(T-273.15+C))/100
y=x.*gamma.*Psat./(phi.'*P); % Transpose phi and check
[BB d] = virial2(T, Tc, Pc, Zc, Vc, w); %virial coefficients needed to calculate fugacity coefficients
phi = virialcoeff(BB, d, T, y, P); % fugacity coefficients
gamma = wilson(x, a, V, R*T);
s = 0; % initialize summation
for i=1:N
s = s + ((gamma(i)*x(i)/phi(i))*(Psat(i)/Psat(1)));
end
Pjsat = P/s;
Tnew = B(1)/(A(1)-log(100*Pjsat))-C(1)+273.15; % T in K
er1 = abs(1-Tnew/T)
T = Tnew;
end
T, y
function phi = virialcoeff(B, delta, T, y, P)
R=83.14; %bar cm^3 mol^-1 K^-1
N=max(size(y));
for k=1:N
s=0;
for i=1:N
for j=1:N
s=s+y(i)*y(j)*(2*delta(i,k)-delta(i,j));
end
end
phi(k) = exp(P/R/T*(B(k,k) + 1/2*s));
end
end
function [BB d] = virial2(T, Tc, Pc, Zc, Vc, w)
% BB contains the virial coefficient of the mixture
% d is delta values (delta = 2*B(j,i) - B(j,j) - B(i,i))
R = 83.14; % bar cm^3 mol^-1 K^-1
N = max(size(Tc));
for i=1:N
for j=1:N
ww(i,j) = (w(i)+w(j))/2;
TTc(i,j) = sqrt(Tc(i)*Tc(j));
VVc(i,j) = ((Vc(i)^(1/3)+Vc(j)^(1/3))/2)^3;
ZZc(i,j) = (Zc(i)+Zc(j))/2;
PPc(i,j) = ZZc(i,j)*R*TTc(i,j)/VVc(i,j);
end
end
TTr = T./TTc;
B0=0.083-0.422./TTr.^1.6;
B1=0.139-0.172./TTr.^4.2;
Bhat = B0+ww.*B1;
BB=Bhat.*TTc*R./PPc;
% calculation for delta(i,j) values
for i=1:max(size(Tc))
for j=1:max(size(Tc))
d(i,j) = 2*BB(i,j) - BB(i,i) - BB(j,j);
end
end
end
function g = wilson(x, a, V, RT)
N=max(size(V)) % number of components
for i=1:N
for j=1:N
L(i,j)=V(j)/V(i)*exp(-a(i,j)/RT);
end
end
for i=1:N
s1=0; % initialize summation
for j=1:N
s1=s1+x(j)*L(i,j);
end
s3=0; % initialize summation
for k=1:N
s2=0; % initialize summation
for j=1:N
s2 = s2 + x(j)*L(k,j);
end
s3 = s3 + x(k)*L(k,i)/s2;
end
lg(i) = 1-log(s1)-s3;
end
g=exp(lg)';
end
Transpose phi in this line
y=x.*gamma.*Psat./(phi.'*P);
to get y as 3x1
More Answers (0)
Categories
Find more on Thermal Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
