image thumbnail

Calculate the equilibrium conversion of a packed bed reactor

by

 

Alkylation reactor producing ethylbenzene and side product di-ethylbenzene

Equilibrium_Conv
function Equilibrium_Conv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written By: George Wu Ting Wei, on 09/Nov/2012  %
% National University of Singapore {NUS}          %
%                                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the equilibrium conversion of a packed bed reactor               %
% Alkylation reactor producing ethylbenzene and side product di-ethylbenzene %
% The conversion of ethylene is of importance as is the minimization of   %
% side product formation                                                  %
% Further details can be found in Analysis, Synthesis and Design of       %
% Chemical Processes. 3rd Edition, Prentice Hall, 2009 by Turton.R,       %
% Bailie R.C., Shaeiwitz J.A                                              %
% Four reactions are considered.                                          %
% Ethylene + Benzene ? Ethylbenzene                                      %
% Ethylbenzene + Ethylene ? Diethylbenzene                               %
% Diethylbenzene + Benzene ? 2 Ethylbenzene                              %
% Toluene + Ethylene ? Ethylbenzene + Propylene                          %
% The code calculates the equilibrium conversion by minimizing Gibbs free %
% energy                                                                  %
% The notation used can be found in                                       %
% Introduction to Chemical Engineering Thermodynamics by Smith, Van Ness, %
% Abbott                                                                  %
% This code is a good example on finding equilibrium conversion when      %
% multiple reactions occur and the feed includes inerts                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear all; close all; tic
global F_ET_0 F_BZ_0 F_EB_0 F_DEB_0 F_ETN_0 F_TOL_0 F_0 F_PRP_0 T P
T_MAX=460+273.15; %Max Temperature, K
T=500; %Temperature, K, starting temperature to evaluate from
P=19.35; %Pressure Absolute, Bar
%Inlet molar flowrate per bed (mol/s)
F_ET_0=86.3*(1000/3600); %Ethylene mol/s, converted from kmol/hr
F_BZ_0=617.04*(1000/3600); %Benzene 
F_EB_0=0.74346*(1000/3600); %Ethylbenzene 
F_DEB_0=8.4097*(1000/3600); %Diethylbenzene
F_ETN_0=7.51166*(1000/3600); %Ethane 
F_TOL_0=9.91648*(1000/3600); %Toluene
F_PRP_0=7.01961*(1000/3600); %Propene
F_0=F_ET_0+F_BZ_0+F_EB_0+F_DEB_0+F_ETN_0+F_TOL_0+F_PRP_0; %Total inlet molar flowrate per bed, (mol/s)
dT=1; % Step size
Tplot=0; 
np=0; % Step counter
%Integration Loop
while T<T_MAX
    x0 = [1442; -1414; -1415; -1.9]; % Make a starting guess at the solution
    [x] = fsolve(@ModelTURTONEQUIL,x0);  % Call solver
    x(1)=x(1,1); %Rxn1 extent
    x(2)=x(2,1); %Rxn2 extent
    x(3)=x(3,1); %Rxn3 extent
    x(4)=x(4,1); %Rxn4 extent
    F_ET=(F_ET_0-x(1)-x(2)-2*x(4)); %Outlet molar flow (mol/s)
    F_BZ=(F_BZ_0-x(1)-x(3));
    F_EB=(F_EB_0+x(1)-x(2)+2*x(3)+x(4));
    F_DEB=(F_DEB_0+x(2)-x(3));
    F_TOL=(F_TOL_0-x(4));
    F_PRP=(F_PRP_0+x(4));
    F=(F_0-x(1)-x(2)-x(4));
    F_ETN=F-(F_ET+F_BZ+F_EB+F_DEB+F_TOL+F_PRP);
    y_ET=F_ET/F; %Mole fractions
    y_BZ=F_BZ/F;
    y_EB=F_EB/F;
    y_DEB=F_DEB/F;
    y_TOL=F_TOL/F;
    y_PRP=F_PRP/F;
    y_ETN=1-(y_ET+y_BZ+y_EB+y_DEB+y_TOL+y_PRP);
    ET_Cnv=(F_ET_0-F_ET)/F_ET_0;
    %Save for plotting
    if T>=Tplot;
        np=np+1; %Add parameters to next column using np as to count
        TT(np)=T-273.15; %Add temperature to T matrix
        %PP(np)=P; %Add pressure to pp matrix
        ET_Cnv1(np)=ET_Cnv*100; %Add Ethylene conversion to ET_Cnv1 matrix
        y_ET1(np)=y_ET*100; %Add composition of y_ET to y_ET1 matrix
        y_BZ1(np)=y_BZ*100; %Add composition of y_BZ to y_BZ1 matrix
        y_EB1(np)=y_EB*100; %Add composition of y_EB to y_EB1 matrix
        y_DEB1(np)=y_DEB*100; %Add composition of y_DEB to y_DEB1 matrix
        y_TOL1(np)=y_TOL*100; %Add composition of y_TOL to y_TOL1 matrix
        y_PRP1(np)=y_PRP*100; %Add composition of y_PRP to y_PRP1 matrix
        y_ETN1(np)=y_ETN*100; %Add composition of y_ETN to y_ETN1 matrix
        AA=ET_Cnv1';
        Tplot=Tplot+1;end
    %Integration
    T=T+dT;
    if P<1;break;end
end
toc %Timer
%Plot Equilibrium Conversion of Ethylene
axes1 = axes('Parent',figure,'FontWeight','bold','FontSize',20);
plot(TT,ET_Cnv1,'linewidth',3);
grid on;
ylabel('% Conversion of Ethylene');
xlabel('Temperature ^oC');
title('Equilibrium Conversion of Ethylene at various inlet temperatures');
end
 function F = ModelTURTONEQUIL(x)
