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 visualizzazione (ultimi 30 giorni)
Tony Mei
Tony Mei il 22 Lug 2022
Commentato: Tony Mei il 22 Lug 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 Commenti
the cyclist
the cyclist il 22 Lug 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 il 22 Lug 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)';

Accedi per commentare.

Accedi per rispondere a questa domanda.

Risposta accettata

VBBV
VBBV il 22 Lug 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

Più risposte (0)

Categorie

Scopri di più su Thermal Analysis in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by