Tanque de Nivel - Modelo Matemático, Simulación y Control

En esta entrada aprenderemos a modelar, linealizar y simular el modelo matemático de un tanque de nivel no lineal, al rededor de un punto de operación.

Si lo deseas puedes ver nuestro curso de análisis de sistema para aprender a modelar otro tipo de sistemas y suscribirte al canal de youtube para seguir aprendiendo más sobre el modelado de procesos.

Si tienes interés en modelar otros sistema hidráulico puedes darle un vistazo a las otras entradas del sitio web:

Modelo matemático de Tanques en Paralelo y simulado en Simulink
Modelo matemático de Tanques Esféricos y Simulados en Matlab
Modelado matemático de Dos Tanques Interconectados en Serie y simulación en matlab y simulink.

Contenidos mostrar

Sistema de un Tanque de Nivel

Consideremos el siguiente sistema de control conformado por un tanque de nivel con accionamiento por válvulas donde el objetivo es mantener constante el nivel a pesar de las variaciones de la demanda

En esta entrada vamos a aprender a realizar el modelado de un tanque y al final de este post veremos un código de implementación adicional para la implementación de un sistema de control de nivel de líquido de un tanque a través de Matlab.

Con base al modelo, podremos realizar diferentes tareas, como por ejemplo diseñar algún tipo de controlador que mantenga las especificaciones de nivel que necesitemos.

Recordemos la importancia de entender el funcionamiento de procesos, que en la teoría de control, son fundamentales para poder diseñar reguladores que mantengan las variables de proceso donde deseamos que ellas esten.

A nivel practico, entender el funcionamiento de un proceso, nos va a tomar el 90% de nuestro tiempo y esfuerzo, ya que una vez obtenido el modelo y una vez, que entendemos como funciona el proceso industrial, realizar el calculo de controladores, se convierte en una tarea mucho mas sencilla.

El proceso hidraulico está conformado por un tanque de nivel de agua junto con dos válvulas. Una será la válvula de control en cuanto la otra será una válvula manual usada como perturbación.

Video en Español

Nota: Las unidades de $k_2$ en el video en el minuto 06:35 son realmente $m^2$ porque en el balance de energias está actuando la energía de presión de altura en este punto.

Video em Português

Modelo Estático del Nivel de un Tanque

Para determinar la velocidade de flujo que sale por una boquilla lisa y redonda en la parte inferior del tanque de nivel como es mostrado en la figura, se aplica el principio de Bernoulli entre el punto de referencia 1 en la superficie y el punto de referencia 2 en la boquilla del tubo.

Si se considera que la presión manométrica entre los puntos 1 y 2 son cero. Que m es la masa de las partículas de fluido entre el punto 1 y 2. Y que v es la velocidad de salida de la particula 1 y 2. se tiene que:

$\dfrac{1}{2}mv^2=mgH$

$v=\sqrt{2gH}$

El flujo de salida del tanque para el caso de una sección transversal es:

$Q=Sv=S\sqrt{2gH}$

Ecuaciones de las Válvulas

En la entrada del elemento final de control explicamos en detalle el modelado matemático de una válvula y como puede expresarse su comportamiento dinámico en relación al tipo de válvula instalada en el proceso.

Considerando el flujo que circula a traves de la válvula se encuentra en estado estacionário, podemos plantear la siguiente ecuación:

$Q_v=K_vA\sqrt{\Delta P}$

$Q_v$ : Flujo a través da válvula

$K_v$ : una constante

$A$ : Área de paso

$\Delta P$ : Presión diferencial a través de la válvula. P2-P1

Se puede concluir que el flujo que pasa por la válvula es proporcional al área de abertura de la válvula en el caso que la diferencia de presión sea constante.

De manera practica tomamos una válvula con un comportamiento inteligente, donde sea posible hacer una aproximación más o menos lineal entre el flujo Qv y la abertura de la válvula.

$Q_s=k_2a_2 \sqrt{2gH}$

Vamos a suponer que el flujo de entrada Qe es proporcional a la abertura de la válvula de entrada considerando un suministro constante.

