Curso básico de modelado de reactores químicos en MATLAB-Octave UPM

Aquí recogemos el material del curso gratuito que tuvo lugar durante la primera y segunda semana de diciembre de 2013 para alumnos de grado de Ingeniería Química. Octave UPM que es la alternativa a MATLAB más interesante para los que busquen máxima compatibilidad.

[box]Accede a los contenidos actualizados de la segunda edición de este curso

Horario y aulas

Horario de 16 a 18 h para los siguientes días de diciembre, y en las aulas de informática de la EPS:

martes 3 (L25);
miercoles 4 (L18);
jueves 5 (L13);

martes 10 (L25);
miercoles 11 (L25);
jueves 12 (L13);

¡Aquellos alumnos que asistan a un 80% del curso, recibirán un certificado de asistencia!
¿Por qué Octave UPM y no MATLAB?

Octave UPM es una adaptación de GNU Octave que se emplea para la docencia de Informática en la Escuela de Ingenieros de Caminos, Canales y Puertos de la UPM y ha sido desarrollado por Israel Herraiz (@herraiz)
Es altamente compatible con MATLAB.
Es software libre y gratuito, disponible para Windows y Linux

Entorno de programación de Octave UPM (similar a MATLAB)

Descargar GNU Octave (Windows), ¡es gratis!

Descargar Octave UPM (versión antígua)

Actualización: La nueva versión de GNU Octave ya incorpora interfaz gráfica (descargar para Windows)

El curso es gratuito y de acceso libre, todo el material utilizado durante el mismo será también públicado para aquellos que no puedan asistir. No obstante, recomendamos este curso ya que su intención es más práctica que teórica y además estaremos a vuestra disposición para ayudaros a encontrar fallos en vuestros scripts así como resolver dudas puntuales que tengáis. Recomendamos que vengáis con vuestro portátil para poder trabajar siempre que queráis con Octave UPM en vuestro ordenador.

¿Por dónde empezar?

Los requisitos del curso serán tener conocimientos básicos de programación y de reactores químicos. Para el que quiera ir repasando:

MateWiki (UPM): Video tutoriales y materiales para aprender a programar en MATLAB-Octave
Introducción a MATLAB y Octave (PDF)

Encuesta de calidad

Día 1:

Iniciacion al modelado de reactores quimicos com MATLAB – Octave (intro) por CAChemE

Día 2 y 3:

Simulación de reactores químicos con octave por CAChemE

Días 4 y 5:

Reactor flujo piston en MATLAB – Octave UPM – Craqueo térmico por CAChemE

Día 6:

Reactor de flujo piston con MATLAB Octave from CAChemE

Archivo principal:

pkg load odepkg

close all
clear all
clc

global k0 Ea R Seccion Qv1 alfa U Tw Tref d DCp DHref Cp Qv2

%Sacarosa(A) medioacido(B) glocosa(C) fructosa(D) indes(E)

%DATOS TERMODINAMICOS Y CINETICOS
%        A  B C D E 
alfa = [-1 -1 1 1 0
        -1 -1 0 0 2];

Cp = [62.0 75.2 80.0 80.0 80.0]; % kJ/kmol.K
PM = [342 18 180 180 180]; % Kg/Kmol

DHref = [-11000 -12000]; % kJ/kmol
Ea = 125300; % J/mol
R = 8.314; % kJ/kmol.K o J/mol.K
k0 = 0.0193/(exp(-Ea/R/(273+50)));
DCp = Cp*alfa'; % Equivalente a sum(Cp.*alfa)
Tref = 25+273; % K

% DATOS DE LOS REACTORES Y REFRIGERACION
Seccion = 0.1; % m2
d = (Seccion*4/pi)^0.5; % m
U = 1.0; % kJ/s·m2·K
Tw = 30+273; % K

% DATOS ALIMENTO
Talim = 50+273; % K
Qv0 = 0.4/60; % m3/s
Qv1 = 0.2/60; % m3/s
Qv2 = 0.3/60; % m3/s
Qv3 = 0.1/60; % m3/s

% EJECUTO PRIMER REACTOR Y TANTEANDO VOY A VER QUE LONGITUD DE REACTOR
% NECESITO PARA LLEGAR A LA SALIDA A C(3)=1.5M

C0 = [2 16.6 0.1 0 0]; % kmol/m3
L1 = 1; % no conocemos este valor, se pone 1 para empezar

y0 = [C0 Talim]; % Condiciones iniciales
[x1,Y1] = ode45(@fdif_rfp1,[0 L1],y0);

% REPRESENTACIÓN DE RESULTADOS
L_r1 = x1; % variable independiente
C_r1 = Y1(:,1:5);
T_r1 = Y1(:,6);
C_salida_r1 = Y1(end,1:5); % variable dependiente (C)
T_salida_r1 = Y1(end,6); % variable dependiente (Tª)

figure(1)
plot(L_r1,C_r1)

xlabel('L / m')
ylabel('C_j / kmol/m^3')
title('Variación de concentraciones en el RFP_1')
legend('Sacarosa', 'Agua acidulada', 'glucosa', 'Fructosa', 'indesada')

figure(2)
plot(L_r1,T_r1)
xlabel('L / m')
ylabel('T / K')
title('Variación de temperatura en el RFP_1')

Archivo de función:

function z=fdif_rfp1(x,y)

global k0 Ea R Seccion Qv1 alfa U Tw Tref d DCp DHref Cp

% Se asginan de nuevo las variables del vector de
% condiciones iniciales "y"

%Sacarosa(A) medioacido(B) glocosa(C) fructosa(D) indes(E)

C=y(1:5);
T=y(6);

n = C*Qv1; % Flujo molar;

% BALANCE DE MATERIA

k1 = k0*exp(-Ea/R/T);
k2 = exp(-5.64-220/T);

k = [k1 k2];

r(1) = k(1)*C(1);
r(2) = k(2)*C(1);

dndL(1) = Seccion*r*alfa(:,1);
dndL(2) = Seccion*r*alfa(:,2);
dndL(3) = Seccion*r*alfa(:,3);
dndL(4) = Seccion*r*alfa(:,4);
dndL(5) = Seccion*r*alfa(:,5);

dCdL = dndL/Qv1;

% % Equivalente pero mediante bucle:
%  for i=1:5;
%    z(i)=(Seccion/Qv1)*r*alfa(:,i);
%  end

% BALANCE DE ENERGIA 
for i=1:2;
    DH(i)=DHref(i)+DCp(i)*(T-Tref);
end

Q = U*pi*d*(Tw-T);

dTdL = (Q - Seccion*sum(r.*DH))/ sum(n'.*Cp);

z = [dCdL, dTdL]';

end

Por CAChemE