# control and simulation of petroleum distillation column

5 views (last 30 days)
Itqan Ismail on 17 Jan 2023
Commented: Itqan Ismail on 31 Jan 2023
% 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);
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:
%y=VLEx2y(x)=anfa*x/(1+(anfa-1)*x);
% find root of the equation: g(x)= VLEx2y(x)-(F*zF-LF*x)/VF = 0
% g'(x)= VLEx2y'(x)+LF/VF
xF=0.1;
while abs(g(xF,D(6),B(6),LF(6),VF(6))/xF)>0.00001
xF=xF-g(xF,D(6),B(6),LF(6),VF(6))/gdot(xF,LF(6),VF(6));
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);
I couldnt run this code since there are error appeared. Anyone know on how to fix this. Thank you in advanced

Alan Stevens on 17 Jan 2023
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)');
STREAM SUMMARY (kmole/h)
disp('Vapor phase rate in the feed: ');disp(VF(6));
Vapor phase rate in the feed: 104.1446
disp('Internal vapor flow rate: ');disp(V(6));
Internal vapor flow rate: 73.9225
disp('Liquid phase rate in the feed: ');disp(LF(6));
Liquid phase rate in the feed: 135.8833
disp('Internal liquid flow rate: ');disp(L(6));
Internal liquid flow rate: 85.0927
disp('Distillate flow rate: ');disp(D(6));
Distillate flow rate: 96.4555
disp('Bottoms flow rate: ');disp(B(6));
Bottoms flow rate: 143.7346
% Flow rates below feed location
disp('FLOW RATES BELOW FEED LOCATION HEAVY COMPONENT');
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
Internal liquid rate across stage 1 is: 220.976
Internal vapor rate across stage 1 is: 73.9225
Internal liquid rate across stage 2 is: 220.976
Internal vapor rate across stage 2 is: 73.9225
Internal liquid rate across stage 3 is: 220.976
Internal vapor rate across stage 3 is: 73.9225
Internal liquid rate across stage 4 is: 220.976
Internal vapor rate across stage 4 is: 73.9225
Internal liquid rate across stage 5 is: 220.976
Internal vapor rate across stage 5 is: 73.9225
Internal liquid rate across stage 6 is: 220.976
Internal vapor rate across stage 6 is: 73.9225
Internal liquid rate across stage 7 is: 220.976
Internal vapor rate across stage 7 is: 73.9225
Internal liquid rate across stage 8 is: 220.976
Internal vapor rate across stage 8 is: 73.9225
disp('-----------------------------------------------------------------');
-----------------------------------------------------------------
% Flow rates above feed location
disp('FLOW RATES ABOVE FEED LOCATION LIGHT COMPONENT');
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
Internal liquid rate across stage 9 is: 85.0927
Internal vapor rate across stage 9 is: 178.0671
Internal liquid rate across stage 10 is: 85.0927
Internal vapor rate across stage 10 is: 178.0671
Internal liquid rate across stage 11 is: 85.0927
Internal vapor rate across stage 11 is: 178.0671
Internal liquid rate across stage 12 is: 85.0927
Internal vapor rate across stage 12 is: 178.0671
Internal liquid rate across stage 13 is: 85.0927
Internal vapor rate across stage 13 is: 178.0671
Internal liquid rate across stage 14 is: 85.0927
Internal vapor rate across stage 14 is: 178.0671
Internal liquid rate across stage 15 is: 85.0927
Internal vapor rate across stage 15 is: 178.0671
disp('------------------------------------------------------------------');
------------------------------------------------------------------
disp('Solving flash equations for feed stream using Newton - Raphson method');
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);
The concentration of equilibrium liquid in the feed is,XF: 642.7001
disp(dispyf);
The concentration of equilibrium vapor in the feed is,YF: 1.2133
Itqan Ismail on 31 Jan 2023
Thank you so much for your response