2 replies to this topic

[catarinajorge7]

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! 

 

[emekar]

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.

[Rachel Tee]

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.