global F_ET_0 F_BZ_0 F_EB_0 F_DEB_0 F_TOL_0 F_0 F_PRP_0 T P
R=8.3144621; %Gas constant
T0=298.15; %Standard temperature
deltaG1=-67235; %Rxn1 delta G at 298.15K, Standard Gibbs
deltaG2=-62650;
deltaG3=-4585;
deltaG4=-65875;
H_0_1=(29920-52510-82930); %Standard Heat of reaction of Reaction 1, J/mol
H_0_2=(-17240-29920-52510); %Standard Heat of reaction of Reaction 2
H_0_3=(2*29920+17240-82930); %Standard Heat of reaction of Reaction 3
H_0_4=(19710+29920-2*52510-50170); %Standard Heat of reaction of Reaction 4
idcph_1=@(T)(9.344936+0.46042932*T-0.000153609464*T.^2)-(11.839136+0.119671716*T-0.000036515088*T.^2)-(-1.712684+0.324778096*T-0.000110584514*T.^2); %J/mol.K
IDCPH_1 = quad(idcph_1,298.15,T); %Following the notation used in Introduction to Chemical Engineering Thermodynamics by Smith, Van Ness, Abbott
idcps_1=@(T)((9.344936+0.46042932*T-0.000153609464*T.^2)-(11.839136+0.119671716*T-0.000036515088*T.^2)-(-1.712684+0.324778096*T-0.000110584514*T.^2))./T; %J/mol.K
IDCPS_1 = quad(idcps_1,298.15,T);
Term1_1=(deltaG1-H_0_1)/(R*T0);
Term2_1=H_0_1/(R*T);
Term3_1=IDCPH_1/T;
Term4_1=IDCPS_1;
K1=exp(-(Term1_1+Term2_1+Term3_1-Term4_1));
idcph_2=@(T)(-33.93+85.5628*(T/100)-543*(T/1000).^2+135.988*(T/1000).^3)-(9.344936+0.46042932*T-0.000153609464*T.^2)-(11.839136+0.119671716*T-0.000036515088*T.^2);
IDCPH_2 = quad(idcph_2,298.15,T);
idcps_2=@(T)((-33.93+85.5628*(T/100)-543*(T/1000).^2+135.988*(T/1000).^3)-(9.344936+0.46042932*T-0.000153609464*T.^2)-(11.839136+0.119671716*T-0.000036515088*T.^2))./T;
IDCPS_2 = quad(idcps_2,298.15,T);
Term1_2=(deltaG2-H_0_2)/(R*T0);
Term2_2=H_0_2/(R*T);
Term3_2=IDCPH_2/T;
Term4_2=IDCPS_2;
K2=exp(-(Term1_2+Term2_2+Term3_2-Term4_2));
idcph_3=@(T)(2*(9.344936+0.46042932*T-0.000153609464*T.^2)-(-1.712684+0.324778096*T-0.000110584514*T.^2)-(-33.93+85.5628*(T/100)-543*(T/1000).^2+135.988*(T/1000).^3));
IDCPH_3 = quad(idcph_3,298.15,T);
idcps_3=@(T)(2*(9.344936+0.46042932*T-0.000153609464*T.^2)-(-1.712684+0.324778096*T-0.000110584514*T.^2)-(-33.93+85.5628*(T/100)-543*(T/1000).^2+135.988*(T/1000).^3))./T;
IDCPS_3 = quad(idcps_3,298.15,T);
Term1_3=(deltaG3-H_0_3)/(R*T0);
Term2_3=H_0_3/(R*T);
Term3_3=IDCPH_3/T;
Term4_3=IDCPS_3;
K3=exp(-(Term1_3+Term2_3+Term3_3-Term4_3));
idcph_4=@(T)((13.610018+0.188777684*T-0.00005749131*T.^2)+(9.344936+0.46042932*T-0.000153609464*T.^2)-2*(11.839136+0.119671716*T-0.000036515088*T.^2)-(2.41106+0.391190328*T-0.000130662824*T.^2));
IDCPH_4 = quad(idcph_4,298.15,T);
idcps_4=@(T)((13.610018+0.188777684*T-0.00005749131*T.^2)+(9.344936+0.46042932*T-0.000153609464*T.^2)-2*(11.839136+0.119671716*T-0.000036515088*T.^2)-(2.41106+0.391190328*T-0.000130662824*T.^2))./T;
IDCPS_4 = quad(idcps_4,298.15,T);
Term1_4=(deltaG4-H_0_4)/(R*T0);
Term2_4=H_0_4/(R*T);
Term3_4=IDCPH_4/T;
Term4_4=IDCPS_4;
K4=exp(-(Term1_4+Term2_4+Term3_4-Term4_4));
y_ET=(F_ET_0-x(1)-x(2)-2*x(4))/(F_0-x(1)-x(2)-x(4));
y_BZ=(F_BZ_0-x(1)-x(3))/(F_0-x(1)-x(2)-x(4));
y_EB=(F_EB_0+x(1)-x(2)+2*x(3)+x(4))/(F_0-x(1)-x(2)-x(4));
y_DEB=(F_DEB_0+x(2)-x(3))/(F_0-x(1)-x(2)-x(4));
y_TOL=(F_TOL_0-x(4))/(F_0-x(1)-x(2)-x(4));
y_PRP=(F_PRP_0+x(4))/(F_0-x(1)-x(2)-x(4));
F = [(K1*P)*(y_ET*y_BZ)-y_EB; (K2*P)*(y_EB*y_ET)-y_DEB; (K3*y_DEB*y_BZ)-y_EB^2; (K4*P)*(y_TOL*y_ET^2)-(y_EB*y_PRP)];
 end

Contact us