notebook.community

Edit and run

¿Cómo funciona la suspensión de un auto?

Una primer aproximación al modelo de la suspensión de un automovil es considerar el oscilador armónico amortiguado.

Referencia:

Un modelo que describe el comportamiento del sistema mecánico anterior es

\begin{equation} m\frac{d^2 x}{dt^2}=-c\frac{dx}{dt}-kx \end{equation}

donde $c$ es la constante de amortiguamiento y $k$ es la constante de elasticidad. Revisar modelado

Documentación de los paquetes que utilizaremos hoy.

En python existe una función llamada _odeint_ del paquete _integrate_ de la libreria _scipy_, que permite integrar sistemas vectoriales de primer orden, del tipo

\begin{equation} \frac{d\boldsymbol{y}}{dt} = \boldsymbol{f}(t,\boldsymbol{y}); \qquad \text{ con }\quad \boldsymbol{y}\in\mathbb{R}^n,\quad \boldsymbol{f}:\mathbb{R}_{+}\times\mathbb{R}^n\to\mathbb{R}^n \end{equation}

con condiciones iniciales $\boldsymbol{y}(0) = \boldsymbol{y}_{0}$. Notar que $\boldsymbol{y}$ representa un vector de $n$ componentes.

Ahora, si nos fijamos bien, el modelo del oscilador armónico amortiguado que obtuvimos es una ecuación diferencial ordinaria (EDO) de segundo orden. No hay problema. La podemos convertir en un sistema de ecuaciones de primer orden de la siguiente manera:

  1. Seleccionamos el vector $\boldsymbol{y}=\left[y_1\quad y_2\right]^T$, con $y_1=x$ y $y_2=\frac{dx}{dt}$.
  2. Notamos que $\frac{dy_1}{dt}=\frac{dx}{dt}=y_2$ y $\frac{dy_2}{dt}=\frac{d^2x}{dt^2}=-\frac{c}{m}\frac{dx}{dt}-\frac{k}{m}x=-\frac{c}{m}y_2-\frac{k}{m}y_1$.
  3. Entonces, el modelo de segundo orden lo podemos representar como el siguiente sistema vectorial de primer orden:
\begin{equation} \frac{d\boldsymbol{y}}{dt}=\left[\begin{array}{c}\frac{dy_1}{dt} \\ \frac{dy_2}{dt}\end{array}\right]=\left[\begin{array}{c}y_2 \\ -\frac{k}{m}y_1-\frac{c}{m}y_2\end{array}\right]=\left[\begin{array}{cc}0 & 1 \\-\frac{k}{m} & -\frac{c}{m}\end{array}\right]\boldsymbol{y}. \end{equation} 




In [1]:



    


# Primero importamos todas las librerias, paquetes y/o funciones que vamos a utlizar
import matplotlib.pyplot as plt         # ¿Para qué sirve?
import numpy as np                      # ¿Para qué sirve?
from scipy.integrate import odeint      # Función para integrar EDO vectoriales de primer orden
%matplotlib inline                      # ¿Para qué sirve?



    


















    






UsageError: unrecognized arguments: # ¿Para qué sirve?

















In [2]:



    


# Definimos los parámetros
k = 3.0                                         # Constante del muelle
m = 1.0                                         # Masa
B = .5                                          # Constante de amortiguación

# Función f(t,y) que vamos a integrar
def armonico(y, t):
    y1, y2 = y
    return [y2, -k * y1 / m - B / m * y2]

y0 = [0.6, 0.4]                                 # Vector de posición inicial y velocidad inicial
                                                # condiciones iniciales y1(0)=0.6 [m] y2(0)=0.4 [m/s]

# Especificamos los puntos de tiempo donde queremos la solución
#tiempo = np.linspace(0, 15)                    # Dominio temporal de 0 a 15 (por defecto 50 puntos)
tt = np.arange(0, 20, .01)                       # Dominio temporal de 0 a 20 con paso de 0.01

yy = odeint(armonico, y0, tt)
#El sistema se resuelve  con 
#odeint(funcion_sistema, condiciones_iniciales, dominio_tiempo)



    











In [3]:



    


# Veamos como nos entrega las soluciones la funcion odeint
np.shape(yy)                                    # shape es una función de numpy que nos entrega las dimensiones del objeto



    


















    
Out[3]:








(2000, 2)

















In [4]:



    


yy



    


















    
Out[4]:








array([[ 0.6       ,  0.4       ],
       [ 0.60389997,  0.37999104],
       [ 0.60759977,  0.35996814],
       ..., 
       [-0.00314526, -0.00501178],
       [-0.00319478, -0.00489192],
       [-0.0032431 , -0.00477119]])
  • $yy$ es una matriz de 2000 filas y 2 columnas.
  • La primer columna de $yy$ corresponde a $y_1$.
  • La segunda columna de $yy$ corresponde a $y_2$.

¿Cómo extraemos los resultados $y_1$ y $y_2$ independientemente?





In [5]:



    


## Una forma
#yy1 = yy[:,0]                                   # tome de yy todas las filas (:) y la primer columna
#yy2 = yy[:,1]                                   # tome de yy todas las filas (:) y la segunda columna

# Otra forma (más compacta)
yy1, yy2 = yy.T                                 # transponemos la matriz yy y extraemos los resultados.



    











In [6]:



    


