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
Resolución de un ejemplo práctico¶
Considere una instalación de bombeo de agua entre dos depósitos abiertos a la atmósfera con una diferencia de cota geométrica entre ellos de 36 metros. Los depósitos se conectan a través de una lÃnea compuesta por dos tramos de 5 m (entre el depósito 1 y la bomba) y 500 m (entre la bomba y el depósito 2) respectivamente. Se dispone de un conducto de sección circular ($D=0.1$ m), para el que asumiremos un factor de fricción constante de $f=0.02$. Se consideran nulas las pérdidas singulares. Para esta instalación se ha escogido una bomba cuyas curvas caracterÃsticas proporcionadas por el fabricante son $H_{B}(Q)=43-40000Q^2$ y $\eta(Q)=180Q-10000Q^2$ (unidades en m y m$^3$/s).
from sympy import * # LibrerÃa para trabajo simbólico
import numpy as np # LibrerÃa para cálculo numérico
import math # LibrerÃa para utilizar sÃmbolos matemáticos como el número pi, que se escribe como math.pi
import matplotlib.pyplot as plt # LibrerÃa para poder dibujar gráficas
from ipywidgets import interact, interactive, fixed, interact_manual
Para calcular en punto de operación igualaremos la curva de la bomba a la curva de la instalación
$$H_B(Q)=H_I(Q)$$es decir
$$43-40000Q^2=z_2-z_1 + \frac{fL}{D} \frac{Q^2}{2gS^2}$$y resolveremos la ecuación.
En Python, definiremos las funciones h_inst(Q) y h_bomba(Q) y resolveremos el punto de operación utilizando la función solve():
##Datos del problema
g=9.81 #aceleración de la gravedad (m/s^2)
f=0.02 #factor de fricción de Darcy
L=5+500 #longitud total conductos (m)
D=0.1 #diámetro conductos (m)
z1=0 #cota depósito 1 (m)
z2=36 #cota depósito 2 (m)
rho=1000 #densidad agua (kg/m^3)
#calculamos la sección de los conductos
S=math.pi*D**2/4
Q=symbols('Q')
def h_inst(Q): #curva de la instalación
return z2-z1+f*L/(D*2*g)*(Q/S)**2
def h_bomba(Q): #curva de la bomba (dato)
return 43-40000*Q**2
#Punto de operación
sol=solve(h_inst(Q)-h_bomba(Q),Q) #para ello utilizamos la función solve(ecuacion,variable), que permite resolver ecuaciones algebraicas
Qp=sol[1] #la ecuación anterior es cuadrática y tiene 2 soluciones, nos quedamos con la segunda (la positiva)
hp=h_inst(Qp) #evaluamos la curva de la instalación (o de la bomba) en Qp
print("La curva resistente de la instalación es:",h_inst(Q))
print("La curva de la bomba es:",h_bomba(Q))
print("Las soluciones de la ecuación son Q =",sol)
print("El punto de operación es: Qp=",Qp,"m^3/s y hp=",hp,"m")
Ahora hacemos una representación gráfica de las curvas caracterÃsticas y del punto de operación:
Lt=0.01 #longitud del eje X (rango de caudales en m^3/s)
N = 100 #numero de puntos a representar
xp = np.linspace(0, Lt, N) #puntos en x
yp1 = h_inst(xp) #puntos en y: curva instalacion
yp2 = h_bomba(xp) #puntos en y: curva bomba
fig, ax = plt.subplots(figsize=(10,7)) #genera el objeto "fig"
ax.plot(xp,yp1,label='Curva resistente de la instalación')
ax.plot(xp,yp2,label='Curva caracterÃstica de la bomba')
ax.plot(Qp,hp,'o')
ax.plot([0,Qp],[hp,hp],'g--')
ax.plot([Qp,Qp],[np.min(yp1),hp],'g--')
ax.set_title('Curvas bomba/instalación')
ax.set_xlabel("$Q(m^3/s)$")
ax.set_ylabel("$H(m)$")
ax.legend()
Otro cálculo de interés ingenieril es la potencia a suministrar a la bomba (potencia en el eje). Para ello, primero evaluaremos el rendimiento para el caudal obtenido (punto de operación)
$$\eta(Q_p)=180Q-10000Q_p^2$$y después hallaremos la potencia en el eje utilizando la definición de rendimiento
$$\eta(Q_p)=\frac{\rho g Q_p H_p}{\dot{W}_{eje}}$$def rendimiento(Q):
return 180*Q-10000*Q**2
rto = rendimiento(Qp)
potencia_fluido = rho*g*hp*Qp
potencia_eje = potencia_fluido/rto
print("El rendimiento es:",rto)
print("Potencia fluido:",potencia_fluido, "W")
print("Potencia eje:",potencia_eje," W")
Y llegados a este punto uno uede preguntarse, ¿es esta bomba adecuada para la instalación?. Para saber si hemos elegido bien una bomba, debemos examinar si el punto de operación está cerca del punto de máxima eficiencia (best efficiency point, BEP). Para ello, podemos calcular este punto numéricamente buscando el máximo de la curva de rendimiento:
sol=solve(diff(rendimiento(Q),Q))
Qmax=sol[0]
rtomax=rendimiento(Qmax)
print("El caudal para el cual el rendimiento es máximo es:",float(Qmax), "m^3/s")
print("El rendimiento máximo es:",float(rtomax))
fig, ax = plt.subplots(figsize=(10,7))
Lt=0.01
N = 100
xp = np.linspace(0, Lt, N)
yp1 = h_inst(xp)
yp2 = h_bomba(xp)
yp3 = rendimiento(xp)
ax.plot([Qp, Qp],[35, 44.8],'--',color='tab:gray')
ax.plot([Qmax, Qmax],[35, 45],'--',color='tab:gray')
ax.plot(xp,yp1,label='Curva de la instalación')
ax.plot(xp,yp2,label='Curva de la bomba')
ax.plot(Qp,hp,'o',label='Punto de operación actual')
ax.plot(Qmax,h_bomba(Qmax),'*',ms=10,color="tab:green",label='Punto de operación nominal')
ax2 = ax.twinx()
color = 'tab:red'
ax2.set_ylabel('$\eta$', color=color)
ax2.plot(xp,yp3, color=color,label='Rendimiento')
ax2.plot(Qp,rto,'o', color=color)
ax2.plot(Qmax,rtomax,'*',ms=10, color=color)
ax2.tick_params(axis='y', labelcolor=color)
ax.set_title('Curvas bomba/instalación')
ax.set_xlabel("$Q(m^3/s)$")
ax.set_ylabel("$H(m)$")
ax.legend(loc='lower right')
ax2.legend()
A la vista del resultado obtenido y de la gráfica inferior, se observa que el punto de operación actual está cerca del punto de operación nominal (de máximo rendimiento) y podrÃamos afirmar que la bomba es adecuada. Estamos trabajando con un rendimiento del 78% y el rendimiento máximo es 81%
Las instalaciones de fluidos pueden requerir de regulación por diversos motivos: para ajustar el caudal, para ajustar el punto de operación e incrementar la eficiencia, etc. Las estrategias más habituales de regulación de caudal son la regulación por estrangulamiento (cierre de una válvula) a velocidad de giro constante, la regulación por variación de la velocidad de giro de la bomba, y la regulación por variación del ángulo de los álabes en el distribuidor en el rodete. En este ejemplo vamos a considerar la primera de estas estrategias.
Consideremos que colocamos una válvula de estrangulamiento en el tramo de impulsión, cuya constante de pérdidas menores la caracterizamos por $k_{v}$. ¿Cómo se modificarÃa el punto de operación al cerrar la válvula?
En este caso la curva resistente de la instalación vendrá dada por:
$$H_I=z_2-z_1+ \left( \frac{fL}{D} +k_{v} \right) \frac{Q^2}{2gS^2}$$donde observamos que al incrementar el valor de $k_{v}$, la curvatura de la curva aumenta, desplazando el punto de operación hacia la izquierda, es decir, reduciendo el caudal. Esto se aprecia en la siguiente figura.
k_valvula=0
def h_inst2(Q):
return z2-z1+(f*L/D+k_valvula)/(2*g)*(Q/S)**2
#Obtenemos el punto de operación
sol=solve(h_inst2(Q)-h_bomba(Q),Q)
Qp=sol[1]
hp=h_bomba(Qp)
print("La nueva curva de la instalación es:",h_inst2(Q))
print("La curva de la bomba es:",h_bomba(Q))
print("El punto de operación es: Qp=",Qp,"m^3/s y hp=",hp,"m")
LL=0.01
N = 100
xp = np.linspace(0, LL, N)
yp1 = h_inst(xp)
yp2 = h_bomba(xp)
yp3 = h_inst2(xp)
fig, ax = plt.subplots(figsize=(10,7))
ax.plot(xp,yp1,'tab:blue',label='Sin válvula ($k_v=0$)',lw=3)
ax.plot(xp,yp2,'tab:orange')
ax.plot(Qp,hp,'o')
for i in range(1,14,1): #hacemos un bucle para pintar varias curvas
k_valvula=20*i**2
yp3 = h_inst2(xp)
ax.plot(xp,yp3,'--',color='tab:blue')
sol=solve(h_inst2(Q)-h_bomba(Q),Q)
Qp=sol[1]
hp=h_bomba(Qp)
ax.plot(Qp,hp,'o',color='tab:blue')
ax.set_title('Curvas bomba/instalación')
ax.set_xlabel("$Q(m^3/s)$")
ax.set_ylabel("$H(m)$")
ax.annotate("", xy=(0.004, 41), xytext=(0.006, 40),arrowprops={'arrowstyle':'->','lw': 3, 'color': 'tab:green', 'alpha': 1.0})
plt.ylim([min(yp1), max(yp1)])
ax.legend()
