The following works, but probably doesn't contain all the right data or equations!
% CALCULATE MATERIAL BALANCES FOR THE DISTILLATION TOWER
% INITIAL DATA
Fmass=156000; %input('Feed flow rate Fmass (thousand ton/year) =');
cF=33; %input('Feed concentration cF (vol%) = ');
Vf=26; %input('Select internal vapor ratio (vol%) = ');
densF=0.670; %feed density in ton/cubic meter
anfa=5.68;
% Variable declaration
F=zeros(1,6);
VF=zeros(1,6);
V=zeros(1,6);
LF=zeros(1,6);
L=zeros(1,6);
D=zeros(1,6);
B=zeros(1,6);
% Feed stream
F(3)=0.670; %ton/m3
F(4)= Fmass/F(3)/350/24; %m3/h
F(5)=76.2; %molar weight (MW)
F(6)=(Fmass/350/24)/F(5)*1000; %kmol/h
% MATERIAL BALANCES
% (1) : volume fraction (vol%)
VF(1) =cF+4;
V(1) = Vf;
LF(1) = 100-cF-4;
L(1) = Vf+4;
D(1) =cF;
B(1) = 100-cF;
% (2): volume flow rates in cubic meter / h
VF(2) =F(4)*VF(1)/100;
V(2) = F(4)*V(1)/100;
LF(2) =F(4)*LF(1) /100;
L(2) =F(4)*L(1)/100;
D(2) = F(4)*D(1)/100;
B(2) = F(4)*B(1)/100;
% (3) : density (ton/m3 liquid)
VF(3)=0.591; %density of equilibrium vapor phase in feed.
V(3)=0.598; %density of internal vapor
LF(3)=0.726; %density of equilibrium liquid phase in feed
L(3)=0.615; % density of internal liquid
% (4): mass flow rates (ton/h)
VF(4) =VF(2)* VF(3);
V(4) = V(2)* V(3);
LF(4) = LF(2)*LF(3);
L(4) = L(2)*L(3);
D(4) = VF(4)+V(4)-L(4);
B(4) = LF(4)+L(4) -V(4);
D(3)=D(4)/D(2); %ton/m3
B(3)=B(4)/B(2); %ton/m3
% (5) molar weight (MW)
VF(5) =58.2;
V(5) = 58.3;
LF(5) = 93.3;
L(5) = 60.1;
D(5) =54.5;
B(5) = 93.8;
% (6) molar flow rates (kmole/h)
VF(6) =VF(4)*1000/VF(5); % Vapor phase rate in the feed
V(6) =V(4)*1000/V(5); %Internal vapor flow rate
LF(6) =LF(4)*1000/LF(5); %Liquid phase rate in the feed
L(6) =L(4)*1000/L(5); %Internal liquid flow rate
D(6) =D(4)*1000/D(5); %Distillate flow rate
B(6) =B(4)*1000/B(5); %Bottoms flow rate
% The liquid holdups (kmole)
MD = 13.12;
M = 5.76;
MB = 25.26;
disp('------------------------------------------------------------------');
% STREAM SUMMARY
disp('STREAM SUMMARY (kmole/h)');
disp('Vapor phase rate in the feed: ');disp(VF(6));
disp('Internal vapor flow rate: ');disp(V(6));
disp('Liquid phase rate in the feed: ');disp(LF(6));
disp('Internal liquid flow rate: ');disp(L(6));
disp('Distillate flow rate: ');disp(D(6));
disp('Bottoms flow rate: ');disp(B(6));
% Flow rates below feed location
disp('FLOW RATES BELOW FEED LOCATION HEAVY COMPONENT');
for n=1:8
Ln(n)= L(6) + LF(6);
Vn(n)=V(6);
displn=['Internal liquid rate across stage ' int2str(n) ' is: ' num2str(Ln(n))];
dispvn=['Internal vapor rate across stage ' int2str(n) ' is: ' num2str(Vn(n))];
disp(displn);
disp(dispvn);
end
disp('-----------------------------------------------------------------');
% Flow rates above feed location
disp('FLOW RATES ABOVE FEED LOCATION LIGHT COMPONENT');
for n=9:15
Ln(n)=L(6);
Vn(n) =V(6) + VF(6);
% In the following two lines make sure there is a space immediately before
% the num2str function
displn=['Internal liquid rate across stage ' int2str(n) ' is: ' num2str(Ln(n))];
dispvn=['Internal vapor rate across stage ' int2str(n) ' is: ' num2str(Vn(n))];
disp(displn);
disp(dispvn);
end
disp('------------------------------------------------------------------');
disp('Solving flash equations for feed stream using Newton - Raphson method');
% vapor liquid equilibrium (VLE) relationship:
VLEx2y = @(x)anfa*x/(1+(anfa-1)*x);
%find root of the equation:
g= @(x) VLEx2y(x)-(F*F'-LF*x)/VF; % In this line you hadn't defined zF
% so I just replaced it by F' (ie the
% transpose of F) to get the code
% working. You will need to replace it
% with the proper value
gdot= @(x) VLEx2y(x)+LF/VF;
xF=0.1;
while abs(g(xF)/xF)>0.00001
xF=xF-g(xF)/gdot(xF);
end
yF=VLEx2y(xF);
dispxf=['The concentration of equilibrium liquid in the feed ' ' is,XF: ' num2str(xF)];
dispyf=['The concentration of equilibrium vapor in the feed ' ' is,YF: ' num2str(yF)];
disp(dispxf);
disp(dispyf);
