Error in PengRobinson (line 21)

1 回表示 (過去 30 日間)
Hamed
Hamed 2020 年 4 月 30 日
コメント済み: Mehmed Saad 2020 年 4 月 30 日
I keep getting this error for both line 21 (P/Pc) and 20 T/Tc, I am trying to calculate using values of P from 10 to 3000 at 200 steps. And then plot the answer. I am new to Matlab so any help would be appreciated ( I got the code from file exchanger) Thank you.
% PengRobinson.m : calculates the compressibility factor,fugacity coefficient and density
% of a pure compound with the Peng Robinson equation of state (PR EOS)
% Authors: Piñero.R, Serna.J.G, Martin. A.
%
% function result = PengRobinson(T,P,Tc,Pc,w,MW,Liquido)
% Parameters: T,P,w,Tc,Pc,w,MW,Liquido
% T: Temperature [=] K
% P: Presure [=] Pa
% Tc: critical temperature [=] K
% Pc: critical presure [=] Pa
% w: accentic factor
% MW: molar weigth [=] kg/mol
% Liquido: if Liquido == 1, then calculates liquid fugacity;
% if Liquido == 0 then calculates vapor fugacity
% Example:
% [Z fhi density] = PengRobinson(273,2*1.013*1e5,304.21,7.382*1e6,0.225,0.044,1)
function [Z,fhi,density] = PengRobinson(T,P,Tc,Pr,w,MW,Liquido)
R = 8.314; % gas constant [=] J/(mol K)
% Reduced variables
Tr = T/Tc;
Pr = P/Pr;
% Parameters of the EOS for a pure component
m = 0.37464 + 1.54226*w - 0.26992*w^2;
alfa = (1 + m*(1 - sqrt(Tr)))^2;
a = 0.45724*(R*Tc)^2/Pr*alfa;
b = 0.0778*R*Tc/Pr;
A = a*P/(R*T)^2;
B = b*P/(R*T);
% Compressibility factor
Z = roots([1 -(1-B) (A-3*B^2-2*B) -(A*B-B^2-B^3)]);
ZR = [];
for i = 1:3
if isreal(Z(i))
ZR = [ZR Z(i)];
end
end
if Liquido == 1
Z = min(ZR);
else
Z = max(ZR);
end
% Fugacity coefficient
fhi = exp(Z - 1 - log(Z-B) - A/(2*B*sqrt(2))*log((Z+(1+sqrt(2))*B)/(Z+(1-sqrt(2))*B)));
if isreal(fhi)
density=P*MW/(Z*R*T);
result = [Z fhi density];
else
'No real solution for "fhi" is available in this phase'
result=['N/A' 'N/A' 'N/A'];
end
% Bibliography:
% - ORBEY. H, SANDLER. I; Modeling Vapor-Liquid Equilibria: cubic equations of state and their mixing rules; Cambridge University Press (1998)
% - Walas,Stanley. M ; Phase Equilibria in Chemical Engineering ; Boston,
% Butterworth Publishers (1984)

回答 (1 件)

Mehmed Saad
Mehmed Saad 2020 年 4 月 30 日
編集済み: Mehmed Saad 2020 年 4 月 30 日
You are running a function directly without passing arguments to it. If you want to run the function PengRobinson you have to call it from another code or command window
Example
Type following i command window
[Z fhi density] = PengRobinson(273,2*1.013*1e5,304.21,7.382*1e6,0.225,0.044,1)
I am assuming that PengRobinson is added in MATLAB path or is in Current folder
Also see the following link
  2 件のコメント
Hamed
Hamed 2020 年 4 月 30 日
Yes, however I dont want a fixed value for P, I want it to be from 10 to 20000 step 200, Do you have any ideas on how to do that?
Thank you
Mehmed Saad
Mehmed Saad 2020 年 4 月 30 日
Use for loop

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeThermal Analysis についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by