Descripción del campo fluidos Asignatura: Mecánica de Fluidos

Departamento: Ciencia y Tecnología de Materiales y Fluidos

Centro: Escuela Universitaria Politécnica de Teruel

Profesor: Adrián Navas Montilla

Movimiento armónico simple de un sistema cilindro-pistón isotermo

Este cuaderno ha sido desarrollado por una estudiante 2º curso del Grado en Ingieniería Electrónica y Automática (EUPT-UNIZAR) como ejercicio voluntario en el marco del PIIDUZ 5494 "Python en ciencias e ingeniería: creación colaborativa"

Autora: Nerea Pilar Yagüe Lorente

Profesor: Adrián Navas Montilla

Asignatura: Temodinámica Técnica y Fundamentos de Transmisión de Calor

Enunciado propuesto:

Considere un dispositivo cilindro pistón en equilibrio con volumen inicial $V_0$, área transversal $A$ y masa del pistón $m_p$, que encierra una masa $m_a$ de aire a temperatura constante, $T$. En el otro lado del pistón, actúa la presión atmosférica. Para pequeños desplazamientos del pistón (respecto a $V_0$/$A$) alrededor de su posición de equilibrio, se puede demostrar que su movimiento sigue un movimiento armónico simple. Utilice la segunda ley de Newton para obtener la ecuación de movimiento del pistón, asumiendo que la longitud total del sistema es mucho mayor que la amplitud de los desplazamientos y considerando que el sistema está en posición horizontal. Prepare un Notebook de Python en el que describa el problema y lo resuelva siguiendo la plantilla proporcionada. Debe incluir una breve derivación explicada de la ecuación de movimiento utilizando sintaxis Latex y Markdown, así como un diagrama de sólido libre del sistema indicando las fuerzas.

Ayuda: Un sistema dinámico de segundo orden no amortiguado viene dado por la ecuación:

$$ \ddot{x} + \omega^2 x =0$$

donde $\omega$ es la frecuencia angular del sistema. Nótese que $x$ es la posición respecto al punto de equilibrio, $dx/dt=\dot{x}$ es la velocidad y $d^2x/dt^2=\ddot{x}$ es la aceleración del sistema.

La ecuación diferencial anterior es de segundo orden y se puede escribir como un sistema de 2 ecuaciones de primer orden:

$$ \frac{d}{dt}X = \left( \begin{array}{c} -\omega^2 x \\ \dot{x} \\ \end{array} \right) $$

con

$$ X = \left( \begin{array}{c} \dot{x} \\ x \end{array} \right)$$

Derivación de la ecuación del movimiento cilindro-pistón

A continuación se va a explicar la derivación de la ecuación del movimiento para un sistema cilindro-pistón cuyos datos son:

  • Masa del pistón: $ m_p $
  • Masa del aire encerrada: $ m_a $
  • Área del pistón: $ A $
  • Volumen inicial: $ V_{0} $
  • Presión atmosférica: $ P_0 $
  • Temperatura constante: $ T $
  • Longitud $ L=\frac{V_{0}}{A} $

La fuerza neta que se ejerce sobre el pistón es igual a la suma de la fuerza que ejerce el gas encerrado más la fuerza que ejerce la presión atmostérica. Dado que $P=\frac{F}{A}$, entonces $F=P\cdot A$. Además, al ejercerse en la misma dirección pero en sentidos contrarios, la fuerza que ejerce la presión atmosférica se pone con signo negativo, dando lugar a:

$$F_{neta} = P(x) \cdot A - P_0 \cdot A$$

Tomando el aire como gas ideal, se puede utilizar la ecuación de los gases ideales ($PV = m_a R T$ con R la constante de los gases ideales para el aire seco) para hallar la presión interna $P(x)$:

$$P(x) = \frac{m_a R T}{V}$$

Si la longitud del cilindro es: $L=\frac{V_{0}}{A}$, el $V$ para cualquier instante será: $V=(L+x)\cdot A$

Sustituyendo, se tiene que:

$$P(x) = \frac{m_a R T}{(L+x)\cdot A}$$

En la situación de equilibrio ($x=0$) el pistón está en reposo, es decir, la fuerza interna iguala a la fuerza atmosférica, por lo que:

$$P_0 = \frac{m_a R T}{L\cdot A}$$

Llevando estas expresiones a la fuerza neta:

