Jump to content



Featured Articles

Check out the latest featured articles.

File Library

Check out the latest downloads available in the File Library.

New Article

Product Viscosity vs. Shear

Featured File

Vertical Tank Selection

New Blog Entry

Low Flow in Pipes- posted in Ankur's blog

0

Designing A Pbr In Matlab For Ethylbenzene Production

reactor matlab fixed bed ode catalyst kynetics alkylation

2 replies to this topic
Share this topic:
| More

#1 catarinajorge7

catarinajorge7

    Brand New Member

  • Members
  • 1 posts

Posted 30 November 2017 - 01:30 PM

Hi, I'm a chemical engineering student doing a project about ethylbenzene production and right now, I'm designing the reactors of the process. 

I'm doing it on matlab to solve the dF/dz at the same time as using Ergun's equation to solve pressure drops but my code is not working and I can't seem to find why. The problem is that my code is based on calculating the conversion of the reaction every step the ode makes but it's not actually doing anything so my results as molar flows are always constant.

 

 

I'll post my code :)

 

clear all
close all
clc
 
%% Definicao variaveis globais
global compostos caudal cat leito T P_i R reactor densidade
 
%% Lista dos compostos
compostos.nome=["Etileno";"Benzeno";"Etilbenzeno";"Dietilbenzeno";"Tolueno";"Metano";"Etano"];
 
%% Massas moleculares
compostos.Mw=[28.05316,...
    78.11184,...   
    106.165,...
    134.2182,...
    92.13842,...
    16.04246,...
    30.06904]; % Massas moleculares  g/mol etileno benzeno etilbenzeno deb
%% Caudais Molares em kmol/h, Massas Moleculares, T's, Press�es
caudal.molar_i=[383.8771,... %Caudais molares da alimentacao em kmol/h;
    1535.508,...
    43.23589,...
    7.513956,...
    0.30148,...
    0.825736,...
    1.006617];
caudal.molar_total=sum(caudal.molar_i); %kmol/h
 
%% Caudais massicos em kg/h, volumetricos e afins
caudal.massico_i=[10768.96,...
    119941.4,...
    4590.139,...
    1008.509,...
    27.77787,...
     13.24684,...
    30.268];
caudal.massico_total=sum(caudal.massico_i); %kg/h
 
 
%% Temperatura da Alimentacao, K (isotermico)
T=172+273.15; 
 
%% Pressao inicial
P_i=28837.1; %Pressao inicial kPa
P_atm=28;
R.atmLKkmol=82.057; %atm.L/(K.kmol)
R.kPaLKkmol=8314.42553; %L.kPa/(K.kmol)
R.kPam3Kkmol=8.314462;%m^3.kPa/(K.kmol)
R.atmLKmol=0.082057;%atm.L/(K.mol)
%% Propriedades do catalisador e do leito
cat.D=4;  %Di�metro do reator em m
cat.r=cat.D/2; %Raio do reator em m
cat.L=29.12; %Comprimento do leito de catalisador em m
cat.Dp=0.034*10^-2; %Di�metro das part�culas de catalisador em m
cat.porosidade=0.5; %Porosidade do particula;
cat.densidade=750; %kg/m3 densidade da particula de catalisador
leito.porosidade=0.25; %valor assumido por 20% vazio no reator (Vvazio/Vleito)
leito.densidade=cat.densidade/(1-leito.porosidade); %Densidade do leito em kg/m^3
 
 
cat.area=pi*cat.r^2; %m^2
cat.massa=(1-leito.porosidade)*cat.area*cat.L*cat.porosidade%Massa de Catalisador em kg
reactor.volume=cat.massa/leito.densidade %volume do reator em m^3
 
%% Constantes para calculo das densidades g/ml
% modified form of the Rackett equation
% DEN = A B ^-(1-Tr)^n
 
