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

1 view (last 30 days)
Tony Mei
Tony Mei on 22 Jul 2022
Commented: Tony Mei on 22 Jul 2022
% 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
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
Tony Mei on 22 Jul 2022
apologies, here is everything!
% 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
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
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
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)';

Sign in to comment.

Sign in to answer this question.

Accepted Answer

VBBV
VBBV on 22 Jul 2022
Ran in:
% 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);
N = 3
% 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
phi = 1×3
0.9740 0.9541 0.9672
gamma = wilson(x, a, V, R*T);
N = 3
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
Psat = 3×1
0.7450 2.9645 2.6766
N = 3
er1 = 6.0074e-04
T, y
T = 364.5202
y = 3×1
0.1132 0.3406 0.5538
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

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!