Matlab help regarding code

2 次查看（过去 30 天）
Nasir Qazi
Nasir Qazi 2012-2-11
编辑： dav 2013-10-9
% Calculate the vapor Pressure of n-heptane using SRK-EOS
% P = RT / V - b - alpha*a / V(V-b)
R = 8.134e-3; % MPa m3 / mole.K
n = 0;
n = n+1;
for i=1:5
Pn = 1.6714;
Tr = 0.9;
Tc = 540.13; % K
Pc = 2.74; % Mpa;
T = Tr*Tc;
w = 0.350;
m = 0.480+1.574*w -0.176*w^2;
al=(1+m*(1-Tr^0.5))^2;
a = 0.42748*(((R*Tc)^2) /Pc);
b = 0.08664*(R*Tc / Pc);
% Z3 - Z2 + (A-B-B^2)Z-AB=0
A = (a*Pn*al)/(R*T)^2;
B = (b*Pn)/(R*T);
Zv=max(roots([1 -1 A-B-B^2 -A*B]));
Zl=min(roots([1 -1 A-B-B^2 -A*B]));
disp(Zv)
disp(Zl)
phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));
phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));
Fv=Pn*phiv;
Fl=Pn*phil;
disp(phiv)
disp(phil);
disp(Fv);
disp(Fl);
A=Fl-Fv;
if (A < 10^-4);
Pnew=Pn;
else
Pn+1 =Pn * (phil/phiv);
end
end
Basically in this code I am trying to assume a pressure 'Pn' to calculate the results which satisfy this relation (A<10^-4) if this satisfy's then our assume pressure is right else if not then we use the relation Pn+1 = Pn*(Phil/Phiv) to caluclate another pressure and update the above pressure (previously assume) with this one. I have error with this and need advice please help.
regards

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2012-2-11
The syntax
Pn+1 =Pn * (phil/phiv);
is not valid. The left side of an assignment must be a valid variable name or array reference. For example,
Pn_plus_1 =Pn * (phil/phiv);
Or
P(n+1) =Pn * (phil/phiv);
In your "while" loop, you do not recalculate "A", so if the loop is entered at all, it is going to be repeated indefinitely. When you assign a variable value, the assignment does not affect any values that were already calculated before that point, including not any assignments that used that same variable name in the formula.
foo = 5;
bar = foo;
foo = 7;
bar
At the end of this, foo will be 7 and bar will still be what it was at the time the assignment for it was executed, namely 5. The update to the variable "foo" that was used in initializing "bar" does not cause an update to any variable already assigned to.
  5 个评论
显示 3更早的评论
Image Analyst
Image Analyst 2012-2-12
I still don't understand. Maybe you can give a Pn where it does get into that part of the if statement and then bails out of the loop. But I see no reason why it should exit the loop and stop the program. There is no break statement anywhere so why should it exit the loop or "stop there" whatever that means. You mean like stop with an error because an exception got thrown? And how/where is Pnew ever used, even if it did get assigned?
Nasir Qazi
Nasir Qazi 2012-2-12
% these are the correct codes, have a look again
% Calculate the vapor Pressure of n-heptane using SRK-EOS
% P = RT / V - b - alpha*a / V(V-b)
R = 8.134e-3; % MPa m3 / mole.K
n = 0;
n = n+1;
for i=1:5
Pn = 1.6714; % Assume Pressure
Tr = 0.9;
Tc = 540.13; % K
Pc = 2.74; % Mpa;
T = Tr*Tc;
w = 0.350;
m = 0.480+1.574*w -0.176*w^2;
al=(1+m*(1-Tr^0.5))^2;
a = 0.42748*(((R*Tc)^2) /Pc);
b = 0.08664*(R*Tc / Pc);
% Z3 - Z2 + (A-B-B^2)Z-AB=0
A = (a*Pn*al)/(R*T)^2;
B = (b*Pn)/(R*T);
Zv=max(roots([1 -1 A-B-B^2 -A*B]));
Zl=min(roots([1 -1 A-B-B^2 -A*B]));
disp(Zv)
disp(Zl)
phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));
phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));
Fv=Pn*phiv;
Fl=Pn*phil;
disp(phiv)
disp(phil);
disp(Fv);
disp(Fl);
error=(Fl-Fv);
if (error < 10^-4); % this statement has to be verify if not
Pn; % |
else % |
Pn+1 =Pn * (phil/phiv); % V
P(n+1) = Pn ( this pressure should update the previous assume pressure )
end
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

尚未输入任何标签。

Community Treasure Hunt

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

Start Hunting!

Translated by