%                  A      B       n   Tc /K
compostos.const=[  0      0       0   282.5;...
                 0.3009 0.2677 0.2818 562.16;...
                 0.2889 0.2644 0.2921 617.17;...
                 0.2747 0.2534 0.2857 657.9;...
                 0.29999 0.2711 0.2989 591.79;...
                   0      0       0    190.85;...
                   0      0       0    305.3 ];
 
 
%%calculo da densidade da mistura inicial
densidade.liquidos_i=calculodens(T, compostos.const(2:5))';
densidade.etileno_i=calculodensg(compostos.Mw(1),caudal.massico_i(1), P_i, T, R.kPam3Kkmol);
densidade.outrosg_i=calculodensg(compostos.Mw(6:7),caudal.massico_i(6:7), P_i, T, R.kPam3Kkmol);
densidade.gas_i=sum([densidade.etileno_i densidade.outrosg_i]);
densidade.mix_i=sum((caudal.massico_i(2:5)./caudal.massico_total).*densidade.liquidos_i+densidade.gas_i) % kg/m^3
 
%caudal volumetrico
caudal.volum_i=caudal.massico_total/densidade.mix_i; 
%% chamar a função ode
cond=[caudal.molar_i(1:4) P_i];
zspan=[0 cat.L];
[Z F]=ode45(@sedo,zspan,cond);
 
figure(1)
subplot(2,1,1)
plot(Z,F(:,(1:4)))
xlabel('L (m)')
ylabel('Fi (kmol/h)')
legend('etileno','benzeno','etilbenzeno','deb')
subplot(2,1,2)
plot(Z,F(:,5))
xlabel('L (m)')
ylabel('P (kPa)')
legend('P')
 
 
And my function for the ode:
 
function [edo]=sedo(z,f)
%% Definicao variaveis globais
global compostos caudal cat leito T P_i R reactor densidade
 
%%
edo=f;
 
%"Etileno";"Benzeno";"Etilbenzeno";"Dietilbenzeno";"Tolueno";"Metano";"Etano"
caudal.molar_iter=[f(1) f(2) f(3) f(4) caudal.molar_i(5) caudal.molar_i(6) caudal.molar_i(7)];
caudal.massico_iter=[f(1)*compostos.Mw(1) f(2)*compostos.Mw(2) f(3)*compostos.Mw(3) f(4)*compostos.Mw(4) caudal.massico_i(5) caudal.massico_i(6) caudal.massico_i(7)];
P=f(5);
% calculo caudal molar total
caudal.molartot_iter=sum(caudal.molar_iter);
% calculo caudal massico total
caudal.massicotot_iter=sum(caudal.massico_iter);
 
%calculo densidade da mistura
densidade.liquidos_iter=calculodens(T, compostos.const(2:5))';
densidade.etileno_iter=calculodensg(compostos.Mw(1),caudal.massico_iter(1), P , T, R.kPam3Kkmol);
densidade.outrosg_iter=calculodensg(compostos.Mw(6:7),caudal.massico_iter(6:7), P , T, R.kPam3Kkmol);
densidade.gas_iter=sum([densidade.etileno_iter densidade.outrosg_iter]);
densidade.mix_iter=sum((caudal.massico_iter(2:5)./caudal.massicotot_iter).*densidade.liquidos_iter+densidade.gas_iter); % kg/m^3
 
%viscosividade
v_etileno=-3.985+3.8726*(10^-1)*T-1.1227*(10^-4)*T^2;
v_metano=3.844+4.0112*(10^-1)*T-1.4303*(10^-4)*T^2;
v_etano=0.514+3.3449*(10^-1)*T-7.1071*(10^-5)*T^2;
visc_gas=((f(1)./caudal.molartot_iter)*v_etileno+(caudal.molar_i(6)./caudal.molartot_iter)*v_metano+(caudal.molar_i(7)./caudal.molartot_iter)*v_etano)*3.6000e-04;%kg.m/h
 
 
%calculo caudal volumétrico total
caudal.volumtot_iter=caudal.massicotot_iter/densidade.mix_iter; %m3/h
 
 
%conversão e concentrações de cada composto i
 
