control and simulation of petroleum distillation column
Afficher commentaires plus anciens
% 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
Réponse acceptée
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);
1 commentaire
Itqan Ismail
le 31 Jan 2023
Plus de réponses (0)
Catégories
En savoir plus sur Distillation Design dans Centre d'aide et File Exchange
le 17 Jan 2023
le 31 Jan 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
