%{
INSTITUCIÓN: Universidad Nacional de Colombia
NOMBRES: - María Paula Bravo Castillo, mbravoc@[Link]
- Juan Pablo Chitiva Arteaga, jchitiva@[Link]
- Paula Alejandra Duarte Salamanca, pduarte@[Link]
- Juan José Guerrero Acosta, juguerrero@[Link]
- Daniel Andrés López Brijaldo, dlopezbr@[Link]
ASIGNATURA: Ingeniería de las Reacciones Químicas
DOCENTE: Hugo Martín Galiendo Valbuena
TÍTULO: Proyecto Final "Diseño de Reactor" - Código en MatLab de Solución de
Ecuaciones
FECHA: 23 de junio de 2023
1) RESUMEN:
Se presenta a continuación el código elaborado en MatLab para el desarrollo y
modelamiento de un reactor tubular de lecho empacado. El código presenta los
siguientes componentes:
- Establecimiento de parámetros de ingreso y cualidades constantes.
- Cálculo de propiedades empleando correcciones con la ecuación de gases
ideales.
- Cálculo de propiedades empleando correlaciones en función de
temperatura.
- Balances de materia, energía y presión en el sistema.
- Establecimiento del ciclo de resolución de ecuaciones diferenciales con
el Solver ODE45 y construyendo una función en MatLab para realizar los
cálculos repetitivos.
- Gráfica de perfiles de concentraciones, temperaturas y presión.
Se buscó que todas las propiedades estuvieran en unidades del sistema
internacional para que los cálculos tuvieran consistencia dimensional.
Así mismo, la nomenclatura empleada para identificar las sustancias químicas a
utilizar es la siguiente:
- A: Metanol (CH3OH)
- B: Oxígeno (O2)
- C: Formaldehído (CH2O)
- D: Agua (H2O)
- E: Monóxido de Carbono (CO)
- F: Nitrógeno (N2)
- G: Fluido de Control
%}
% 2) DECLARACIÓN DE VALORES CONSTANTES E INICIALES (Pre - Función):
% Pesos moleculares [kg / mol]: Valores consultados y convertidos a kg.
PM_A = 32.04 / 1000; PM_B = 32.00 / 1000; PM_C = 30.03 / 1000; PM_D =
18.00 / 1000; PM_E = 28.01 / 1000; PM_F = 28.01 / 1000;
1
% Parámetros cinéticos obtenidos con las expresiones de Cinética B
(Denominador más corto), [Eact = J/mol]:
A_k1 = 121.5539; Eact_k1 = 62661.3851;
A_K1A = 62661.3851; Eact_K1A = 133441.5273;
A_k2 = 0.017971; Eact_k2 = 32820.4093;
A_K2D = 1.1538; Eact_K2D = -16817.5167;
% Parámetros de ingreso de fluido de reacción (En x = 0):
% Condiciones generales iniciales para fluido de reacción:
T_rin = 573; % [K]
PR_in = 2.5 * 101.325 * 1000; % [Pa]
ConstR = 8.31446261815324; % [J / mol K] o [Pa * m3 / mol K]
% Concentraciones molares iniciales:
y_Ain = 0.2; y_Cin = 0.001; y_Din = 0.02; y_Ein = 0.008;
y_Bin = (1 - y_Ain - y_Cin - y_Din - y_Ein) * 0.21; y_Fin = (1 - y_Ain -
y_Cin - y_Din - y_Ein) * 0.79;
% Flujos de reacción másicos (fr [kg / s]) y molares (FR [mol / s])
iniciales:
fR_in = 0.6;
FR_in = fR_in / (y_Ain * PM_A + y_Bin * PM_B + y_Cin * PM_C + y_Din *
PM_D + y_Ein * PM_E + y_Fin * PM_F);
FR_Ain = FR_in * y_Ain; FR_Bin = FR_in * y_Bin; FR_Cin = FR_in * y_Cin;
FR_Din = FR_in * y_Din; FR_Ein = FR_in * y_Ein;
FR_Fin = FR_in * y_Fin;
fr_Ain = FR_Ain * PM_A; fr_Bin = FR_Bin * PM_B; fr_Cin = FR_Cin * PM_C;
fr_Din = FR_Din * PM_D; fr_Ein = FR_Ein * PM_E; fr_Fin = FR_Fin * PM_F;
% Caudales (Q) iniciales [m3]: Se utiliza la ecuación de gas ideal con
los parámetros de entrada.
Q_in = ConstR * FR_in * T_rin / (PR_in);
Q_Ain = Q_in * y_Ain; Q_Bin = Q_in * y_Bin; Q_Cin = Q_in * y_Cin;
Q_Din = Q_in * y_Din; Q_Ein = Q_in * y_Ein; Q_Fin = Q_in * y_Fin;
% Concentraciones molares (CR) [mol / m3] iniciales:
CR_in = FR_in / Q_in;
CR_Ain = CR_in* y_Ain; CR_Bin = CR_in * y_Bin; CR_Cin = CR_in * y_Cin;
CR_Din = CR_in * y_Din; CR_Ein = CR_in * y_Ein; CR_Fin = CR_in * y_Fin;
% Parámetros de Fluido de Control:
% Referencia: HP Hytherm 600 ([Link]
specialties/thermic-fluids/hytherm-500-and-600-thermic-fluid-oil)
f_G = 0.2; % [kg / s];
T_cin = 300; % [K]
CondTerm_G = 0.090 / 1.1629999998093; % [W / m K]
Mu_G = 0.0279; % [kg / m s]
% Características de Diseño de Reactor: Aspectos físicos brindados y/o
elegidos para el reactor.
di = 2.0574 / 100; % [m]
2
A_transi = pi() * di^2 / 4; % [m2]
do = 2.54 / 100; % [m]
d = do - di; % [m]
MLDd = (do - di) / log(do / di); % [m]
Censu = 0.01 / 5.67826334106816; % [W / m2 K]
CondTerm_Tubo = 16.20; % [W / m2 K]
% Características del Sistema del Lecho: Incluye propiedades asociadas a las
partículas catalizadoras y al lecho:
dp = 0.4 / 100; % [m]
Ro_p = 1.892 * 1000; % [kg / m3]
ap = (6 / dp) * 1/Ro_p; % [m2 particula / Kg cat]
Vacio = 0.390 + (1.740 / ((di/dp) + 1.140)^2);
nR_G = 4.5; % [kg / m2 s]
f_s = nR_G * A_transi; % [kg / s]
T_sin = 300; % [K]
Ro_Lecho = (1 - Vacio) * Ro_p; % [kg / m3]
% 3) ESTABLECIMIENTO RESOLUCIÓN DE SISTEMA DE ECUACIONES CON SOLVER ODE45:
% Configuración general:
y = [0, 0, 0, 0, 0];
a = 48; % 24 48 96 192
Long = 0.6;
Lmax = 2.5; % [m]
tspan = linspace(0, Lmax, 1000); % Intervalo de evaluación
% CI = [FR_Ain; FR_Cin; T_rin; T_cin; PR_in]; % Condiciones Iniciales
para CasoRelEsteq1 (Metanol y Formaldehído)
CI = [FR_Ain; FR_Ein; T_rin; T_cin; PR_in]; % Condiciones Iniciales para
CasoRelEsteq2 (Metanol y Monóxido de Carbono)
y = CI;
M = [1 0 0 0 0; 0 1 0 0 0 ; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1]; % Matriz de
masa
options = odeset('RelTol', 1e-2, 'AbsTol', [1e-2 1e-2 1e-2 1e-2 1e-2]); %
Configuración propuesta para el solver
[J Y] = ode45(@Derivadas, tspan, CI); % Ejecución del Solver
% 4) MANIPULACIÓN DE VECTORES RESULTANTES Y ELABORACIÓN DE GRÁFICAS:
CasoMani = 2;
switch(CasoMani)
case 1 % Empleo de Metanol (A) y Formaldehido (C)
for i = 1 : length(Y')
VFR_A(i) = Y(i, 1); % Flujo de Metanol
VFR_C(i) = Y(i, 2); % Flujo de Formaldehído
VE1(i) = (FR_Ain - VFR_A(i)); % Avance de Reacción 1 en
función de Metanol (A)
VE2(i) = (FR_Cin - VFR_C(i) + VE1(i)); % Avance de
Reacción 2 en función de Formaldehído (C)
VFR_B(i) = (FR_Bin - 0.5 * VE1(i) - 0.5 * VE2(i)); % Flujo
Molar de Oxígeno (B)
VFR_D(i) = (FR_Din + VE1(i) + VE2(i)); % Flujo Molar de
Agua (D)
VFR_E(i) = (FR_Ein + VE2(i)); % Flujo Molar de Monóxido de
Carbono (C)
3
VFR_F(i) = FR_Fin; % Flujo Molar de Nitrógeno (F)
end
case 2 % Empleo de Metanol (A) y Monóxido de Carobono (E)
for i = 1 : length(Y') % Valores obtenidos en la Resolución de
Ecuaciones
VFRA(i) = Y(i, 1); % Flujo de Metanol
VFRE(i) = Y(i, 2); % Flujo de CO
end
for i = 1 : (length(Y') - 1) % Moles que reaccionaron
VE1(i) = (VFRA(i) - VFRA(i + 1)) * a; % Diferencial de
Moles de Reacción 1
VE2(i) = (VFRE(i + 1) - VFRE(i)) * a; % Diferencial de
Moles de Reaccion 2
X(i) = 0.5 * (J(i) + J(i + 1)); % Valores Promedio
end
% Creación de Vectores de Concentración y Planteamiento de sus
Valores Iniciales:
VFR_A = zeros(length(Y') - 1, 1); VFR_B = zeros(length(Y') - 1,
1); VFR_C = zeros(length(Y') - 1, 1);
VFR_D = zeros(length(Y') - 1, 1); VFR_E = zeros(length(Y') - 1,
1); VFR_F = zeros(length(Y') - 1, 1); VFR_R = zeros(length(Y') - 1, 1);
VFR_A(1) = FR_Ain; VFR_B(1) = FR_Bin; VFR_C(1) = FR_Cin;
VFR_D(1) = FR_Din; VFR_E(1) = FR_Ein; VFR_F(1) = FR_Fin; VFR_R(1) = FR_in;
% Variación de Concentraciones:
for i = 2 : (length(Y') - 1)
VFR_A(i) = VFR_A(i - 1) - VE1(i);
VFR_B(i) = VFR_B(i - 1) - 0.5 * VE1(i) - 0.5 * VE2(i);
VFR_C(i) = VFR_C(i - 1) + VE1(i) - VE2(i);
VFR_D(i) = VFR_D(i - 1) + VE1(i) + VE2(i);
VFR_E(i) = VFR_E(i - 1) + VE2(i);
VFR_F(i) = VFR_F(i - 1);
VFR_R(i) = VFR_A(i) + VFR_B(i) + VFR_C(i) + VFR_D(i) +
VFR_E(i) + VFR_F(i);
end
% Vectores de Concentraciones, Temperaturas y Presión:
for i = 1 : length(Y') - 1
VT_r(i) = Y(i, 3); % Temperatura de R
VT_c(i) = Y(i, 4); % Temperatura de C
VPR(i) = Y(i, 5) / 1000; % Presión del FLuido de Reaccion
VQ(i) = Q_in * (VFR_R(i) / FR_in) * (VT_r(i) / T_rin) *
(PR_in / VPR(i)) / 1000;
VConv_A(i) = (FR_Ain - VFR_A(i)) / FR_Ain;
VRend_C(i) = (VE1(i) - VE2(i)) / (VE1(i));
VConc_A(i) = VFR_A(i) / VQ(i);
VConc_B(i) = VFR_B(i) / VQ(i);
VConc_C(i) = VFR_C(i) / VQ(i);
4
VConc_D(i) = VFR_D(i) / VQ(i);
VConc_E(i) = VFR_E(i) / VQ(i);
VConc_F(i) = VFR_F(i) / VQ(i);
VConc_R(i) = VFR_R(i) / VQ(i);
end
end
% Construcción de gráficas:
subplot(2, 3, 1);
plot(X, VConv_A, '-', 'linewidth', 2);
hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Conversión")
title("Conversión de Metanol")
grid on;
axis([0 Long 0 1])
subplot(2, 3, 2);
plot(X, VRend_C, '-', 'linewidth', 2);
hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Rendimiento")
title("Rendimiento de Formaldehído")
grid on;
axis([0 Long min(VRend_C) * 0.999 max(VRend_C) * 1.001])
subplot(2, 3, 3);
plot(X, VFR_A, '-', 'linewidth', 2); hold on;
plot(X, VFR_B, '-', 'linewidth', 2); hold on;
plot(X, VFR_C, '-', 'linewidth', 2); hold on;
plot(X, VFR_D, '-', 'linewidth', 2); hold on;
plot(X, VFR_E, '-', 'linewidth', 2); hold on;
plot(X, VFR_F, '-', 'linewidth', 2); hold on;
plot(X, VFR_R, '-', 'linewidth', 2); hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Flujo molar [mol /s]")
title("Flujos Molares de Reacción")
legend("Metanol", "Oxígeno", "Formaldehído", "Agua", "Monóxido de
Carbono", "Nitrógeno", "Flujo Total")
axis([0.01 Long 0 VFR_R(length(VFR_R)) * 1.2])
grid on;
subplot(2, 3, 4);
plot(X, VConc_A, '-', 'linewidth', 2); hold on;
plot(X, VConc_B, '-', 'linewidth', 2); hold on;
plot(X, VConc_C, '-', 'linewidth', 2); hold on;
plot(X, VConc_D, '-', 'linewidth', 2); hold on;
plot(X, VConc_E, '-', 'linewidth', 2); hold on;
plot(X, VConc_F, '-', 'linewidth', 2); hold on;
plot(X, VConc_R, '-', 'linewidth', 2); hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Concentracion molar [mol / m3]")
title("Concentraciones de Flujo de Reacción")
5
legend("Metanol", "Oxígeno", "Formaldehído", "Agua", "Monóxido de
Carbono", "Nitrógeno", "Flujo Total")
axis([0.01 Long 0 VConc_R(1) * 1.2])
grid on;
subplot(2, 3, 5)
plot(X * Long / Lmax, VT_r, 'linewidth', 2); hold on;
plot(X * Long / Lmax, VT_c, 'linewidth', 2); hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Temperatura [K]")
title("Temperaturas de Flujos de Reacción y de Control")
legend("Fluido de Reacción", "Fluido de Control")
axis([0.01 Long T_cin T_rin])
grid on;
subplot(2, 3, 6)
plot(X, VPR, 'linewidth', 2); hold on;
xlabel("Distancia de Reactor [m]")
ylabel("Presión [kpa]")
title("Presión en las Tuberías de Flujo de Reacción")
grid on;
% 5. PLANTEAMIENTO DE FUNCIÓN PARA SOLUCIÓN DE ECUACIONES:
function Ec = Derivadas(x, y)
% 5.A) DECLARACIÓN DE VALORES CONSTANTES E INICIALES (Dentro de la Función):
Se realiza esto con el fin de no requerir el llamado de todos los valores
% anteriormente utilizados como constantes.
% Pesos moleculares [kg / mol]: Valores consultados y convertidos a kg.
PM_A = 32.04 / 1000; PM_B = 32.00 / 1000; PM_C = 30.03 / 1000; PM_D =
18.00 / 1000; PM_E = 28.01 / 1000; PM_F = 28.01 / 1000;
% Parámetros cinéticos obtenidos con las expresiones de Cinética B
(Denominador más corto), [Eact = J/mol]:
A_k1 = 121.5539; Eact_k1 = 62661.3851;
A_K1A = 62661.3851; Eact_K1A = 133441.5273;
A_k2 = 0.017971; Eact_k2 = 32820.4093;
A_K2D = 1.1538; Eact_K2D = -16817.5167;
% Parámetros de ingreso de fluido de reacción (En x = 0):
% Condiciones generales iniciales para fluido de reacción:
T_rin = 573; % [K]
PR_in = 2.5 * 101.325 * 1000; % [Pa]
ConstR = 8.31446261815324; % [J / mol K] o [Pa * m3 / mol K]
% Concentraciones molares iniciales:
y_Ain = 0.2; y_Cin = 0.001; y_Din = 0.02; y_Ein = 0.008; y_Bin = (1 -
y_Ain - y_Cin - y_Din - y_Ein) * 0.21; y_Fin = (1 - y_Ain - y_Cin - y_Din -
y_Ein) * 0.79;
% Flujos de reacción másicos (fr [kg / s]) y molares (FR [mol / s])
iniciales:
fR_in = 0.6;
6
FR_in = fR_in / (y_Ain * PM_A + y_Bin * PM_B + y_Cin * PM_C + y_Din *
PM_D + y_Ein * PM_E + y_Fin * PM_F);
FR_Ain = FR_in * y_Ain; FR_Bin = FR_in * y_Bin; FR_Cin = FR_in * y_Cin;
FR_Din = FR_in * y_Din; FR_Ein = FR_in * y_Ein; FR_Fin = FR_in * y_Fin;
fr_Ain = FR_Ain * PM_A; fr_Bin = FR_Bin * PM_B; fr_Cin = FR_Cin * PM_C;
fr_Din = FR_Din * PM_D; fr_Ein = FR_Ein * PM_E; fr_Fin = FR_Fin * PM_F;
% Caudales (Q) iniciales [m3]: Se utiliza la ecuación de gas ideal con
los parámetros de entrada.
Q_in = ConstR * FR_in * T_rin / (PR_in);
Q_Ain = Q_in * y_Ain; Q_Bin = Q_in * y_Bin; Q_Cin = Q_in * y_Cin; Q_Din =
Q_in * y_Din; Q_Ein = Q_in * y_Ein; Q_Fin = Q_in * y_Fin;
% Concentraciones molares (CR) [mol / m3] iniciales:
CR_in = FR_in / Q_in;
CR_Ain = CR_in* y_Ain; CR_Bin = CR_in * y_Bin; CR_Cin = CR_in * y_Cin;
CR_Din = CR_in * y_Din; CR_Ein = CR_in * y_Ein; CR_Fin = CR_in * y_Fin;
% Parámetros de Fluido de Control:
% Referencia: HP Hytherm 600 ([Link]
specialties/thermic-fluids/hytherm-500-and-600-thermic-fluid-oil)
f_G = 0.2; % [kg / s];
T_cin = 300; % [K]
CondTerm_G = 0.090 / 1.1629999998093; % [W / m K]
Mu_G = 0.0279; % [kg / m s]
% Características de Diseño de Reactor: Aspectos físicos brindados y/o
elegidos para el reactor.
di = 2.0574 / 100; % [m]
A_transi = pi() * di^2 * 0.25; % [m2]
do = 2.54 / 100; % [m]
d = do - di; % [m]
MLDd = (do - di) / log(do / di); % [m]
Censu = 0.001 * 5.67826334106816; % [m2 K / W]
CondTerm_Tubo = 16.20; % [W / m2 K]
% Características del Sistema del Lecho: Incluye propiedades asociadas a las
partículas catalizadoras y al lecho:
dp = 0.4 / 100; % [m]
Ro_p = 1.892 * 1000; % [kg / m3]
ap = (6 / dp) * (1 / Ro_p); % [m2 particula / Kg cat]
Vacio = 0.390 + (1.740 / ((di/dp) + 1.140)^2);
nR_Rsup = 4.5; % [kg / m2 s] % Flux másico superficial
f_Rsup = nR_Rsup * A_transi; % [kg / s] % Flujo másico Superficial
T_sin = 600; % [K]
Ro_Lecho = (1 - Vacio) * Ro_p; % [kg / m3]
% 5.B ESTABLECIMIENTO DE VARIABLES:
% Variables Dependientes:
% - y(1): Flujo Molar 1
% - y(2): Flijo Molar 2
% - y(3): Tr
7
% - y(4): Tc
% - y(5): PR
% y = [FR_Ain; FR_Ein; T_rin; T_cin; PR_in];
T_r = y(3);
T_c = y(4);
PR = y(5);
T_s = T_sin;
% Relaciones Estequiométricas y Flujos Molares (FR) [mol / s]: Se plantean dos
alternativas: Emplear como variables independientes 1) Metanol y Formaldehído
% 2) Metanol y Monóxido de Carbono:
% A: Metanol (CH3OH), B: Oxígeno (O2), C: Formaldehído (CH2O), D: Agua
(H2O), E: Monóxido de Carbono (CO), F: Nitrógeno (N2)
CasoRelEstq = 2;
switch(CasoRelEstq)
case 1
FR_A = y(1); % Variable 1: Flujo Molar de Metanol (A)
FR_C = y(2); % Variable 2: Flujo Molar de Formaldehído (C)
E1 = FR_Ain - FR_A; % Avance de Reacción 1 en función de
Metanol (A)
E2 = FR_Cin - FR_C + E1; % Avance de Reacción 2 en función de
Formaldehído (C)
FR_B = FR_Bin - 0.5 * E1 - 0.5 * E2; % Flujo Molar de Oxígeno
(B)
FR_D = FR_Din + E1 + E2; % Flujo Molar de Agua (D)
FR_E = FR_Ein + E2; % Flujo Molar de Monóxido de Carbono (C)
FR_F = FR_Fin; % Flujo Molar de Nitrógeno (F)
FR = FR_A + FR_B + FR_C + FR_D + FR_E + FR_F; % Flujo Molar
Total
case 2
FR_A = y(1); % Variable 1: Flujo Molar de Metanol (A)
FR_E = y(2); % Variable 2: Flujo Molar de Monóxido de Carbono
(C)
E1 = FR_Ain - FR_A; % Avance de Reacción 1 en función de
Metanol (A)
E2 = FR_E - FR_Ein; % Avance de Reacción 2 en función de
Monóxido de Carbono (E)
FR_B = FR_Bin - 0.5 * E1 - 0.5 * E2; % Flujo Molar de Oxígeno
(B)
FR_C = FR_Cin + E1 - E2; % Flujo Molar de Formaldehído (C)
FR_D = FR_Din + E1 + E2; % Flujo Molar de Agua (D)
FR_F = FR_Fin; % Flujo Molar de Nitrógeno (F)
FR = FR_A + FR_B + FR_C + FR_D + FR_E + FR_F; % Flujo Molar
Total
end
% Flujos másicos de reacción (fR) [kg / s]: Producto entre flujo molar y peso
molecular.
fR_A = FR_A * PM_A; fR_B = FR_B * PM_B; fR_C = FR_C * PM_C; fR_D = FR_D *
PM_D; fR_E = FR_E * PM_E; fR_F = FR_F * PM_F;
fR = fR_A + fR_B + fR_C + fR_D + fR_E + fR_F;
8
% Fluxes Molares (NR) [kg / m2 s] y Másicos de reacción (NR) [kg / m2 s]:
Cociente entre flujos molares o másicos con el área transversal interna.
NR_A = FR_A / A_transi; NR_B = FR_B / A_transi; NR_C = FR_C / A_transi;
NR_D = FR_D / A_transi; NR_E = FR_E / A_transi; NR_F = FR_F / A_transi;
NR_R = NR_A + NR_B + NR_C + NR_D + NR_E + NR_F;
nR_A = fR_A / A_transi; nR_B = fR_B / A_transi; nR_C = fR_C / A_transi;
nR_D = fR_D / A_transi; nR_E = fR_E / A_transi; nR_F = fR_F / A_transi;
nR_R = nR_A + nR_B + nR_C + nR_D + nR_E + nR_F;
% Concentraciones másicas de reacción (x) [kg i / kg total] y molares de
reacción (y): Cocicientes entre flujos molares y másicos parciales con los
flujos
% totales.
x_A = fR_A / fR; x_B = fR_B / fR; x_C = fR_C / fR; x_D = fR_D / fR; x_E =
fR_E / fR; x_F = fR_F / fR;
y_A = FR_A / FR; y_B = FR_B / FR; y_C = FR_C / FR; y_D = FR_D / FR; y_E =
FR_E / FR; y_F = FR_F / FR;
PM_R = y_A * PM_A + y_B * PM_B + y_C * PM_C + y_D * PM_D + y_E * PM_E +
y_F * PM_F; % Peso molecular promedio [kg / mol].
% Caudales de reacción (Q) [m3 / s]: Cálculo de caudales empleando una
corrección de gases ideales (Las variables FR, T_r y PR son de iteración, por
lo que
% continuamente corregirán el caudal), así como el cálculo de caudales
parciales con las concentraciones molares:
% (PR / PR_in) * (Q / Q_in) = (FR / FR_in) * (T_r / T_rin)
Q = Q_in * (FR / FR_in) * (T_r / T_rin) * (PR_in / PR);
Q_A = Q * y_A; Q_B = Q * y_B; Q_C = Q * y_C; Q_D = Q * y_D; Q_E = Q *
y_E; Q_F = Q * y_F;
Vel_R = Q / A_transi; % Velocidad de Flujo de Reacción [m /s]:
% Presiones parciales de reacción (PR) [Pa]: Cálculo de presiones parciales
con la presión total y la concentración molar parcial.
PR_A = PR * y_A; PR_B = PR * y_B; PR_C = PR * y_C; PR_D = PR * y_D; PR_E
= PR * y_E; PR_F = PR * y_F;
% Densidades de reacción (Ro) [kg / m3]: Cocientes entre los flujos másicos y
los caudales volumétricos parciales de cada especie. La densidad promedio
% se cálcula con la fórmula de inversos:
Ro_A = fR_A / Q_A; Ro_B = fR_B / Q_B; Ro_C = fR_C / Q_C; Ro_D = fR_D /
Q_D; Ro_E = fR_E / Q_E; Ro_F = fR_F / Q_F;
Ro_R = 1 / (x_A / Ro_A + x_B / Ro_B + x_C / Ro_C + x_D / Ro_D + x_E /
Ro_E + x_F / Ro_F);
Vel_Rsup = nR_Rsup / Ro_R;
% Concentraciones molares [mol / m3]: Cociente entre los flujos molares y el
caudal total del flujo de reacción:
9
CR_A = FR_A / Q; CR_B = FR_B / Q; CR_C = FR_C / Q; CR_D = FR_D / Q; CR_E
= FR_E / Q; CR_F = FR_F / Q;
% Viscosidades de Flujo de Reacción [Pa * s]: Se emplean correlaciones
en función de la temperatura de reacción, las cuales estarán cambiando
constantemente.
% Sus valores bajos se deben a que las sustancias se encuentran en
fase gaseosa. Se calculó su valor promedio empleando la fórmula de valores
inversos.
Mu_A = 3 * 10^-8 * T_r + 5 * 10^-7;
Mu_B = 2 * 10^-9 * T_r^2 - 2 * 10-7 * T_r + 1 * 10^-5;
Mu_C = 3 * 10^-8 * T_r + 5 * 10^-7;
Mu_D = 4 * 10^-9 * T_r - 2 * 10^-6;
Mu_E = 1 * 10^-9 * T_r^2 - 2 * 10^-7 * T_r + 9 * 10^-6;
Mu_F = 2 * 10^-9 * T_r^2 - 2 * 10^-7 * T_r + 1 * 10^-5;
Mu_R = 1 / (x_A / Mu_A + x_B / Mu_B + x_C/ Mu_C + x_D / Mu_D + x_E / Mu_E
+ x_F / Mu_F);
% Conductividades Térmicas de Flujo de Reacción [J / m K s]: Se emplean
correlaciones en función de la temperatura de reacción, las cuales estarán
cambiando
% constantemente. Se calculó su valor promedio empleando la fórmula de
valores inversos.
CondTerm_A = 5 * 10^-7 * T_r^2 - 0.0002 * T_r + 0.0284;
CondTerm_B = 1 * 10^-5 * T_r^2 - 0.0018 * T_r + 0.0758;
CondTerm_C = 5 * 10^-7 * T_r^2 - 0.0002 * T_r + 0.0284;
CondTerm_D = 1 * 10^-6 * T_r^2 - 0.0011 * T_r + 0.2011;
CondTerm_E = 9 * 10^-6 * T_r^2 - 0.0014 * T_r + 0.0588;
CondTerm_F = 7 * 10^-6 * T_r^2 - 0.0010 * T_r + 0.0462;
CondTerm_R = 1 / (x_A / CondTerm_A + x_B / CondTerm_B + x_C / CondTerm_C
+ x_D / CondTerm_D + x_E / CondTerm_E + x_F / CondTerm_F);
% Régimen Hidrodinámico y Factor de Caída de Presión: El número de Reynolds se
calcula en torno a la partícula. Luego se calcula el factor de fricción.
Re_R = (dp * nR_R / Mu_R); % [m] * [kg / m2 s] / [kg / m s] = Ad
f = (1 - Vacio) / (Vacio^3) * (1.75 + (150 * (1 - Vacio)) / Re_R);
% Coeficientes de Capacidad Calorífica, Calores Específicos y Entalpías de
Reacción: Dos fuentes de información serán empleadas para los valores de las
% especies en el fluido de reacción. Para el fluido de control se utilizó
una correlación empleando su información dada:
CasoCapCal = 2;
switch(CasoCapCal)
case 1 % Obtenidos de "Principios Elementales de los Procesos
Químicos - Felder" - Tabla B.2 (Todos menos G)
% Coeficientes: Para ser manejados en grados Celsius, en unidades de
[KJ / mol K] y de la forma "C1 + C2 * T + C3 * T^2 + C4 * T^3":
C1_A = (+42.93 * 10^-3); C2_A = (+8.3010 * 10^-5); C3_A = (-1.8700 *
10^-8); C4_A = (-8.030 * 10^-12);
C1_B = (+29.10 * 10^-3); C2_B = (+1.1580 * 10^-5); C3_B = (-0.6076 *
10^-8); C4_B = (+1.311 * 10^-12);
C1_C = (+34.28 * 10^-3); C2_C = (+4.2680 * 10^-5); C3_C = (+0.0000 *
10^-8); C4_C = (-8.694 * 10^-12);
10
C1_D = (+33.46 * 10^-3); C2_D = (+0.6880 * 10^-5); C3_D = (+0.7604 *
10^-8); C4_D = (-3.593 * 10^-12);
C1_E = (+28.95 * 10^-3); C2_E = (+0.4110 * 10^-5); C3_E = (+0.3548 *
10^-8); C4_E = (-2.220 * 10^-12);
C1_F = (+29.00 * 10^-3); C2_F = (+0.2199 * 10^-5); C3_F = (+0.5723 *
10^-8); C4_F = (-2.871 * 10^-12);
% Calores Específicos y deltas CP: Corregidos para ser manejados en
[J /mol K]:
Cp_A = (C1_A + C2_A * (T_r - 273.15) + C3_A * (T_r - 273.15)^2 +
C4_A * (T_r - 273.15)^3) * 1000;
Cp_B = (C1_B + C2_B * (T_r - 273.15) + C3_B * (T_r - 273.15)^2 +
C4_B * (T_r - 273.15)^3) * 1000;
Cp_C = (C1_C + C2_C * (T_r - 273.15) + C3_C * (T_r - 273.15)^2 +
C4_C * (T_r - 273.15)^3) * 1000;
Cp_D = (C1_D + C2_D * (T_r - 273.15) + C3_D * (T_r - 273.15)^2 +
C4_D * (T_r - 273.15)^3) * 1000;
Cp_E = (C1_E + C2_E * (T_r - 273.15) + C3_E * (T_r - 273.15)^2 +
C4_E * (T_r - 273.15)^3) * 1000;
Cp_F = (C1_F + C2_F * (T_r - 273.15) + C3_F * (T_r - 273.15)^2 +
C4_F * (T_r - 273.15)^3) * 1000;
DCP_1 = (C1_D + C1_C - 0.5 * C1_B - C1_A) * (T_r - 273.15)^0 + (C2_D
+ C2_C - 0.5 * C2_B - C2_A) * (T_r - 273.15)^1 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A) * (T_r - 273.15)^2 + (C4_D
+ C4_C - 0.5 * C4_B - C4_A) * (T_r - 273.15)^3;
DCP_2 = (C1_E + C1_D - 0.5 * C1_B - C1_C) * (T_r - 273.15)^0 + (C2_E
+ C2_D - 0.5 * C2_B - C2_C) * (T_r - 273.15)^1 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C) * (T_r - 273.15)^2 + (C4_E
+ C4_D - 0.5 * C4_B - C4_C) * (T_r - 273.15)^3;
% Entalpía de Reacción: Se indica el valor estándar a 25 C, y la
integral entre temperaturas con la antiderivada desarrollada y
% evaluada, teniendo unidades de [J / mol]:
EntR1 = ((-115.90 - 241.83) - (-201.20 + 0.5 * 0)) * 1000 + ...
(((C1_D + C1_C - 0.5 * C1_B - C1_A)/1 * (T_r - 273.15)^1 +
(C2_D + C2_C - 0.5 * C2_B - C2_A)/2 * (T_r - 273.15)^2 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A)/3 * (T_r - 273.15)^3 +
(C4_D + C4_C - 0.5 * C4_B - C4_A)/4 * (T_r - 273.15)^4) - ...
((C1_D + C1_C - 0.5 * C1_B - C1_A)/1 * (25)^1 + (C2_D +
C2_C - 0.5 * C2_B - C2_A)/2 * (25)^2 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A)/3 * (25)^3 + (C4_D +
C4_C - 0.5 * C4_B - C4_A)/4 * (25)^4));
EntR2 = ((-110.52 - 241.83) - (-115.90 + 0.5 * 0)) * 1000 + ...
(((C1_E + C1_D - 0.5 * C1_B - C1_C)/1 + (T_r - 273.15)^1 +
(C2_E + C2_D - 0.5 * C2_B - C2_C)/2 * (T_r - 273.15)^2 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C)/3 * (T_r - 273.15)^3 +
(C4_E + C4_D - 0.5 * C4_B - C4_C)/4 * (T_r - 273.15)^4) - ...
((C1_E + C1_D - 0.5 * C1_B - C1_C)/1 + (25)^1 + (C2_E +
C2_D - 0.5 * C2_B - C2_C)/2 * (25)^2 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C)/3 * (25)^3 + (C4_E +
C4_D - 0.5 * C4_B - C4_C)/4 * (25)^4));
11
case 2 % Obtenidos de la base de datos de NIST
% Coeficientes: Para ser manejados en [J / mol K] y de la forma "C1
+ C2 * T + C3 * T^2 + C4 * T^3 + C5 / T^-2":
C1_A = (+26.84900); C2_A = (+0.075000); C3_A = (-0.000010); C4_A =
(0); C5_A = (0); % [J / mol K]
C1_B = (+31.32234); C2_B = (-20.23531); C3_B = (+57.86644); C4_B =
(-36.50624); C5_B = (-0.007374); %
C1_C = (+5.193767); C2_C = (+93.23249); C3_C = (-44.85457); C4_C =
(+7.882279); C5_C = (+0.551175); % [J / mol K]
C1_D = (+30.09200); C2_D = (+6.832514); C3_D = (+6.793435); C4_D =
(-2.534480); C5_D = (+0.082139); % [J / mol K]
C1_E = (+25.56759); C2_E = (+6.096130); C3_E = (+4.054656); C4_E =
(-2.671021); C5_E = (+0.131021);% [J / mol K]
C1_F = (+19.50583); C2_F = (+19.88705); C3_F = (-8.598535); C4_F =
(+1.369784); C5_F = (+0.527601); % [J / mol K]
% Calores Específicos y deltas CP: Corregidos para ser manejados en
[J /mol K]:
Cp_A = (C1_A + C2_A * (T_r / 1000) + C3_A * (T_r / 1000)^2 + C4_A *
(T_r / 1000)^3 + C5_A * (T_r / 1000)^-2);
Cp_B = (C1_B + C2_B * (T_r / 1000) + C3_B * (T_r / 1000)^2 + C4_B *
(T_r / 1000)^3 + C5_B * (T_r / 1000)^-2);
Cp_C = (C1_C + C2_C * (T_r / 1000) + C3_C * (T_r / 1000)^2 + C4_C *
(T_r / 1000)^3 + C5_C * (T_r / 1000)^-2);
Cp_D = (C1_D + C2_D * (T_r / 1000) + C3_D * (T_r / 1000)^2 + C4_D *
(T_r / 1000)^3 + C5_D * (T_r / 1000)^-2);
Cp_E = (C1_E + C2_E * (T_r / 1000) + C3_E * (T_r / 1000)^2 + C4_E *
(T_r / 1000)^3 + C5_E * (T_r / 1000)^-2);
Cp_F = (C1_F + C2_F * (T_r / 1000) + C3_F * (T_r / 1000)^2 + C4_F *
(T_r / 1000)^3 + C5_F * (T_r / 1000)^-2);
DCP_1 = (C1_D + C1_C - 0.5 * C1_B - C1_A) * (T_r / 1000)^0 + (C2_D +
C2_C - 0.5 * C2_B - C2_A) * (T_r / 1000)^1 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A) * (T_r / 1000)^2 + (C4_D +
C4_C - 0.5 * C4_B - C4_A) * (T_r / 1000)^3 + ...
(C5_D + C5_C - 0.5 * C5_B - C5_A) * (T_r / 1000)^-2;
DCP_2 = (C1_E + C1_D - 0.5 * C1_B - C1_C) + (T_r / 1000)^0 + (C2_E +
C2_D - 0.5 * C2_B - C2_C) * (T_r / 1000)^1 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C) * (T_r / 1000)^2 + (C4_E +
C4_D - 0.5 * C4_B - C4_C) * (T_r / 1000)^3 + ...
(C5_E + C5_D - 0.5 * C5_B - C5_C) * (T_r / 1000)^-2;
% Entalpía de Reacción: Se indica el valor estándar a 25 C, y la
integral entre temperaturas con la antiderivada desarrollada y
% evaluada, teniendo presente la necesidad de multiplicar por
1000 debido a la sustitución que se realizó en el desarrollo de la integral,
% y teniendo unidades de [J / mol]:
EntR1 = ((-115.90 - 241.83) - (-201.20 + 0.5 * 0)) * 1000 + ...
(((C1_D + C1_C - 0.5 * C1_B - C1_A)/1 * (T_r / 1000)^1 +
(C2_D + C2_C - 0.5 * C2_B - C2_A)/2 * (T_r / 1000)^2 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A)/3 * (T_r / 1000)^3 +
(C4_D + C4_C - 0.5 * C4_B - C4_A)/4 * (T_r / 1000)^4 + ...
(C5_D + C5_C - 0.5 * C5_B - C5_A)/-1 * (T_r / 1000)^-1)
- ...
12
((C1_D + C1_C - 0.5 * C1_B - C1_A)/1 * (298.15 / 1000)^1 +
(C2_D + C2_C - 0.5 * C2_B - C2_A)/2 * (298.15 / 1000)^2 + ...
(C3_D + C3_C - 0.5 * C3_B - C3_A)/3 * (298.15 / 1000)^3 +
(C4_D + C4_C - 0.5 * C4_B - C4_A)/4 * (298.15 / 1000)^4 + ...
(C5_D + C5_C - 0.5 * C5_B - C5_A)/-1 * (298.15 / 1000)^-1))
* 100;
EntR2 = ((-110.52 - 241.83) - (-115.90 + 0.5 * 0)) * 1000 + ...
(((C1_E + C1_D - 0.5 * C1_B - C1_C)/1 + (T_r / 1000)^1 +
(C2_E + C2_D - 0.5 * C2_B - C2_C)/2 * (T_r / 1000)^2 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C)/3 * (T_r / 1000)^3 +
(C4_E + C4_D - 0.5 * C4_B - C4_C)/4 * (T_r / 1000)^4 + ...
(C5_E + C5_D - 0.5 * C5_B - C5_C)/-1 * (T_r / 1000)^-1)
- ...
((C1_E + C1_D - 0.5 * C1_B - C1_C)/1 + (298.15 / 1000)^1 +
(C2_E + C2_D - 0.5 * C2_B - C2_C)/2 * (298.15 / 1000)^2 + ...
(C3_E + C3_D - 0.5 * C3_B - C3_C)/3 * (298.15 / 1000)^3 +
(C4_E + C4_D - 0.5 * C4_B - C4_C)/4 * (298.15 / 1000)^4 + ...
(C5_E + C5_D - 0.5 * C5_B - C5_C)/-1 * (298.15/ 1000)^-1))
* 1000;
end
Cp_R = 1 / (y_A / Cp_A + y_B / Cp_B + y_C / Cp_C + y_D / Cp_D + y_E /
Cp_E + y_F / Cp_F);
C1_G = (+1.74 * 10^-5); C2_G = (+2.99 * 10^-7); % [J / kg K]
Cp_G = C1_G + C2_G * (T_r);
% Primer cálculo de Cinéticas de Reacción: Para una primera aproximación, se
emplean los parámetros cinéticos indicados desde el comienzo y se calculan
% con la temperatura del flujo de reacción. Se verifican las unidades con
el fin de establecer los resultados correctos.
k1 = A_k1 * exp(-Eact_k1 / (ConstR * T_r)) / 101.32; % [mol / kgcat Pa]
K1A = A_K1A * exp(-Eact_K1A / (ConstR * T_r)) / (1000 * 101.32); % [1 /
Pa]
k2 = A_k2 * exp(-Eact_k2 / (ConstR * T_r)) / 101.32; % [mol / kgcat Pa]
K2D = A_K2D * exp(-Eact_K2D / (ConstR * T_r)) / (1000 * 101.32); % [1 /
Pa]
% Aspectos de Calor Se calculan coeficientes asociados a la conductividad
cinética:
Pr_R = (Mu_R * Cp_R) / (CondTerm_R * PM_R);
h_CalorR = (0.458 / Vacio) * (Cp_R * nR_R) * (dp * nR_R / Mu_R)^ -0.407 /
(Mu_R * Cp_R / CondTerm_R) ^ (2/3); % [W / m2 s]
JD = (0.458 / Vacio) * (dp * nR_R / Mu_R) ^ -0.407;
JH = (h_CalorR / (Cp_R * nR_R)) * (Mu_R * Cp_R / CondTerm_R) ^ (2/3);
h_CalorPelI = 2.26 * (Re_R ^ 0.8) * (Pr_R ^ 0.33) * CondTerm_R * exp(-6 *
dp / di) / di; % [W / m2 s]
h_CalorPelE = JH * (Cp_G * Mu_G / CondTerm_G) ^ 0.33 * (CondTerm_G /
do); % [W / m2 s]
U = 1 / ((1 / h_CalorPelI) + (d / CondTerm_Tubo) * (di / MLDd) + (1 /
h_CalorPelE) * (di / do) + Censu); % [W / m2 s]
13
% Resolución de expresión algebráica:
n1 = 1; n2 = 1;
syms Balance(T_s)
%
Balance(T_s) = (h_CalorR * ap * (T_s - T_r)) - ...
(n1 * EntR1 * (k1 * (PR_in * (Q_Ain / Q_A) * (FR_A /
FR_Ain) * (T_r / T_rin))) / ...
(1 + K1A * (PR_in * (Q_Ain / Q_A) * (FR_A / FR_Ain) *
(T_r / T_rin)))) - ...
(n2 * EntR2 * (k2 * (PR_in * (Q_Cin / Q_C) * (FR_C /
FR_Cin) * (T_r / T_rin))) / ...
(1 + K2D * (PR_in * (Q_Din / Q_D) * (FR_D / FR_Din) *
(T_r / T_rin))));
%{
Balance(T_s) = (h_CalorR * ap * (T_s - T_r)) - ...
(n1 * EntR1 * (k1 * FR_A) / ((Q_in * (FR / FR_in) * (T_r /
T_rin) * (PR_in / PR) * (1 / (ConstR * T_r)) + K1A * FR_A))) - ...
(n2 * EntR2 * (k2 * FR_C) / ((Q_in * (FR / FR_in) * (T_r /
T_rin) * (PR_in / PR) * (1 / (ConstR * T_r)) + K2D * FR_D)));
%}
T_s = vpasolve(Balance);
T_sup = double(T_s);
% Segundo cálculo de Cinéticas de Reacción: Para una segunda aproximación, se
emplean los parámetros cinéticos indicados desde el comienzo y se calculan,
% ahora, con la temperatura superficial. Se verifican las unidades con el
fin de establecer los resultados correctos:
k1 = A_k1 * exp(-Eact_k1 / (ConstR * T_sup)) / 101.32; % [mol / kgcat Pa]
K1A = A_K1A * exp(-Eact_K1A / (ConstR * T_sup)) / (1000 * 101.32); % [1 /
Pa]
k2 = A_k2 * exp(-Eact_k2 / (ConstR * T_sup)) / 101.32; % [mol / kgcat Pa]
K2D = A_K2D * exp(-Eact_K2D / (ConstR * T_sup)) / (1000 * 101.32); % [1 /
Pa]
% Paquete de Ecuaciones: Cada diferencial se encuentra igualado a la
expresión:
% Variables Dependientes:
% - y(1): Flujo Molar 1
% - y(2): Flijo Molar 2
% - y(3): Tr
% - y(4): Tc
% - y(5): PR
Ec = zeros(5,1);
Tubos = 1;
% Ecuación 1: Balance de Materia de Metanol (A):
%
Ec1 = -((k1 * FR_A) / ((Q_in * (FR / FR_in) * (T_r / T_rin) * (PR_in /
PR)) + K1A * FR_A) * (1 / (ConstR * T_r))) ...
* Ro_Lecho * A_transi * Tubos; % A
%
14
Ec1 = -(k1 * (PR_in * (Q_Ain / Q_A) * (FR_A / FR_Ain) * (T_r /
T_rin))) / ...
(1 + K1A * (PR_in * (Q_Ain / Q_A) * (FR_A / FR_Ain) * (T_r /
T_rin))) * Ro_Lecho * A_transi * Tubos; % B
% Ecuación 2: Balance de Materia de Formaldehído (C) o de Monóxido de
Carbono (E)
CasoEc2 = 2;
switch(CasoEc2)
case 1 % En torno a Formaldehído (C):
%
Ec2 = ((k1 * FR_A) / ((Q_in * (FR / FR_in) * (T_r / T_rin) * (PR_in / PR)
* (1 / (ConstR * T_r))) + K1A * FR_A) - ...
(k2 * FR_C) / ((Q_in * (FR / FR_in) * (T_r / T_rin) * (PR_in / PR)
* (1 / (ConstR * T_r))) + K2D * FR_D)) ...
* Ro_Lecho * A_transi * Tubos; % A
Ec2 = ((k1 * (PR_in * (Q_Ain / Q_A) * (FR_A / FR_Ain) * (T_r /
T_rin))) / ...
(1 + K1A * (PR_in * (Q_Ain / Q_A) * (FR_A / FR_Ain) * (T_r /
T_rin))) - ...
(k2 * (PR_in * (Q_Cin / Q_C) * (FR_C / FR_Cin) * (T_r /
T_rin))) / ...
(1 + K2D * (PR_in * (Q_Din / Q_D) * (FR_D / FR_Din) * (T_r /
T_rin)))) * Ro_Lecho * A_transi * Tubos; % B
%
%
case 2 % En torno a Monóxido de Carbono (E):
%
Ec2 = (k2 * (PR_in * (Q_Cin / Q_C) * (FR_C / FR_Cin) * (T_r /
T_rin))) / ...
(1 + K2D * (PR_in * (Q_Din / Q_D) * (FR_D / FR_Din) * (T_r /
T_rin))) * Ro_Lecho * A_transi * Tubos;
end
%
% Ecuación 3: Balance de Energía de Fluido de Reacción - Se encuentra
dTr / dx igualado a la expresión:
Ec3 = (-h_CalorR * ap * Ro_Lecho * (T_r - T_sup) * A_transi + U * pi() *
di * (T_c - T_r)) * Tubos / ...
(FR_A * Cp_A + FR_B * Cp_B + FR_C * Cp_C + FR_D * Cp_D + FR_E *
Cp_E + FR_F * Cp_F);
% Ecuación 4: Balance de Energía de Fluido de Control: - Se encuentra
dTc / dx igualado a la expresión:
Ec4 = -(U * pi() * do * (T_c - T_r)) * Tubos / (f_G * Cp_G);
% Ecuación 5: Balance de Presión - Se encuentra dPR / dz igualado a la
expresión
Ec5 = -f * Ro_R * Vel_Rsup^2 / dp;
15
% Conjunto de Ecuaciones:
Ec = [Ec1; Ec2; Ec3; Ec4; Ec5];
end
% Fin :)
Published with MATLAB® R2023a
16