$Q_e=k_1a_1$

Sistemas y Modelos

Error en Estado Estacionario

Modelado de un Reactor CSTR de Van de Vusse

Modelo del Tanque

Para ir entendiendo el modelo del tanque, vamos a asignar algunos valores numéricos al proceso.

Este tranque puede trabajar de dos formas diferentes, ya que posee dos entradas que influencian el nivel del mismo. Las dos entradas son las dos válvulas del proceso y este puede operar así:

Tanque controlado por la válvula de entrada, considerando $a_2$ como perturbación
Tanque controlado por la válvula de salida, considerando $a_1$ como perturbación

Consideremos los siguientes parámetros para nuestro modelo:

$A=0.5\ {\rm m^2}$ (tanque cilíndrico)
$H_{max}=1\ {\rm m}$
$k_1=0.05\ {\rm \dfrac{m^3}{s}}$
$k_2=0.015\ {\rm m^2}$
$a_1=0.6$ (abertura de equilíbrio válvula de entrada)
$a_2=0.5$ (abertura de equilíbrio válvula de salida)
$g=10\ {\rm \dfrac{m}{s^2}}$ gravedad usada para los cálculos

Modelado Matemático de un Tanque de Nivel

Para obtener las ecuaciones matemáticas del tanque de nivel partimos inicialmente del balance de masas del sistema para obtener las ecuaciones diferenciales (EDO) del tanque:

En este balance, básicamente todo lo que le entra al tanque menos todo lo que le sale al tanque es igual a su capacitancia:

A\dfrac{dH}{dt}=q_e-q_s

Reemplazando en $Q_e$ y $Q_s$ :

A\dfrac{dH}{dt}=k_1a_1-k_2 a_2 \sqrt{2gH}

De la ecuacion del modelo del tanque de nivel anterior se nota que el sistema hidráulico corresponde a un sistema NO LINEAL por causa de que la variable de interés H (Altura) se encuentra dentro de una raíz cuadrada.

Da click en el siguiente enlace si deseas saber que es un sistema lineal y no lineal.

En este caso como queremos representar el sistema hidráulico del tinaco de nivel por funciones de transferencia, vamos a tener que linealizar las edo del tanque utilizando las Series de Taylor.

Linealización del Tanque de Nivel

Procedemos a encontrar el punto de equilibrio de la ecuación diferencial del tanque de nivel volviendo para eso todas las derivadas iguales a cero.

$0=k_1a_1^*-k_2a_2^* \sqrt{2gH_o}$

$k_1a_1^*=k_2a_2^* \sqrt{2gH_o}$

Eso quiere decir que todo lo que le entra al tanque será igual a todo lo que le sale al tanque. Note que en este caso, hemos cambiado la variable $H$ por $H_o$ y $a_1$ , $a_2$ por $a_1^*$ , $a_2^*$ para identificar que son variables de equilibrio del tanque.

Despejando $H_o$ vamos a saber cual es la altura en estado estacionario o de equilibrio:

$H_o=\left(\dfrac{k_1a_1^*}{k_2a_2^*\sqrt{2g}} \right)^2$

Reemplazando los valores de los parámetros tenemos:

$H_o = 0.8$ m

Linealización con relación a la válvula de entrada

Teniendo el punto de equilibrio. Linealizamos el tanque de nivel por series de Taylor considerando la abertura de la válvula de salida como constante.

La función de Taylor va a depender entonces solo de $a_1$ y $h$ igual a la capacitancia del tanque:

f(a_1,H)=A\dfrac{dH}{dt}=k_1a_1-k_2a_2 \sqrt{2gH}

La ecuación de Taylor es una serie infinita la cual expandimos únicamente hasta la primera derivada, por lo tanto viene dado por la siguiente expresión:

f(a_1,H)\approx f(a_1^*,H_o)+\left.\dfrac{\partial f}{\partial a_1}\right|_{a_1^*,H_o}(a_1-a_1^*)+\left.\dfrac{\partial f}{\partial H}\right|_{a_1^*,H_o}(H-H_o)

