matlab 2019b, Peng Robinson for loop - repeat y values, Project

how come when I run my script the y values are repeated?

input script:

clc; clear; close all; clear all

%methane is denoted with 1, n-pentane is denoted with 2

T = 37.78 + 273.15;

R = 0.00008314; %m.^3.*Bar./K.*mol

%accentric factors:

w1 = 0.008; %methane

w2 = 0.251; %n-pentane

%Critical temps and pressures:

Tc1 = 190.6; % K, methane

Tc2 = 469.6; % K, n-pentane

Pc1 = 46.00; % Bar, methane

Pc2 = 33.74; %Bar, n-pentane

Tr1 = T./Tc1;

Tr2 = T./Tc2;

%Antoine coefficients:

%Methane:

A1 = 8.6041;

B1 = 897.84;

C1 = -7.16;

%n-Pentane

A2 = 9.2131;

B2 = 2477.07;

C2 = -39.94;

%Peng-Robinson EOS parameters:

a1 = 0.45724.*((R.^2.*Tc1.^2)./Pc1);

b1 = 0.07780.*((R.*Tc1)./Pc1);

kapa1 = 0.37464 + 1.5422.*w1 - 0.26992.*w1.^2;

alfa1 = (1+kapa1.*(1-sqrt(Tr1))).^2;

a2 = 0.45724.*((R.^2.*Tc2.^2)./Pc2);

b2 = 0.07780.*((R.*Tc2)./Pc2);

kapa2 = 0.37464 + 1.5422.*w2 - 0.26992.*w2.^2;

alfa2 = (1+kapa2.*(1-sqrt(Tr2))).^2;

tol = 0.00001;

xx = linspace(0,.99, 50);

PP = zeros(4,50); % 4 values of k, 50 values of x

Y1 = zeros(4,50);

Y2 = zeros(4,50);

kk = [0, 0.025, 0.05, 0.10];

Psat1 = exp(A1 - (B1./(T + C1)));

Psat2 = exp(A2 - (B2./(T + C2)));

%-------------Values used to test while loop

%y1 = 0.849;

%y2 = 0.14;

%x1 = 0.735;

%x2 = 0.265;

%P = 172.7;

for j = 1:4

k = kk(j);

%j = 1;-------------used to test "i" for loop

for i = 1:length(xx)

x1 = xx(i);

x2 = 1- x1;

P = x1.*Psat1 + x1.*Psat2;

y1 = (x1.*Psat1)./P %----right here is where im having problems, if you run the code you'll see that the initial values here already equal 1 before going into the next loop.

y2 = (x1.*Psat2)./P %--If you un-mute the vectors for y1 and y2 at the bottom of the code you'll see that there's an issue with how those matrices are being made as well, but i think that's happening because of this first problem up here.

err = 1-(y1+y2);

while err > tol

aalfa_v = y1.^2.*(a1.*alfa1) + 2.*y1.*y2.*sqrt((a1.*alfa1).*(a2.*alfa2)).*(1-k) + y2.^2.*(a2.*alfa2);

bmix_v = y1.*b1 + y2.*b2;

aalfa_L = x1.^2.*(a1.*alfa1) + 2.*x1.*x2.*sqrt((a1.*alfa1).*(a2.*alfa2)).*(1-k) + x2.^2.*(a2.*alfa2);

bmix_L = x1.*b1 + x2.*b2;

Av = (aalfa_v.*P)./(R.^2.*T.^2);

Bv = (bmix_v.*P)./(R.*T);

AL = (aalfa_L.*P)./(R.^2.*T.^2);

BL = ((bmix_L.*P)./(R.*T));

%Try to make Z=Pv./RT to solve for v rather than Z used pg 254

%for an example and layout of Cubic EOS on pgs 250-251

% syms v

% a1EOS = (.45724.*R.^2.*Tc1.^2)./Pc1

% b1EOS = .07780.*R.*Tc1./Pc1

% a2EOS = .45724.*R.^2.*Tc2.^2./Pc2.^2

% b2EOS = .07780.*R.*Tc2.^2./Pc2

% aEOSmix = y1.^2.*a1EOS+2.*y1.*y2.*sqrt(a1EOS.*a2EOS)+y2.^2.*a2EOS

% bEOSmix = y1.*b1EOS+y2.*b2EOS

% vroots = fzero((R.*T)./(v-bEOSmix)-(aEOSmix.*(aalfa1+aalfa2)./(v.^2+2.*bEOSmix.*v-bEOSmix.^2)))

% v = roots((1+Bv./v+Cv./v.^2+Dv./v.^3).*((R.*T)./P))