# Comandos para hacer los números de los ejes más grandes
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size 
mpl.rcParams['ytick.labelsize'] = label_size



    











In [7]:



    


# Graficamos
plt.plot(tt, yy1, c = 'r', label="Posicion") 
plt.plot(tt, yy2, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('$t$', fontsize = 14)
plt.ylabel('$y_1$, $y_2$', fontsize = 14)
plt.show()



    


















    
























In [8]:



    


# Frecuencia natural
omega0 = k/m
# Graficamos el espacio fase
plt.figure(figsize = (6,6))
plt.scatter(yy1, yy2/omega0, lw = 0, s = 3, cmap = 'viridis', c = yy1)
plt.show()

 Dependiendo de los parámetros, 3 tipos de soluciones

Teníamos \begin{equation} m\frac{d^2 x}{dt^2} + c\frac{dx}{dt} + kx = 0 \end{equation}

si recordamos que $\omega_0 ^2 = \frac{k}{m}$ y definimos $\frac{c}{m}\equiv 2\Gamma$, tendremos

\begin{equation} \frac{d^2 x}{dt^2} + 2\Gamma \frac{dx}{dt}+ \omega_0^2 x = 0 \end{equation}

Subamortiguado

Si $\omega_0^2 > \Gamma^2$ se tiene movimiento oscilatorio subamortiguado.





In [13]:



    


omega0 = k/m
Gamma = B/(2*m)



    











In [14]:



    


omega0**2, Gamma**2



    


















    
Out[14]:








(9.0, 0.0625)

















In [15]:



    


omega0**2 > Gamma**2



    


















    
Out[15]:








True

Entonces, el primer caso que ya habíamos presentado corresponde a movimiento amortiguado.





In [18]:



    


plt.figure()
plt.plot(tt, yy1, c = 'r', label="Posicion") 
plt.plot(tt, yy2, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('t', fontsize = 14) 
plt.show()

Sobreamortiguado

Si $\omega_0^2 < \Gamma^2$ se tiene movimiento oscilatorio sobreamortiguado.





In [19]:



    


k = .1 # Constante del muelle
m = 1.0   # Masa
B = .5 # Constante de amortiguación



    











In [20]:



    


omega0 = k/m
Gamma = B/(2*m)



    











In [21]:



    


omega0**2, Gamma**2



    


















    
Out[21]:








(0.010000000000000002, 0.0625)

















In [22]:



    


omega0**2 < Gamma**2



    


















    
Out[22]:








True

















In [24]:



    


y0 = [0.6, 0.4] 
tt = np.arange(0, 20, .01)
yy = odeint(armonico, y0, tt)



    











In [25]:



    


yy1s, yy2s = yy.T  # extraer posición y velocidad.



    











In [26]:



    


plt.plot(tt, yy1s, c = 'r', label="Posicion") 
plt.plot(tt, yy2s, c = 'k', label="Velocidad")
plt.legend(loc = 'best', prop={'size': 14})    
plt.xlabel('tiempo', fontsize = 14) 
plt.show()

Amortiguamiento crítico

Si $\omega_0^2 = \Gamma^2$ se tiene movimiento críticamente amortiguado.





In [27]:



    


k = np.sqrt(.0625) # Constante del muelle
m = 1.0   # Masa
B = .5 # Constante de amortiguación



    











In [28]:



    


omega0 = k/m
Gamma = B/(2*m)



    











In [29]:



    


omega0**2, Gamma**2



    


















    
Out[29]:








(0.0625, 0.0625)

















In [30]:



    


omega0**2 == Gamma**2



    


















    
Out[30]:








True

















In [31]:



    


y0 = [0.6, 0.4] 
tt = np.arange(0, 20, .01)
yy = odeint(armonico, y0, tt)



    











In [33]:



    


yy1c, yy2c = yy.T  # extraer posición y velocidad.



    











In [35]:



    


plt.plot(tt, yy1c, c = 'r', label="Posicion") 
plt.plot(tt, yy2c, c = 'k', label="Velocidad")
plt.legend(loc = 'best',prop={'size': 14})  
plt.xlabel('tiempo', fontsize = 14) 
plt.show()

En resumen, se tiene entonces:





In [36]:



    


fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col',
                                                       sharey='row',figsize =(10,6))
ax1.plot(tt, yy1, c = 'k') 
ax1.set_title('Amortiguado', fontsize = 14)
ax1.set_ylabel('Posición', fontsize = 14)
ax2.plot(tt, yy1s, c = 'b') 
ax2.set_title('Sobreamortiguado', fontsize = 14)
ax3.plot(tt, yy1c, c = 'r') 
ax3.set_title('Crítico', fontsize = 16)
ax4.plot(tt, yy2, c = 'k')
ax4.set_ylabel('Velocidad', fontsize = 14)
ax4.set_xlabel('tiempo', fontsize = 14)
ax5.plot(tt, yy2s, c = 'b')
ax5.set_xlabel('tiempo', fontsize = 14)
ax6.plot(tt, yy2c, c = 'r')
ax6.set_xlabel('tiempo', fontsize = 14)
plt.show()

Actividad. ¿Cómo se ve el espacio fase para los diferentes casos así como para diferentes condiciones iniciales?

Created with Jupyter by Lázaro Alonso. Modified by Esteban Jiménez Rodríguez. 




In [ ]:

Content source: Jesusomar97/Simulacion2017

Similar notebooks:

notebook.community | gallery | about