- 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

# Undefined function or variable C

3 views (last 30 days)

Show older comments

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;

disp("Dew point not found");

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;

disp("bad initial guess")

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

##### 1 Comment

Rik
on 23 May 2019

There are many issues with the code you posted.

And this is already from the first few lines. That last one will already cause the error you got.

### Answers (1)

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.

##### 0 Comments

### See Also

### Community Treasure Hunt

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

Start Hunting!