Saltar al contenido
Control Automático Educación

4. Linealización de un Reactor de Van de Vusse NO Isotérmico (Proceso MIMO)

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:

Reactor de Van de Vusse No Isotermico

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:

  •  C_{a_0} concentración de A en la entrada del reactor;
  •  T_0 la temperatura de la carga;
  • C_a y C_b son las concentraciones de A y B en el reactor;
  • T es la temperatura del reactor;
  • F es el flujo de alimentación.

y es considerado:

  • que el reactor es perfectamente agitado;
  • que la densidad del líquido \rho y capacidad calorífica C_p constantes;
  • Flujos iguales a F en la entrada y en la salida del reactor;
  • reacción A\rightarrow B de 1a orden en relación a A;
  • reacción B\rightarrow C de 1a orden en relación a B;
  • reacción 2A\rightarrow D de 2a orden en relación a A;

Los parámetros de las tasas de reacción k_i,i = 1,2,3, dependen de la temperatura por medio de la ley de Arrhenius

K_i(T)=k_{i0} exp\left( \dfrac{-E_i/R}{T(C)+273.15} \right)

Donde E_i,i = 1,2,3, son las energías de activación de las tres diferentes reacciones que ocurren en el sistema y R 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;
  • T_k temperatura de la camisa del reactor constante;
  • Calor transferido del reactor para la camisa como Q = K_w A_R (T - T_k), donde K_w es el coeficiente de transferencia de calor e A_R la área de superficie del reactor;
  • Las temperaturas de la reacción son llamados como (-\Delta H_{RAB}),(-\Delta H_{RBC}),(-\Delta H_{RAD}).

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.

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1(T)C_A-K_3(T)C_A^2

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1(T)C_A-K_2(T)C_B

\dfrac{d(T)}{dt}=\dfrac{1}{\rho C_p} \left[k_1(T)C_a(-\Delta H_{RAB})+k_2(T)C_b(-\Delta H_{RBC})+k_3(T)C_a^2(-\Delta H_{RAD})\right]+\dfrac{F}{V} (T_0-T)+\dfrac{K_w A_R}{\rho C_p V} (T_k-T)

Los valores de los parámetros que describen este sistema fueron extraídos de (TRIERWEILER, 1997) e están descritos en la siguiente Tabla.

parametros reactor cstr

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 \dfrac{F}{V}=72h^{-1} y la entrada de temperatura T_k=128.95K.

Para este punto de Operación se llega a la siguiente matriz de funciones de transferencia MIMO 2×2:

P_n=\begin{bmatrix} \dfrac{-1.1033 (s+26.52) (s+94.31)}{(s+162.2) (s^2 + 203s + 1.284e04)} & \dfrac{142.23 (s+31.3)}{(s+162.2) (s^2 + 203s + 1.284e04)}\\ \\ \dfrac{-6.7606 (s+159.2) (s-7.671)}{(s+162.2) (s^2 + 203s + 1.284e04)} & \dfrac{30.829 (s+131) (s+168.8)}{(s+162.2) (s^2 + 203s + 1.284e04)} \end{bmatrix}

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, \dfrac{F}{V}=77h^{-1}.

Reactor CSTR Lineal vs no lineal
Reactor CSTR Lineal vs no lineal

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.

Entradas relacionadas

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.

Comentarios (2)

Hay problemas al correr el programa principal
Index exceeds matrix dimensions.

Error in VanVusseModel (line 9)
Tk = u(2); %ºC (Temperatura de la camisa)

Error in @(x)VanVusseModel(t,x,u)

Error in fsolve (line 230)
fuser = feval(funfcn{3},x,varargin{:});

Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.
Igual sucede en la funcion:
function dydt = VanVusseModel(~,x,u)

Error: Function definitions are not permitted in this context.
Teng instalado MATLAB2017a

Responder

Lo rodé aquí y funciona sin problemas tanto en Matlab 2015 como en el 2020. Debes guardar los dos archivos en la misma carpeta, y solo ejecutar el programa principal. El programa principal llamalo como mainVV.m y la función del modelo como VanVusseModel.

Responder