# Undefined function or variable C

3 views (last 30 days)
Answered: Guillaume on 23 May 2019
I've been given a function by my professor that calculates PT flash for non -ideal systems. I'm trying to use it to simulate an experiment. But continously get error.
I think it may be the way i've set my input up.
The sub-functions are given below:
function [PD,xDew]=PDew_Raoult(Ps,y,T)
% compute dew pressure at defined T for a mixture with a user given vapour
% pressure using the Raoult Law
c=length(y);
for i=1:c
Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
PD=1/sum(y./Ps);
xDew=y*PD./Ps;
end
function[PB,Ybub]=Pbubble_Raoult(Ps,x)
%this is a function that calculates the bubble point of a fluid with given
%vapour pressure and composition
Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)))
i=length(x);
PB=Ps.*x;
Ybub=Ps.*x/PB;
end
function [ alphaV ] = RACHFORDRICE_BISECT( K,Z )
%This function applies the bisection method to the Rachford-Rice equation
%for alphaV in (0,1).
%Input/output are on molar basis.
a=0;
b=1;
%psi_0=sum(z.*(K-1))
%psi_1=sum((z.*(K-1))./K)
while b-a>0.000001
alphaV=(a+b)/2;
psi_alphaV=sum((Z.*(K-1))./((1-alphaV)+alphaV*K));
psi_b=sum((Z.*(K-1))./((1-b)+b*K));
if psi_alphaV*psi_b<0
a=alphaV;
else
b=alphaV;
end
end
function[PB,Ybub]=PB_VLE_Wilson(Z,C,MW,rhoL,T,BIP)
% this is a function that calculates the bubble point using the wilson
% activity coeefficient model.
c=length(Z);
for i=1:c
Psi=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
[gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z);
PB=sum(Psi.*Z.*gamma);
Ybub=Psi.*gamma/PB;
end
function [PD,xDew]=PD_VLE_wilson (C,MW,rhoL,BIP,T,Z)
% this is a function that calculates the dewpoint of a mixture with the
% wilson correlation for generating activity coefficients for non ideal
%liquid phase
c=length(Z);
for i=1:c
Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
end
% initial guess of Pdew using PD_VLE_ID(C,T,y)
[PD,xDew]=PDew_Raoult(Ps,Z,T); %initial guess of PD(1) and calculation of x1
[gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z);% calculation of gamma at the initial guess of x1
iter=0;
while max(abs((PD*Z-Ps.*x.*gamma)./(PD*Z)))>0.001 &&iter<1000 %isofugacity conditions and control on the number of iterations
PD=1/sum(Z./Ps.*gamma);
xDew=(PD*Z)./(Ps.*gamma);
[gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,x);% initial guess of gamma
iter =iter +1;
end
if iter<1000
PD=1/sum(Z./Ps.*gamma);
xDew=(PD*Z)./(Ps.*gamma);
else
PD=0;
xDew=0;
end
end
function [gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z)
%This program calculates the activity coefficients (gamma) and the
%activities (a) of each component of a mixture of c components using the
%Wilson model.
%INPUT PARAMETERS: MW: vector (1xc) reporting the molecular weights of the
%c components; rhoL, vector (1xc) reporting the liquid density of the c
%pure components at temperature T; BIP is a matrix cxc reporting the energy
%interaction parameters (BIP(i,j)=lambda_ij-lambda_ii, J/mol). The energy
%interaction parameters are temperature independent; T: temperature of the
%system; x vector (1xc) reporting the mole fractions of the components of
%the mixture.
%OUTPUT PARAMETERS: gamma: vector 1xc reporting the activity
%coefficients of the components of the mixture; a: vector 1xc reporting the
%activities of the components of the mixture.
%Unless otherwise stated, all input/output parameters are expressed
%according to MKS.
R=8.314;
c=length(Z);
%Molar volumes of the pure liquid components composing the mixture
Vl=1./((rhoL*1000)./MW);
%Lambda terms (dimensionless) of the Wilson formula
for i=1:c
for j=1:c
Lambda(i,j)=(Vl(j)/Vl(i))*exp(-BIP(i,j)/(R*T));
end
end
for i=1:c
for j=1:c
A=sum(Z.*Lambda(j,:));
C(j)=Z(j)*Lambda(j,i)/A;
end
lngamma(i)=1-log(sum(Z.*Lambda(i,:)))-sum(C);
gamma(i)=exp(lngamma(i));
a(i)=gamma(i)*Z(i);
end
end
%% the main function to be ran is given as
function[Xeq, Yeq, alphaV, Fl, Fv]=PTFLASH_VLE_Wilson(C,MW,rhoL,BIP,P,T,Z)
%this is a function that calculates PT flash for non -ideal systems where K
%cannot be explicit since it depends on X which is also a variable
c=length(Z);
for i=1:c
Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
end
% validating the possibility of having VLE
[PB,y]=PB_VLE_Wilson(Z,C,MW,rhoL,T,BIP);
[PD,x]=PD_VLE_wilson (C,MW,rhoL,BIP,T,Z);
% checking if there is VLE
if P<=PD %% y here represents the mole fraction of componenets
Xeq=0;
Yeq=y;
alphaV=1;
Fl=0;
Fv=P*Z;
elseif P>=PB %% x here represents the mole fraction of components
Xeq=x;
Yeq=0;
alphaV=0;
[gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Xeq);
Fl=Ps.*Xeq.*gamma;
Fv=0;
else
Xeq=(x+y)/2 ; % first guess on x
[gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Xeq); %first guess on gamma
K=(Ps.*gamma)/P; % first guess on K
Psi_0=sum(z.*(K-1));
Psi_1=sum(z.*(K-1)./K);
if Psi_0*Psi_1>=0 % bad initial guess
Xeq=0;
Yeq=0;
alphaV=0;
Fl=0;
Fv=0;
else %supposedly good initial guess
Fl=zeros(1,c); % trick to enter the while cycle at the begininng
Fv=Fl+1;
iter=0; % initialization of iteration count
while max(abs((Fv-Fl)./Fv))>0.00001 && iter<10000 && Psi_0*Psi_1<0
[ alphaV ] = RACHFORDRICE_BISECT( K,Z );
xeq=Z./((1-alphaV)+alphaV*K);
Yeq=K.*xeq;
[gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z); %Gamma new guess
K=(Ps.*gamma)/P;
psi_0=sum(z.*(K-1));
psi_1=sum(z.*(K-1)./K);
Fv=P*Yeq;
Fl=Ps.Xeq.*gamma;
iter=iter+1;
end
end
end
end
Rik on 23 May 2019
There are many issues with the code you posted.
• There are some extra or missing end keywords, so we can't run this code at all
• You didn't provide reasonable inputs for PDew_Raoult
• There are still several mlint warnings
• The indentation suggests that some code was intended to be in a loop
• You assing a value to c (lower case), but then also use C (upper case) without defining it
And this is already from the first few lines. That last one will already cause the error you got.

Guillaume on 23 May 2019
If this exactly the code you've been given without you having made any modification:
function [PD,xDew]=PDew_Raoult(Ps,y,T)
% compute dew pressure at defined T for a mixture with a user given vapour
% pressure using the Raoult Law
c=length(y);
for i=1:c
Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
then it was never working and you need to go back to your professor to ask for code that actually work.
C is indeed undefined (and clearly not meant to be the same as the lowercase c). Until you know what C is meant to be, there's nothing you can do to fix it. I would suspect that C is meant to be y with y a Nx5 matrix (with N>5 otherwise, you'll hit a bug in the code!)
Also completely absurd is that the function takes Ps as an input and then immediately discard that input to replace it with a calculated value.
Also, completely absurd is tha the function calculates values in a loop and at each step of the loop, completely throws away the values calculated in the previous iteration. In the end, it only keeps the values calculated in the last, the previous steps having been a waste of time.
So, again, that code never worked, even if it had worked, it had bugs, made no sense and was poorly written.