note que $\Delta a_1= a_1-a_1^*$ y $\Delta H= H-H_o$ que permite expresar el problema en variable de desviación. Es decir estamos trasladando el problema al punto donde queremos linealizar ( $a_1^*,H_o*$ ).

Sustituyendo la ecuación en sus puntos de equilibrio (primer término de la ecuación anterior) y aplicando la derivada parcial de la función (dos últimos términos de la expresión anterior) tenemos:

f(a_1,H)\approx k_1a_1^*-k_2a_2\sqrt{2gH_o} + k_1\Delta a_1-\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\Delta H

Sabemos que partiendo del modelo matemático del tanque:

\underset{A\dfrac{dH}{dt}}{\underbrace{f(a_1,H)}}\approx \underset{\left.A\dfrac{dH}{dt}\right|_{a_1^*,H_o}}{\underbrace{k_1a_1^*-k_2a_2\sqrt{2gH_o}}} + k_1\Delta a_1-\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\Delta H

por lo tanto despejando:

A\dfrac{dH}{dt}-\left.A\dfrac{dH}{dt}\right|_{a_1^*,H_o}\approx k_1\Delta a_1-\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\Delta H

A\dfrac{d\Delta H}{dt}\approx k_1\Delta a_1-\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\Delta H

Notemos que la ecuación anterior ahora es una ecuación lineal que depende de la variable desvío $\Delta H$ y que ésta no se encuentra dentro de ninguna raíz cuadrada.

Podemos Aplicar transformada de Laplace para llevar la ecuación lineal al plano transformado S. Para simplificar la notación, renombramos:

\Delta H = \bar{h}

\Delta a_1 = \bar{a_1}

Transformada de Laplace:

As\bar{h}(s)= k_1\bar{a_1}(s)-\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\bar{h}(s)

\left(As+\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}\right)\bar{h}(s)= k_1\bar{a_1}(s)

La función de transferencia del Tanque de Nivel manteniendo la válvula de salida constante viene dado por:

\dfrac{\bar{h}}{\bar{a_1}} = \dfrac{k_1}{As+\dfrac{k_2a_2\sqrt{2g}}{2\sqrt{H_o}}}

\dfrac{\bar{h}}{\bar{a_1}} = \dfrac{k_1/A}{s+\dfrac{k_2a_2\sqrt{2g}}{2A\sqrt{H_o}}}

sustituyendo los parámetros del modelo:

\dfrac{\bar{h}}{\bar{a_1}} = \dfrac{0.1}{s+0.0375}

Linealización con relación a la válvula de salida

Aqui se aplica exactamente el mismo procedimiento, solo que en este caso se considera la válvula de entrada como constante.

f(a_2,H)\approx f(a_2^*,H_o)+\left.\dfrac{\partial f}{\partial a_2}\right|_{a_2^*,H_o}(a_2-a_2^*)+\left.\dfrac{\partial f}{\partial H}\right|_{a_1^*,H_o}(H-H_o)

f(a_2,H)\approx \left.A\dfrac{dH}{dt}\right|_{a_2^*,H_o} - k_2\sqrt{2gH_o}\Delta a_2-\dfrac{k_2a_2^*\sqrt{2g}}{2\sqrt{H_o}}\Delta H

A\dfrac{d\Delta H}{dt}\approx - k_2\sqrt{2gH_o}\Delta a_2-\dfrac{k_2a_2^*\sqrt{2g}}{2\sqrt{H_o}}\Delta H

renombrando:

\Delta H = \bar{h}

\Delta a_2 = \bar{a_2}

Y aplicando transformada de Laplace, llegamos a la función de transferencia del tanque de nivel con relación a la válvula de salida:

\dfrac{\bar{h}}{\bar{a_2}} = \dfrac{-\dfrac{k_2}{A}\sqrt{2gH_o}}{s+\dfrac{k_2a_2^*\sqrt{2g}}{2A\sqrt{H_o}}}