X=(caudal.molar_i(1)-f(1))/caudal.molar_i(1)
% para saber a conversão de etileno na reação concorrente que dá deb,
% calculou-se o valor de etileno convertido com base nos BMs vai se assumir
%uma seletividade de 95% para a reação desejada
 
Cet=caudal.molar_i(1)*(1-X)/caudal.volumtot_iter*1000; %mol/m^3
Cb=(caudal.molar_i(2)-caudal.molar_i(1)*X)/caudal.volumtot_iter*1000;
Ceb=(caudal.molar_i(3)+caudal.molar_i(2)*X)/caudal.volumtot_iter*1000;
Cdeb=(caudal.molar_i(4)+caudal.molar_i(2)*X)/caudal.volumtot_iter*1000;
%dividir por 1000 para a concentraçao entrar na equaçao de r com mol/m3
% P=f(5);
 
%
%cinetica
 
r1=0.084*exp(-9502/R.kPam3Kkmol*T)*Cet*(Cb^0.32)*3600/1000; %kmol/h.m3
r2=03603*exp(-15396/R.kPam3Kkmol*T)*(Cet^1.3)*(Cb^0.33)*3600/1000;
r3=0.00085*exp(-20643/R.kPam3Kkmol*T)*(Cet^1.77)*(Ceb^0.35)*3600/1000;
%*3600/1000 para a velocidade de reação passar para h e por kmol, de forma
%a bater certo no balanço de massa
% km1=1.672*exp(-21.735/(R.kPaLKkmol*1000*T));
% km1_inv=0.348*exp(-37.304/(R.kPaLKkmol*1000*T));
% km2=0.354*exp(-26.919/(R.kPaLKkmol*1000*T));
% km2_inv=7.3*10^-2*exp(-40.117/(R.kPaLKkmol*1000*T));
% n1=0.795;
% n2=0.110;
% n3=0.276;
% n4=0.842;
% n5=0.403;
% n6=0.548;
% r1=km1*(Cet^n1)*(Cb^n2)-km1_inv*(Ceb^n3);
% r2=km2*(Cet^n4)*(Ceb^n5)-km2_inv*(Cdeb^n6);
 
%edos para dar os caudais novamente 
edo(1)=-r1-2*r2-r3
edo(2)=-r1-r2;
edo(3)=r1-r3;
edo(4)=r2+r3;
 
gc=1;
u=caudal.massicotot_iter/cat.area;
G=densidade.gas_iter*u;
AA=G/(gc*densidade.mix_i*cat.Dp);
BB=(1-leito.porosidade)/(leito.porosidade^3);
CC=150*visc_gas*(1-leito.porosidade)/cat.Dp;
beta=AA*BB*(1.75.*G+CC); 
y=P/P_i;
coef=-1-1+1;
ezinho=coef*((f(1)+caudal.molar_i(6)+caudal.molar_i(7))/caudal.molartot_iter);
ergun=(-beta./(cat.area*(1-leito.porosidade)*cat.densidade*P_i*(f(5)))*(caudal.molartot_iter/caudal.molar_total))*(1+ezinho*X);
edo(5)=ergun-y;
end
 
If you can help me, I would be much apreciate it :) so thank you for being cool! 

 



#2 emekar

emekar

    Junior Member

  • Members
  • 14 posts

Posted 19 December 2017 - 04:49 AM

Why not just use matlab ode  solver to solve the set of coupled differential equation for material, energy and momentum balance it will be a lot shorter code and it would most likely run i have done a couple of them using matlab ode solver.



#3 Rachel Tee

Rachel Tee

    Junior Member

  • Members
  • 11 posts

Posted 23 February 2021 - 10:12 AM

Hi, is it okay if you share your code? :') I tried doing matlab also and having same issue. 


Edited by Rachel Tee, 23 February 2021 - 10:13 AM.





Similar Topics