Práctica 13: Simulación de sistemas dinámicos

Un ente matemático que sirve para representar la relación que existe entre las entradas y las salidas en los sistemas es una ecuación diferencial. A lo largo de nuestra formación académica como ingenieros, estamos en contacto con el término ecuación, cuya definición formal es: Una ecuación se define como una igualdad establecida entre dos expresiones, en la cual puede haber una o más incógnitas o términos cuyo valor numérico es desconocido.

A diferencia de las ecuaciones con incógnitas reales o complejas, las ecuaciones diferenciales tienen como incógnita a una función. Una ecuación diferencial en el contexto de los sistemas dinámicos, involucra a las señales de entrada y salida, así como diversas derivadas de cualquiera de ellas o de ambas. Por ejemplo, la ecuación diferencial:

(12)\[\begin{equation} 2\frac{d^2y(t)}{dt^2}+300y(t)=4.5x(t)+5\frac{dx(t)}{dt} \end{equation}\]

La solución a la ecuación diferencial es la salida del sistema \(y(t)\), es decir, es aquella función que satisface la igualdad establecida en la transformación o estructura (12). En el caso de los sistemas mecánicos \(x(t)\) representará la evolución en el tiempo de una fuerza de entrada, en forma de par o empuje. En sistemas en los que existan almacenes de energía cinética en forma de masas, es común usar también aceleración. En el caso de los sistemas eléctricos, los almacenes de energía son los inductores o bobinas y los capacitores y la entrada son variables como los voltajes o corrientes. Es posible encontrar análogos muy útiles entre los sistemas mecánicos y eléctricos debido a su comportamiento y las relaciones constitutivas que se usan para construir la estructura de la ecuación diferencial que los modelan. Ejemplos de relaciones constitutivas para los sistemas mecánicos, de acuerdo con la ley de Newton son:

(13)\[\begin{equation} m\frac{d^2x(t)}{dt^2}=ma~~F_b=b\frac{dx(t)}{dt}~~F_r=kx(t) \end{equation}\]

Donde \(m\) es la masa en \([kg]\), \(b\) la constante de amortiguador con unidades \([Ns/m]\) y \(k\) en unidades \([N/m]\). Las fuerzas \(F_b\) y \(F_r\) son las fuerzas del amortiguador y la fuerza del resorte respectivamente. Además, de acurdo con los conceptos de la cinemática, se sabe que la aceleración es la segunda derivada de la posición y la velocidad es la razón de cambio o derivada de la posición \(a=\frac{d^2x(t)}{dt}\) y \(v=\frac{dx(t)}{dt}\). La Figura 10 muestra un típico sistema masa resorte amortiguador:

Figura 10 Sistema masa resorte amortiguador.

Aplicando la segunda ley de Newton al sistema masa resorte amortiguador se obtiene la ecuación diferencial o el modelo del sistema masa resorte amortiguador:

(14)\[\begin{equation} m\frac{d^2x(t)}{dt^2}+b_1\frac{dx(t)}{dt}+k_1x(t)=f(t) \end{equation}\]

En este caso, la salida se nombra \(x(t)\) dado que es un desplazamiento causado por la fuerza de entrada \(f(t)\), que sería equivalente a la función estímulo o de entrada, se dice entonces que para este sistema, la entrada es una fuerza y la salida un desplazamiento.

Nota

Es común usar la notación abreviada para las ecuaciones diferenciales:

\[ \frac{d^2x(t)}{dt^2}=\ddot{x}~~\frac{x(t)}{dt}=\dot{x} \]

Usando esta notación se obtiene la forma alternativa de la ecuación diferencial de movimiento o modelo matemático del sistema masa resorte amortiguador como:

\[ m\ddot{x}+b\dot{x}+kx=f \]

En MATLAB®, la solución de las ecuaciones diferenciales se obtiene usando el comando dsolve(edo,ci), donde edo es la ecuación diferencial que se desea resolver y ci son las condiciones iniciales \(x(0)\) y \(\frac{dx}{dt}|_{t=0}\). Por ejemplo, para un sistema masa-resorte-amortiguador con la ecuación diferencial de movimiento:

\[ 2.5\frac{dx^2}{dt}+3\frac{dx}{dt}+350x(t)=0 \]

con condiciones iniciales \(x(0)=0.01\) y \(\frac{dx}{dt}|_{t=0}=0\) se resuelve con los siguientes comandos:

