En esta entrada analizaremos el modelado fenomenológico de un reactor de Van de Vusse NO isotérmico, el cual será un proceso multivariable MIMO y será de gran ayuda para ejemplificar varios controladores Multivariables que trataremos en el Blog Control Automático Educación.
El modelado del Reactor de Van de Vusse No Isotérmico es básicamente el mismo modelado explicado en TOTAL DETALLE en la entrada pasada (Por favor dar Click Aquí para ver la entrada pasada si aún no la has visto) con la diferencia que vamos adicionar una tercera ecuación que representará el comportamiento de la temperatura del reactor, lo que nos permitirá ejercer control tanto en la concentración del producto B cuanto en la temperatura del reactor (Sistema MIMO). El esquema del proceso se muestra en la siguiente figura:
Ecuaciones Diferenciales del Reactor de Van de Vusse
A continuación son descritas las ecuaciones analíticas del balance de masa para los componentes A e B y el balance de energía en el reactor. Para esto es adoptada la siguiente nomenclatura:
-
concentración de A en la entrada del reactor;
-
la temperatura de la carga;
y
son las concentraciones de A y B en el reactor;
es la temperatura del reactor;
es el flujo de alimentación.
y es considerado:
- que el reactor es perfectamente agitado;
- que la densidad del líquido
y capacidad calorífica
constantes;
- Flujos iguales a
en la entrada y en la salida del reactor;
- reacción
de 1a orden en relación a A;
- reacción
de 1a orden en relación a B;
- reacción
de 2a orden en relación a A;
Los parámetros de las tasas de reacción ,i = 1,2,3, dependen de la temperatura por medio de la ley de Arrhenius
Donde ,i = 1,2,3, son las energías de activación de las tres diferentes reacciones que ocurren en el sistema y
es la constante universal de los gases.
- En el flujo de entrada solo se considera el reagente A;
- Se considera la dinámica de la camisa de refrigeración despreciable;
temperatura de la camisa del reactor constante;
- Calor transferido del reactor para la camisa como
, donde
es el coeficiente de transferencia de calor e
la área de superficie del reactor;
- Las temperaturas de la reacción son llamados como
.
El proceso es descrito por el modelo no-adiabático representado por las ecuaciones de balance de masa por componente y energía en el reactor.
Los valores de los parámetros que describen este sistema fueron extraídos de (TRIERWEILER, 1997) e están descritos en la siguiente Tabla.
Linealización del Reactor
Aplicando el concepto de la Jacobiana, explicado en detalle en la entrada anterior (Click aquí para ver la entrada pasada) linealizamos el sistema sobre un punto de operación. El punto de operación escogido se da con la entrada de y la entrada de temperatura
.
Para este punto de Operación se llega a la siguiente matriz de funciones de transferencia MIMO 2×2:
A continuación vemos una comparación entre el sistema Lineal y el Sistema NO Lineal haciendo una variación del tipo escalón igual a 5 sobre la variable de flujo, .
Código del modelo del Reactor en Matlab
A continuación se muestra el código en MATLAB que simula el modelo fenomenológico del reactor de Van de Vusse, realiza la linealización sobre el punto de operación y calcula la matriz de funciones de transferencia del sistema.
Archivo principal del programa, llamar como “mainVV.m”
%% Sergio Castaño % Rio de Janeiro - 2017 % https://controlautomaticoeducacion.com/ % Reactor de Van de Vusse NO isotermico proceso MIMO %% clc close all clear all %Entradas u.Q=128.95; %Entrada do sistema 1 (Fluxo) 128.95 u.F=72; %Entrada do sistema 2 (Temp) 72 %Condiciones Iniciales CA0d = 5.1; % mol/l CB0d = 0; % mol/l T0d = 130; % C x0 = [CA0d CB0d T0d]; %Periodo de Muestreo Ts=0.1; nm=100; t = 0:Ts:(nm-1)*Ts; %Tiempo de Simulación tsim = 0:Ts:0.1; %% Linealizacion del sistema por Jacobiana %DADOS k10 = 1.287e12; %h^-1 k20 = 1.287e12; %h^-1 k30 = 9.043e9; %L/molA.h E1 = -9758.3; %K E2 = -9758.3; %K E3 = -8560.0; %K deltaAB = -4.20; %kJ/molA deltaBC = 11.00; %kJ/molB deltaAD = 41.85; %kJ/molA p = 0.9342; %Kg/L cp = 3.01; %Kj/Kg.K Kw = 4032; %Kj/h.K.m^2 Ar = 0.215; %m^2 V = 10; %l %Valores iniciales T0 = x0(3); %ºC Ca0 = x0(1); %molA/L %Entradas Tk = u.Q; %ºC FV = u.F; %h^-1 par = [k10 k20 k30 E1 E2 E3 deltaAB deltaBC deltaAD p cp Kw Ar V Tk T0 Ca0 FV]; % Encuentra el Estado Estacionario para las entradas dadas X = fsolve(@(x)VanVusseModel(t,x,u),x0); %Establesco los Estados en el Estado Estacionario Cae = X(1) Cbe = X(2) Te = X(3) %Cálculo das constantes cinéticas no E.E. K1e = k10*exp(E1/(Te+273.15)); K2e = k20*exp(E2/(Te+273.15)); K3e = k30*exp(E3/(Te+273.15)); %% Linealización usando Jacobiana % NOTA: En el código de hoy vamos a hacer la Jacobiana Manualmente solo para % mostrar un alternativa de como hacer este tipo de linealización % manualmente. Pueden usar la misma funcion de MATLAB si lo desean tal y % como vimos en la entrada pasada del Reactor de Van de Vusse, para eso % deben colocarlo en el mismo formato del ejemplo de la entrada pasada, que % puenden encontrar en: % % https://controlautomaticoeducacion.com/modelado-de-un-reactor-cstr-de-van-de-vusse/ % % A = jacobian(fx,x); % B = jacobian(fx,u); %Cálculo da matriz A a11 = -FV-K1e-2*K3e*Cae; a12 = 0; a13 = K1e*Cae*E1./(Te+273.15)^2 + K3e*(Cae^2)*E3./(Te+273.15)^2; a21 = K1e; a22 = -FV-K2e; a23 = -K1e*Cae*E1./(Te+273.15)^2 + K2e*Cbe*E2./(Te+273.15)^2; a31 = (1/(p*cp))*(K1e*deltaAB+2*K3e*Cae*deltaAD); a32 = (1/(p*cp))*(K2e*deltaBC); a33 = (1/(p*cp))*(-K1e*Cae*deltaAB*E1./(Te+273.15)^2 - K2e*Cbe*deltaBC*E2./(Te+273.15)^2 - K3e*(Cae^2)*deltaAD*E3./(Te+273.15)^2)-FV-(Kw*Ar)/(p*cp*V); A = [a11 a12 a13 a21 a22 a23 a31 a32 a33]; %Cálculo da matriz B B = [Ca0-Cae 0 -Cbe 0 T0-Te (Kw*Ar)/(p*cp*V)]; %Cálculo da matriz C C = [0 0 0 0 1 0 0 0 1]; %Cálculo da matriz D D = [0 0 0 0 0 0]; % %Determino la Funcion de Transferencia para la entrada 1 Cb % display('Representación en Función de Transferencia para Cb') % [num1,den1]=ss2tf(A,B,C,D,1); % ft1=tf(num1,den1) H = ss(A,B,C,D); [z,p,k] = zpkdata(H); K = zpk(z,p,k) %Comprobar Modelo Lineal VS Modelo NO Lineal du=5; u.F=u.F+du; [t,x] = ode45(@(t,x)VanVusseModel(t,x,u), tsim , X, u); tsm=0:0.001:0.1; in=zeros(length(tsm),2); in(:,1)=du*ones(length(tsm),1); in(:,2)=0*ones(length(tsm),1); [ylin]=lsim(H,in,tsm); %Grafica figure plot(t,x(:,2),'-b',tsm,ylin(:,2)+X(2),'--r','LineWidth',2),ylabel('C_b'); title('Comparación modelo Lineal vs No Lineal (C_b)'); xlabel('tiempo (min)'); ylabel('C_b mol/litro'); legend('No Lineal','Lineal','Location','SouthEast') figure plot(t,x(:,3),'-b',tsm,ylin(:,3)+X(3),'--r','LineWidth',2),ylabel('Q'); title('Comparación modelo Lineal vs No Lineal (T)'); xlabel('tiempo (min)'); ylabel('T (K)'); legend('No Lineal','Lineal','Location','SouthEast')
Archivo función del modelo del Reactor, nombrar el archivo como “VanVusseModel.m”
function dydt = VanVusseModel(~,x,u) %% Función que representa la dinámica del reactor de Van de Vusse % x = Estados del sistema % u = Estructura con las dos entradas del sistema %% Datos % Entradas Tk = u.Q; %ºC (Temperatura de la camisa) fov = u.F; %h^-1 (Flujo) %Constantes de equlibrio k10 = 1.287e12; %h^-1 k20 = 1.287e12; %h^-1 k30 = 9.043e9; %L/molA.h E1 = -9758.3; %K E2 = -9758.3; %K E3 = -8560.0; %K deltaAB = -4.20; %kJ/molA deltaBC = 11.00; %kJ/molB deltaAD = 41.85; %kJ/molA p = 0.9342; %Kg/L cp = 3.01; %Kj/Kg.K Kw = 4032; %Kj/h.K.m^2 Ar = 0.215; %m^2 V = 10; %l T0 = 130; %ºC (Temperatura a la entrada) Ca0 = 5.1; %molA/L (Concentracion a la entrada) % % Notación para las Variables de Estado % ca = x(1); %Concentración A cb = x(2); %Concentración B T = x(3); %Temperatura %Constantes cinéticas k1 = k10*exp(E1./(T+273.15)); k2 = k20*exp(E2./(T+273.15)); k3 = k30*exp(E3./(T+273.15)); % % Ecuaciones Diferenciales del Sistema % dydt = zeros(3,1); dydt(1) = fov*(Ca0 - ca) -k1*ca - k3*ca*ca; dydt(2) = -fov*cb + k1*ca - k2*cb; dydt(3) = (1/(p*cp))*(k1*ca*deltaAB + k2*cb*deltaBC + k3*ca^2*deltaAD)... + fov*(T0-T) + (Kw*Ar/(p*cp*V))*(Tk-T);
Bibliografia
- TRIERWEILER, J. (1997). A Systematic Approach to Control Structure Design, Tesis de Doctorado. Alemania: Universität Dortmund.
Mi nombre es Sergio Andres Castaño Giraldo, y en este sitio web voy a compartir una de las cosas que mas me gusta en la vida y es sobre la Ingeniería de Control y Automatización. El sitio web estará en constante crecimiento, voy a ir publicando material sobre el asunto desde temas básicos hasta temas un poco más complejos. Suscríbete al sitio web, dale me gusta a la página en Facebook y únete al canal de youtube. Espero de corazón que la información que comparto en este sitio, te pueda ser de utilidad. Y nuevamente te doy las gracias y la bienvenida a control automático educación.