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

% 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

[the cyclist]

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]

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)';

Answer 1

[VBBV]

% 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

[Tony Mei]

May I ask why we transpose that ?@VBBV