clear
close all
clc
syms x(t);
Dx=diff(x,t);
ode = 2.5*diff(x,t,2)+3*diff(x,t)+350*x==0;
cond1=x(0)==0.01;
cond2=Dx(0)==0;
ci=[cond1 cond2];
xSol(t) =dsolve(ode,ci)

\[xSol(t) =\frac{{\mathrm{e}}^{-\frac{3\,t}{5}} \,{\left(3491\,\cos \left(\frac{\sqrt{3491}\,t}{5}\right)+3\,\sqrt{3491}\,\sin \left(\frac{\sqrt{3491}\,t}{5}\right)\right)}}{349100}\]

La solución de la ecuación diferencial es \(x(t)\) y describe el movimiento de la masa como una función del tiempo. La ecuación diferencial anterior, al estar igualada a cero, es decir, cuando \(f(t)=0\), se denomina ecuación diferencial homogénea, en este tipo de ecuaciones diferenciales, la única manera en la que la solución es diferente de cero, es cuando al menos una de las condiciones iniciales \(x(0)\) y \(\dot{x}(0)\) sea diferente de cero. Un ejemplo de aplicación para cuando \(f(t)\) es diferente de cero, es el caso de la ecuadión diferencial no homogénea, En MATLAB® se plantea y resuelve con las siguientes instrucciones:

clear 
close all
clc

syms x(t);
Dx=diff(x,t);
ode = 2.5*diff(x,t,2)+3*diff(x,t)+350*x==4;
cond1=x(0)==0.0;
cond2=Dx(0)==0.0;
ci=[cond1 cond2];
xSol(t) =dsolve(ode,ci)

\[xSol(t) =\frac{2}{175}-\frac{6\,\sqrt{3491}\,{\mathrm{e}}^{-\frac{3\,t}{5}} \,\sin \left(\frac{\sqrt{3491}\,t}{5}\right)}{610925}-\frac{2\,{\mathrm{e}}^{-\frac{3\,t}{5}} \,\cos \left(\frac{\sqrt{3491}\,t}{5}\right)}{175}\]

La ecuación diferencial que modela a un sistema masa-resorte-amortiguador se puede escribor de forma alternativa si se hace un despeje de la función de movimiento y se escribe de forma alternativa:

\[ \frac{d^2x}{dt^2}= \frac{1}{m}\left ( f(t) - b\frac{dx}{dt}-kx(t)\right ) \]

Si se integra dos veces con respecto al tiempo la ecuación anterior se obtiene:

\[ \int^{t}_0\int^{\tau_1}_0\frac{d^2x}{dt^2}d\tau_1d\tau_2= \int^{t}_0\int^{\tau_1}_0\frac{1}{m}\left ( f(\tau_1) - b\frac{dx(\tau_1)}{d\tau_1}-kx(\tau_1)\right )d\tau_1d\tau_2 \]

se obtiene la siguiente expresión:

\[ x(t)= \frac{1}{m}\int^{t}_0\int^{\tau_1}_0f(\tau_1)d\tau_1d\tau_2-\frac{b}{m}\int^{t}_0\frac{dx(\tau_1)}{d\tau_1}d\tau_1-\frac{k}{m}\int^{t}_0\int^{\tau_1}_0x(\tau_1)d\tau_1d\tau_2 \]

El operador integral en MATLAB® Simulink® es el bloque marcado con la etiqueta \(\frac{1}{s}\), un diagrama de bloques de ejmeplo, para una ecuación diferencial lineal y de coeficientes constantes se muestra en la Figura 11

Figura 11 Diagrama de bloques para la solución de una ecuación diferencial.

Para el caso de un sistema eléctrico del tipo resistencia capacitor inductancia, como el mostrado en la Figura 12

Figura 12 Circuito RLC.

Las relaciones constitutivas del sistema eléctrico son:

(15)\[\begin{equation} v_L=L\frac{di(t)}{dt}~~i_C=C\frac{dv_c(t)}{dt}~~v_R=i(t)R \end{equation}\]

Con \(L\) la inductancia en \([H]\), \(C\) la capacitancia en \([F]\) y \(R\) la resistencia en \([\Omega]\). Además, por la ley de Kirchoff, la corriente en el capacitor \(i_C(t)=i(t)\). En este caso, la salida se nombra \(v_c(t)\) el voltaje en el capacitor, debido a la presencia de la fuente de alimentación \(V_s(t)\), que sería equivalente a la función estímulo o de entrada. La ecuación diferencial que modela a este sistema eléctrico es:

(16)\[\begin{equation} LC\frac{d^2v_c(t)}{dt^2}+RC\frac{dv_c(t)}{dt}+v_c(t)=V_s(t) \end{equation}\]

El diagrama de bloques de esta última ecuación diferencial se realiza de forma similar al del sistema masa-resorte-amortiguador. En la siguiente sección se detalla la manera de realizar la simulación usando diagramas de bloques en el ambiente gráfico Simulink®.

El entorno de simulación gráfica Simulink®

El entorno gráfico de simulación de MATLAB® permite la construcción de modelos de sistemas dinámicos en forma de diagramas de bloques. La representación en diagramas de bloques es común dada su facilidad para representar la interacción entre las variables y los elementos de un sistema. Existen muchas maneras de repesentar el mismo sistema dinámico; en forma de ecuación matricial o en espacio de estados, la representación haciedo uso de la función de transferencia y la reperesentación como ecuación diferencial en la que se presentan los integradores como los bloques básicos de solución al problema de modelado matemático. En este curso no ahondaremos en el tema de los modelos matemáticos, esta práctica tiene por objetivo ser una introducción al ambiente gráfico de simulación de MATLAB® conocido como Simulink®. El entorno gráfico de simulink® se encuentra en el menú de la parte superior de la página principal de MATLAB® como se muestra en la Figura 13 y la Figura 14:

Figura 13 Localización del botón de inicio de simulink®.

Figura 14 Interfaz de inicio de simulink®.

Al hacer click en ese ícono, se abrirá la interfaz gráfica para la creación de un modelo en forma de diragrama de bloques como se muestra en el siguiente video:

La creación de un diagrama de bloques se logra arrastrando los bloques necesarios desde la biblioteca de symulink, que se muestran en la Figura 15

Figura 15 Biblioteca de bloques de simulink®

La solución de la ecuación diferencial es \(x(t)\) y describe el movimiento de la masa como una función del tiempo. La herramienta simulink® permite crear simulaciones en forma de diagramas de bloques como se muestra en los siguientes videos.

La creación del diagrama de bloques y los detalles de la simulación se muestran en el siguiente video:

En el ambiente gráfico de Simulink®, existe un bloque que permite almacenar el resultado de la simulación en un arreglo para poder hacer una gráfica o hacer calculos estadísticos como el promedio, desviación estandar, etc. Los bloques de exportación de datos al workspace se muestra en el siguiente video:

Ejercicio de la práctica 13

1.-Simule el sistema masa resorte amortiguador mostrado en la Figura 10, con los parámetros de la Tabla 15, para los siguientes casos de \(f(t)\):

1. Escalón con amplitud \(10~[N]\) y retardo de \(0.5\) segundos.

2. Escalón con amplitud de \(15~[N]\) sin retardo.

3. Fuerza definida por la función \(f(t)=5sen(4\pi t)~[N]\).

Tabla 15 Tabla de parámetros

\[Parámetro\]	\[Valor\]
\[m\]	\[3.5~[Kg]\]
\[b\]	\[0.15~[Ns/m]\]
\[k\]	\[430~[N/m]\]

Añada los bloques necesarios para la exportación de datos a workspace y genere una gráfica de la fuerza de entrada contra el tiempo y una del desplazamiento contra el tiempo en la misma figura.

2.-Simule el circuito RLC mostrado en la Figura 12, con los parámetros de la Tabla 16, para los siguientes casos de \(V_s(t)\):

1. Escalón con amplitud \(2.5~[V]\) y retardo de \(0.5\) segundos.

2. Escalón con amplitud de \(1.63~[V]\) sin retardo.

3. Voltaje definido por la función \(f(t)=5sen(2\pi t)~[V]\).

Tabla 16 Tabla de parámetros

\[Parámetro\]	\[Valor\]
\[R\]	\[1600~[\Omega]\]
\[L\]	\[15~[H]\]
\[C\]	\[1~[\mu F]\]

Añada los bloques necesarios para la exportación de datos a workspace y genere una gráfica del voltaje de entrada contra el tiempo y una del voltaje en el capacitor contra el tiempo en la misma figura.