Reemplazando los valores de los parámetros:

\dfrac{\bar{h}}{\bar{a_2}} = \dfrac{-0.12}{s+0.0375}

Función de Transferencia del Tanque de Nivel

Finalmente si sumamos las dos funciones de transferencia (dado que son un sistema lineal y obedecen al principio de la superposición) vamos a obtener la ecuación de la altura que representa la dinámica del tanque:

\bar{h}(s) = \dfrac{0.1}{s+0.0375}\bar{a_1}-\dfrac{0.12}{s+0.0375}\bar{a_2}

Código de Simulación del Modelo y Control del tanque de Nivel

A continuación de dejo el código de implementación del tanque en MATLAB para que lo puedas copiar y pegar en tu computador.

También te dejo un segundo ejemplo donde se aplica un control PI sobre este tanque.

Antes que nada, si te interesa el modelado y control de sistemas, te invito a realizar mi Curso de Simulink desde Cero con más de 10 horas aprendiendo a utilizar esta herramienta para la simulación de procesos y posteriormente para el control de los mismos (Certificado de Aprobación al finalizar el curso):

Comenzar Curso de Simulink con Certificado

Código del Modelado del Tanque de Nivel en Simulink/Matlab

Para correr este código, es necesario que tengas el archivo de Simulink en la misma carpeta donde grabaste el código.

Primero debes descargar el modelo del tanque en Simulink dando click en el siguiente botón:

Descargar modelo del Tanque en Simulink

Luego debes copiar el siguiente código y grabarlo en la misma carpeta donde tienes el archivo de simulink descargado en el paso anterior:

% Sistemas Realimentados
% Por: Sergio Castaño
%% Dinamica de um tanque de Nível
clc
clear all
close all
%Variaveis do tanque
k1=0.05;k2=0.015;g=10;
A=0.5;
a1=0.0:0.1:1;   %Avertura da valvula de entrada
a2=0:0.1:1;     %Avertura da valvula de Saida
k=1;

A1o=0.6; %Avertura da Valvula no Pto de Equilibrio
A2o=0.5; %Avertura da valvula de saida no equilibrio

%Acho o Pto de Equilibrio do Tanque
Ho=(k1^2*A1o^2)/(k2^2*A2o^2*2*g);

%Grafica nao linear do comportamento do Nível manipulando a valvula de
%entrada e deixando a valvula de saída constante e viceverza
for i=1:11   
       hA1(i)=(k1^2*a1(i)^2)/(k2^2*A2o^2*2*g);
       hA2(i)=(k1^2*A1o^2)/(k2^2*a2(i)^2*2*g);
       k=k+1;
end
figure
t=0:0.1:1;
subplot(2,1,1);
plot(t,hA1),grid
axis([0 1 0 1]);
title('Variação em valvula de entrada a1');
ylabel('Altura (h)');
xlabel('Abertura da Valvula a1');
subplot(2,1,2);
plot(t,hA2),grid
axis([0 1 0 1]);
title('Variação em valvula de entrada a2');
ylabel('Altura (h)');
xlabel('Abertura da Valvula a2');

%% Processo de Simulação
tsim=300;       %tempo de simulação
tA1=10;         %Tempo em que a variação de a1 acontece
tA2=150;        %Tempo em que a variação de a2 acontece

dA1=0.02;       %Variação em a1
dA2=0.03;       %Variação em a2
%Simulação
sim('tanque');
%Grafica dos resultados
figure
subplot(3,1,1);
plot(t,rY(:,1),t,rY(:,3)),grid
title('Nível do Tanque');
ylabel('Altura (h)');
xlabel('Tempo (s)');
legend('Não Linear','Linear');

subplot(3,1,2);
plot(t,rY(:,2)),grid
title('Abertura em a1');
ylabel('% a1');
xlabel('Tempo (s)');

subplot(3,1,3);
plot(t,rY(:,4)),grid
title('Abertura em a2');
ylabel('% a2');
xlabel('Tempo (s)');
clc

Código del sistema de control PI de nivel de líquido de un tanque

