problem here
1 次查看(过去 30 天)
显示 更早的评论
% Note :- The code working perfect with single value of of P but %when I want to use multi values like 1:0.5:8 , its says (??? Error %using ==> mpower %Inputs must be a scalar and a square matrix.) help plz
clear all
L = 1;
P = 5.0 % <<<--
T(L) = 270 ;%K
R = 8.314e-3;% MPa m3 / mole.K
Tc = [126.19 304.1 190.4 305.3 369.8 408.05 425.15 460.4 469.8... 507.6 540.2 568.7]';
Tr = T(L)./Tc;
Pc = [3.3978 7.38 4.6 4.8839 4.87 3.648 3.796 3.385326 3.36 3.02... 2.81 2.49]';
omega = [0.0377 0.2236 0.0115 0.0995 0.1523 0.177 0.2002 0.2270...
0.2515 0.3013 0.3495 0.3996]';
s = 0.480 + 1.574.*omega - 0.176.*omega.^2;
alpha = (1 + s.*(1-sqrt(Tr))).^2;
a = 0.42747.*(((R*Tc).^2)./Pc);
b = 0.08664.*(R.*Tc./Pc);
%-----------------------------------
%For Gas phase , using y
%--------------------------------------
q=0;
y = [0.05651 0.00284 0.833482 0.07526 0.02009 0.00305 0.0052 0.0012...
0.000144 0.00068 0.000138 0.00011]';
x = [0.025 0.001 0.8819 0.078 0.01 0.000525 0.0019 0.0006 0.00062...
0.00034 0.000067 0.000054]';
while q==0
% for Vapor Phase using y
sumofA = 0;
rou = zeros(12,1);
for i = 1:12
for j = 1:12
sumofA = sumofA +
y(i).*y(j)*sqrt(a(i).*a(j)*alpha(i).*alpha(j));
rou(i) = rou(i) +
y(j).*sqrt(a(i).*a(j).*alpha(i).*alpha(j));
end
end
sumofB = sum(y.*b);
%----------------------------------------
% For Liquid Phase, using x
%---------------------------------------
sumofAA = 0;
mew = zeros(12,1);
for i=1:12
for j=1:12
sumofAA = sumofAA +
x(i).*x(j)*sqrt(a(i).*a(j)*alpha(i).*alpha(j));
mew(i) = mew(i) +
x(j).*sqrt(a(i).*a(j).*alpha(i).*alpha(j));
end
end
sumofBB = sum(x.*b);
% Calculate the Values of A and B for gas Phase
AV = (sumofA .* P)./(R.^2.*T(L).^2);
BV = (sumofB.*P)./(R.*T(L));
% Calculate the Values of A and B for Liquid Phase
AL = (sumofAA .* P)./(R.^2.*T(L).^2);
BL = (sumofBB.*P)./(R.*T(L));
% Find the compressibility factor for gas Phase
ZV=roots([1 -1 AV-BV-BV^2 -AV*BV]);
Zv=max(ZV);
Zvc=isreal(Zv);
if Zvc == false
Zv=abs(Zv);
disp('error 1')
end
% Find the compressibilty factor for Liquid
ZL=roots([1 -1 AL-BL-BL^2 -AL*BL]);
Zl=min(ZL);
Zlc=isreal(Zl);
if Zlc == false
Zl=abs(Zl);
disp('error 2')
end
%Fagucity of Liquid
phil =exp((b.*(Zl-1)./sumofBB) - log(Zl-BL)-(AL/BL).*
((2*mew/sumofAA)-(b/sumofBB)).*log(1 + BL/Zl));
%disp(phil)
%Fagucity of Vapor
phiv =exp((b.*(Zv-1)./sumofB) - log(Zv-BV)-(AV/BV).*
((2*rou/sumofA)-(b/sumofB)).*log(1 + BV/Zv));
%disp(phiv);
% equilibirum calculation Ki
K = (phil./phiv);
% disp(K)
raw(1)= sum(y./K)-1;
if (abs(raw(1)) < 1e-5 & abs((y./K)-x) < 1e-5)
% disp(K)
break
else
if L==1
L=L+1;
T(L) = 1.01 * T(L-1);
raw(2)=raw(1);
else
L=L+1;
T(L)=T(L-1)-raw(1)/(raw(1)-raw(2)).*(T(L-1)-T(L-2));
raw(2)=raw(1);
end
%adjustment
N = y./K;
x = (1-0.5).*N + 0.5.*x;
end
n=0;
n=n+1;
if n>1000
break
end
end disp(T(L));
1 个评论
回答(1 个)
Jan
2012-3-4
At first your program fails in :
ZV = roots([1 -1 AV-BV-BV^2 -AV*BV]);
This can be fixed - I've inserted commas to improve the readability:
ZV = roots([1, -1, AV-BV-BV.^2, -AV.*BV]);
You need the elementwise power ".^", not the matrix power "^".
But later on the script fails in:
phil = exp((b.*(Zl-1)./sumofBB), ...
-log(Zl-BL)-(AL/BL) .* ...
((2*mew./sumofAA)-(b./sumofBB)) .* log(1 + BL/Zl));
I do not understand, what you want to calculate. This line contains vectors of length 12 and 15. Perhaps you want to use bsxfun to get a [12 x 15] matrix as result, but I cannot know this detail.
You can use the debugger to find such problems by your own:
dbstop if error
Then Matlab stops, when the error occurs and you can inspect the values of the variables in the command window.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!