Linealización de un Reactor de Van de Vusse NO Isotérmico

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:

Contenidos mostrar

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.

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}$ .

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.

Sergio A. Castaño Giraldo

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.