A continuación te dejo los archivos en Matlab y Simulink, donde al modelo del tanque anterior, se le aplica un CONTROL Proporcional y un CONTROL PI.

Dicho control PI esta sintonizado de dos formas. La primera, es por una sintonizacion por cancelamiento de polos, donde el termino $\tau_i$ se hace igual al polo del sistema linealizado del tanque.

La segunda sintonia, se hace con un Control PI por Asignación de Polos, donde ubico dos polos reales e iguales.

Descargar Archivos en Matlab y Simulink

% Sistemas Realimentados
% Por: Sergio Castaño
% Prof. Julio Elias Normey Rico
%% Dinamica de um tanque de Nível
clc
clear all
close all
%Variaveis do tanque
k1=0.05;k2=0.015;g=10;
A=0.5;
a1=0.0:0.1:1;   %Avertura da valvula de entrada
a2=0:0.1:1;     %Avertura da valvula de Saida
k=1;

A1o=0.6; %Avertura da Valvula no Pto de Equilibrio
A2o=0.5; %Avertura da valvula de saida no equilibrio

%Encuentro o Pto de Equilibrio do Tanque
Ho=(k1^2*A1o^2)/(k2^2*A2o^2*2*g);

%Grafico nao linear do comportamento do Nível manipulando a valvula de
%entrada e deixando a valvula de saída constante e viceverza
for i=1:11   
       hA1(i)=(k1^2*a1(i)^2)/(k2^2*A2o^2*2*g);
       hA2(i)=(k1^2*A1o^2)/(k2^2*a2(i)^2*2*g);
       k=k+1;
end
figure
t=0:0.1:1;
subplot(2,1,1);
plot(t,hA1),grid
axis([0 1 0 1]);
title('Variação em valvula de entrada a1');
ylabel('Altura (h)');
xlabel('Abertura da Valvula a1');
subplot(2,1,2);
plot(t,hA2),grid
axis([0 1 0 1]);
title('Variação em valvula de entrada a2');
ylabel('Altura (h)');
xlabel('Abertura da Valvula a2');

%% Processo de Simulação
tsim=400;       %tempo de simulação
tA1=100;         %Tempo em que a variação de a1 acontece
tA2=250;        %Tempo em que a variação de a2 acontece

dA1=0.5;       %Variação em a1
dA2=0.05;       %Variação em a2
%Simulação
sim('tanque1');
%Grafica dos resultados
figure
subplot(3,1,1);
plot(t,rY(:,1),t,rY(:,3)),grid
title('Nível do Tanque');
ylabel('Altura (h)');
xlabel('Tempo (s)');
legend('Não Linear','Linear');

subplot(3,1,2);
plot(t,rY(:,2)),grid
title('Abertura em a1');
ylabel('% a1');
xlabel('Tempo (s)');

subplot(3,1,3);
plot(t,rY(:,4)),grid
title('Abertura em a2');
ylabel('% a2');
xlabel('Tempo (s)');
clc

%% Projeto dos Controladores
%% Processo de Simulação

%Controle Proporcional
tsim=120;       %tempo de simulação

ref= [1 0.6 0.2 0.9 0.9]'; %Referencias

per=[0.5 0.85]'; %Perturbacoes
tper=90;

Kc=100;
%Simulação
sim('tanque_Kc');
%Grafica dos resultados
figure
subplot(2,1,1);
plot(t,rY(:,1),t,rY(:,2)),grid
title('Controle Kc=100');
ylabel('Altura (h)');
xlabel('Tempo (s)');
subplot(2,1,2);
plot(t,rY(:,3)),grid
ylabel('Abertura A2');
xlabel('Tempo (s)');

%Controle PI Cancelamento
%Por cancelamento se faz Ti igual ao polo do modelo linear

Tau=A/((k2*A2o*sqrt(2*g))/(2*sqrt(Ho)));
Kp=k1/((k2*A2o*sqrt(2*g))/(2*sqrt(Ho)));