comp_v = roots([1 -(1-Bv) (Av-2*Bv-3*Bv^2) -(Av*Bv-Bv^2-Bv^3)]); %need to fix polynomial to solve for v instead of Z

Z(1) = max(comp_v); %vapor --Aryana's suggestion

comp_L = roots([1 -(1-BL) (AL-2.*BL-3.*BL.^2) -(AL.*BL-BL.^2-BL.^3)]);

Z(2) = min(comp_L); %liquid

%phi's with volume

%phi1v = exp((b1./bmix).*(z-1) - log(((vL-bmix).*P)./(R.*T)) + aalfa./(2.*sqrt(2).*bmix.*R.*T).*((b1./bmix)-(2./aalfa).*(y1.*aalfa + y2.*aalfa)).*log((vL + (1+sqrt(2)).*bmix)./(vL + (1-sqrt(2)).*bmix)));

%phi1L = exp((b1./bmix).*(z-1) - log(((vv-bmix).*P)./(R.*T)) + aalfa./(2.*sqrt(2).*bmix.*R.*T).*((b1./bmix)-(2./aalfa).*(y1.*aalfa + y2.*aalfa)).*log((vv + (1+sqrt(2)).*bmix)./(vv + (1-sqrt(2)).*bmix)));

%phi2v = exp((b2./bmix).*(z-1) - log(((vL-bmix).*P)./(R.*T)) + aalfa./(2.*sqrt(2).*bmix.*R.*T).*((b2./bmix)-(2./aalfa).*(y1.*aalfa + y2.*aalfa)).*log((vL + (1+sqrt(2)).*bmix)./(vL + (1-sqrt(2)).*bmix)));

%phi2L = exp((b2./bmix).*(z-1) - log(((vv-bmix).*P)./(R.*T)) + aalfa./(2.*sqrt(2).*bmix.*R.*T).*((b2./bmix)-(2./aalfa).*(y1.*aalfa + y2.*aalfa)).*log((vv + (1+sqrt(2)).*bmix)./(vv + (1-sqrt(2)).*bmix)));

%phi's with compressibility factor

phi1v = exp((b1./bmix_v).*(Z(1)-1) - log(Z(1)-Bv) + aalfa_v./(2.*sqrt(2).*bmix_v.*R.*T).*((b1./bmix_v)-(2./aalfa_v).*((y1.*a1)+(y2.*sqrt(a1.*a2).*(1-k)))).*log((Z(1)+(1+sqrt(2)).*Bv)./(Z(1)+(1-sqrt(2)).*Bv)));

phi1L = exp((b1./bmix_L).*(Z(2)-1) - log(Z(2)-BL) + aalfa_L./(2.*sqrt(2).*bmix_L.*R.*T).*((b1./bmix_L)-(2./aalfa_L).*((y1.*a1)+(y2.*sqrt(a1.*a2).*(1-k)))).*log((Z(2)+(1+sqrt(2)).*BL)./(Z(2)+(1-sqrt(2)).*BL)));

phi2v = exp((b2./bmix_v).*(Z(1)-1) - log(Z(1)-Bv) + aalfa_v./(2.*sqrt(2).*bmix_v.*R.*T).*((b2./bmix_v)-(2./aalfa_v).*((y1.*a1)+(y2.*sqrt(a1.*a2).*(1-k)))).*log((Z(1)+(1+sqrt(2)).*Bv)./(Z(1)+(1-sqrt(2)).*Bv)));

phi2L = exp((b2./bmix_L).*(Z(2)-1) - log(Z(2)-BL) + aalfa_L./(2.*sqrt(2).*bmix_L.*R.*T).*((b2./bmix_L)-(2./aalfa_L).*((y1.*a1)+(y2.*sqrt(a1.*a2).*(1-k)))).*log((Z(2)+(1+sqrt(2)).*BL)./(Z(2)+(1-sqrt(2)).*BL)));

%remember to double check phi equations

%update y

y1 = (x1.*phi1L)./phi1v;

y2 = (x1.*phi2L)./phi2v;

%update P

P = P.*(y1+y2);

err = 1 - (y1+y2);

end

%these allow storing of values to create plots

PP(j,i) = P;

y1(j,i) = y1;

y2(j,i) = y2;

end

end

% figure(1)

% plot(y1,P)

% hold on

% plot(y2,P)

% hold off

output:

y1 =

NaN

y2 =

NaN

y1 =

0.9962

y2 =

0.0038

y1 =

0.9962

y2 =

0.0038

y1 =

0.9962

y2 =

0.0038

y1 =

0.9962

…

[A Emmert]

the KK are binary coeffecients

[Rik]

Are you going to delete this question as well? Also, please use the layout tools to make your code more readable.