$$F_{neta} = \frac{m_a R T}{(L+x)\cdot A} \cdot A - \frac{m_a R T}{L\cdot A} \cdot A = \frac{m_a R T}{L+x} - \frac{m_a R T}{L}=m_a R T \left(\frac{1}{L+x}-\frac{1}{L}\right)=-m_a R T\frac{x}{L(L+x)}$$

Cuando el pistón se mueve, la fuerza neta no depende directamente de la presión atmosférica $P_0$, sino de la diferencia en la presión entre el gas y la atmósfera. Para pequeños desplazamientos $x$, la diferencia en la presión se vuelve pequeña y aproximadamente lineal. En otras palabras, la fuerza neta es aproximadamente proporcional a $x$ cuando $x$ es pequeño.

Esto indica que:

$$F_{neta} \approx -\frac{m_a R T}{L^2} x$$

Aplicando la segunda Ley de Newton ($\sum F=m_p·\ddot{x}$) se obtiene que:

$$-\frac{m_a R T}{L^2} x = m_p \ddot{x}$$

Por tanto, la ecuación del movimiento es la ecuación diferencial del oscilador armónico simple: $$\ddot{x} + \omega^2 x = 0$$

donde la frecuencia angular es: $$\omega^2 = \frac{m_a R T}{m_p L^2}$$ $$\omega=\sqrt\frac{m_a R T}{m_p L^2}=\sqrt\frac{m_a R T}{m_p}·\frac{1}{L}$$

Siendo R la constante de los gases ideales para el aire seco, un valor conocido de $R=287\frac{J}{kgK}$.

Diagrama de cuerpo libre

A continuación se muestra el diagrama de cuerpo libre del pistón con las fuerzas involucradas:

  • $ P(x)A $: presión del gas actuando hacia la izquierda.
  • $ P_0A $: presión atmosférica actuando hacia la derecha.
In [ ]:
# Representación del diagrama de cuerpo libre
fig, ax = plt.subplots(figsize=(8, 3))  # se abre una figura vacía de dimesiones 8·3

# Dibujo del pistón (se ha necesitado importar la librería patches para añadir los rectángulos)
#  para representar la masa del aire:
ax.add_patch(patches.Rectangle((3.9, 1.25), 1, 0.5, edgecolor='k', facecolor='lightgreen'))
#  para representar el pistón:
ax.add_patch(patches.Rectangle((3.7, 1.25), 0.2, 0.5, edgecolor='k', facecolor='k'))
ax.add_patch(patches.Rectangle((2.7, 1.43), 1, 0.13, edgecolor='k', facecolor='k'))
ax.text(3.6, 1.9, 'PISTÓN')

# Flecha hacia la derecha (presión atmosférica)
ax.arrow(3, 1.65, 0.6, 0, head_width=0.1, head_length=0.1, fc='b', ec='b')  # coordenadas x y de inicio, vector dirección de la flecha, ancho de la punta, largo de la punta y colores
ax.text(2.9, 1.75, r'$P_0 A$', color='b') # coordenadas del texto de la flecha, texto y color

# Flecha hacia la izquierda (presión del gas)
ax.arrow(4.62, 1.5, -0.6, 0, head_width=0.1, head_length=0.1, fc='r', ec='r')  # coordenadas x y de inicio, vector dirección de la flecha, ancho de la punta, largo de la punta y colores
ax.text(4.5, 1.6, r'$P(x) A$', color='r') # coordenadas del texto de la flecha, texto y color

# Se establecen unos límites para que la imagen quede centrada
ax.set_xlim(1.5, 6)
ax.set_ylim(0.8, 2.2)

# Se quitan los ejes
ax.axis('off')

ax.set_title('DIAGRAMA DE CUERPO LIBRE DEL PISTÓN')
plt.show()

EJEMPLO 1: aprender a dibujar una gráfica de la posición frente al tiempo de un movimiento armónico simple para un valor de omega dado

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from scipy.integrate import odeint
In [ ]:
def rhs(X, t, omega): #Esto es el lado derecho de la ecuación, que equivale a la derivada temporal de X
    xdot, x = X
    return [-omega**2 * x, xdot] #omega es velocidad angular, x es posición y xdot es velocidad

Vamos a definir un vector de tiempos en el que resolver la ecuación:

In [ ]:
dt = 0.1 #paso de tiempo para la integración temporal, cuanto más pequeño mejor resolución
tiempo_final=12.0 #tiempo final para resolver el sistema
ts = np.arange(0, tiempo_final, dt)

Y a añadir las condiciones iniciales:

In [ ]:
X0 = [0, 1] #donde la primera componente es xdot y la segunda x: en este caso hacemos que comience con velocidad cero y posición x=1

Y finalmente se resuelve la ecuación utilizando un integrador temporal predefinido odeint

In [ ]:
omega = 2
scipysol = odeint(rhs, X0, ts, args=(omega,)) #esta función sirve para obtener la solución mediante integración numérica (como hace Matlab y Simulink) en los instantes definidos en el vector ts.
#La solución se almacena en scipysol, que será un array de 2 dimensiones
plt.plot(ts, scipysol[:, 1]); #esto realiza una representación de la segunda componente del vector X, es decir, de la posición, frente al tiempo

EJEMPLO 2: aprender, a partir de los valores de los parámetros físicos, a calcular omega y dibujar con ella una gráfica de la posición frente al tiempo de un movimiento armónico simple

In [ ]:
# Se estiman valores para los parámetros físicos
mp = 0.1       # kg
ma = 0.01      # kg
A = 0.01       # m^2
V0 = 0.001     # m^3
R = 287        # J/kg·K
T = 300        # K
L= 0.5         # m        Se elige con idea de que sea mucho mayor que los desplazamientos de x

# Cálculo de omega
omega = np.sqrt((ma * R * T) / (mp * L**2))
print(f"Frecuencia angular ω = {omega:.2f} rad/s")

# Se define la ecuación
def rhs(X, t, omega): #Esto es el lado derecho de la ecuación, que equivale a la derivada temporal de X
    xdot, x = X
    return [-omega**2 * x, xdot] #omega es velocidad angular, x es posición y xdot es velocidad

# Se define el vector de tiempos
dt = 0.1 #paso de tiempo para la integración temporal, cuanto más pequeño mejor resolución
tiempo_final=12.0 #tiempo final para resolver el sistema
ts = np.arange(0, tiempo_final, dt)

# Condiciones iniciales: velocidad 0, desplazamiento 1 mm (para que sea lo suficientemente más pequeño que L)
X0 = [0, 0.001] #donde la primera componente es xdot y la segunda x: en este caso hacemos que comience con velocidad cero y posición x=0.001m
scipysol = odeint(rhs, X0, ts, args=(omega,)) #esta función sirve para obtener la solución mediante integración numérica (como hace Matlab y Simulink) en los instantes definidos en el vector ts.
#La solución se almacena en scipysol, que será un array de 2 dimensiones

# Graficar desplazamiento
plt.plot(ts, scipysol[:, 1]*1000)  #esto realiza una representación de la segunda componente del vector X, es decir, de la posición, frente al tiempo. Se multiplica por 1000 para representarlo en mm.
plt.title("Desplazamiento del pistón (Movimiento Armónico Simple)")
plt.xlabel("Tiempo [s]")
plt.ylabel("Desplazamiento [mm]")
plt.grid()

Frecuencia angular ω = 185.58 rad/s

Conclusión

El pistón realiza un movimiento armónico simple para pequeños desplazamientos. Su frecuencia angular depende de:

  • Masa del aire encerrada $m_a$
  • Constante de los gases ideales para el aire seco $R=287\frac{J}{kgK}$
  • Temperatura del aire $T=cte$
  • Masa del pistón $m_p$
  • Longitud del pistón $L$

El análisis es válido bajo la suposición de condiciones isotérmicas y amplitudes pequeñas.

Bibliografía

Para llegar a la ecuación del movimiento:

Pistón oscilante. (s. f.). http://tesla.us.es/wiki/index.php/Pist%C3%B3n_oscilante (Este enlace proporciona el estudio del pistón en vertical)

Para representar el diagrama de cuerpo libre:

matplotlib.patches.Rectangle — Matplotlib 3.10.1 documentation. (s. f.). https://matplotlib.org/stable/api/_as_gen/matplotlib.patches.Rectangle.html

matplotlib.patches.Arrow — Matplotlib 3.10.1 documentation. (s. f.). https://matplotlib.org/stable/api/_as_gen/matplotlib.patches.Arrow.html