%Por cancelamento só devo dar o tempo desejado da MF (Tr)
Tr=Tau/5;

Ti=Tau;
Kc=Ti/(Kp*Tr);

tsim=500;       %tempo de simulação

ref= [0.8 0.6 0.7 0.65 0.65]'; %Referencias

per=[0.5 0.80]'; %Perturbacoes
tper=350;

%Simulação
sim('tanque_PI_Can');
%Grafica dos resultados
figure
subplot(3,1,1);
plot(t,rY(:,1),t,rY(:,2)),grid
title('Controle PI Cancelamento');
ylabel('Altura (h)');
xlabel('Tempo (s)');
subplot(3,1,2);
plot(t,rY(:,4)),grid
ylabel('Controle');
xlabel('Tempo (s)')
subplot(3,1,3);
plot(t,rY(:,3)),grid
ylabel('Abertura A2');
xlabel('Tempo (s)')

%Por cancelamento se evidencia que a rejeicao de perturbacao é muito lenta

C=tf(Kc*[Ti 1],[Ti 0]);
P=tf(0.1,[1 0.0375]);
Q=tf(-0.12,[1 0.0375]);

YR=minreal((C*P)/(1+C*P));
YQ=minreal((Q)/(1+C*P));
figure
subplot(2,2,1)
pzmap(YR);
title('DPZ_R');
subplot(2,2,2)
rlocus(YR);
subplot(2,2,3)
pzmap(YQ);
title('DPZ_Q');
subplot(2,2,4)
rlocus(YQ);

%Controle por Alocacao de polos
%% Calculo do controlador (Alocação de Polos)

TsMA=(Tau*3)/5;       %Tempo de acomodação MA (Sistema de Primera Ordem)
ep=1;                             %Coeficente de amortecimento
Wn=3/(ep*TsMA);                     %Frequência Natural
Sd=[-ep*Wn+1i*Wn*sqrt(1-ep^2), -ep*Wn-1i*Wn*sqrt(1-ep^2)]; %Alocação de Polos
Pds=poly(Sd);

Kc=(Pds(2)*Tau-1)/Kp;        %Calculo de Kc
Ti=(Kp*Kc)/(Pds(3)*Tau);     %Calculo de ti

%Simulação
sim('tanque_PI_Can');
%Grafica dos resultados
figure
subplot(3,1,1);
plot(t,rY(:,1),t,rY(:,2)),grid
title('Controle PI Polos Reais Iguais');
ylabel('Altura (h)');
xlabel('Tempo (s)');
subplot(3,1,2);
plot(t,rY(:,4)),grid
ylabel('Controle');
xlabel('Tempo (s)')
subplot(3,1,3);
plot(t,rY(:,3)),grid
ylabel('Abertura A2');
xlabel('Tempo (s)')

C=tf(Kc*[Ti 1],[Ti 0]);
P=tf(0.1,[1 0.0375]);
Q=tf(-0.12,[1 0.0375]);

YR=minreal((C*P)/(1+C*P));
YQ=minreal((Q)/(1+C*P));
figure
subplot(2,2,1)
pzmap(YR);
title('DPZ_R');
subplot(2,2,2)
rlocus(YR);
subplot(2,2,3)
pzmap(YQ);
title('DPZ_Q');
subplot(2,2,4)
rlocus(YQ);

Simulación y Código del tanque en LABVIEW

A continuación se muestra la simulación dinámica del llenado y vaciado del tanque representado por una ecuación diferencial y discretizada por el método de EULER.

Para descargar los VI puedes hacerlo a través del siguiente enlace:

Descargar los VIs Aquí

Eso es todo por la entrada del dia de hoy, espero les haya gustado y hayan aprendido algo nuevo. Si te ha servido el contenido de esta entrada, de los videos y los códigos de implementación y deseas apoyar mi trabajo invitandome a un café super barato, puedes hacerlo en el siguiente link:

👉 Invitar a Sergio a un Café ☕️

Que esten muy bien, nos vemos en la siguiente entrada.

Volver al Curso de Análisis de Sistemas

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.