Los
términos HP (Q) y HL (Q) representan la cabeza agregada al sistema por
una bomba, y la cabeza perdida en el sistema de efectos del mundo real
(como fricción, viscosidad, etc.) respectivamente. Tenga en cuenta que
ambos términos son funciones de la tasa de flujo de fluido del sistema,
Q. (Como consecuencia interesante del párrafo anterior que describe la
interpretación de la cabeza, la cabeza de la bomba le dice cuán alta
podría una bomba teóricamente empujar un fluido). Examinaremos los
términos de bomba y pérdida más a fondo, pero antes de que lo hagamos,
simplifiquemos la ecuación anterior para nuestro problema específico.
Mirando
el sistema anterior nuevamente, elegiremos convenientemente nuestras
dos ubicaciones para la ecuación de Bernoulli, de modo que la mayoría de
los términos se cancelen. Podemos hacer esto eligiendo las ubicaciones 1
y 2 para estar en la superficie libre del agua de cada tanque,
respectivamente, donde la presión es constante e igual a la presión
atmosférica (P1 = P2), y la velocidad es aproximadamente constante y
cero (v1 = V2 = 0). También asumiremos que la altura del agua en los dos
tanques es la misma en el instante en que estamos analizando el sistema
de tal manera que Z1 = Z2. Después de simplificar el álgebra, vemos que
casi todos los términos se cancelan, y nos queda el hecho de que la
cabeza producida por la bomba debe igualar la cabeza perdida en el
sistema debido a las no idealidades. Dicho de otra manera, la bomba está
compensando cualquier pérdida de energía en el sistema.
Esta situación se puede ver
cualitativamente en la figura a continuación. La cabeza producida por
una bomba disminuye al aumentar la velocidad de flujo, mientras que las
pérdidas en un sistema de tuberías aumentan al aumentar la velocidad de
flujo. El punto donde las dos curvas se cruzan (cabeza de la bomba =
pérdida de cabeza) determina el punto de operación (caudal) del sistema.
El último paso antes de que
podamos saltar al código es plantear el problema como un problema de
búsqueda de raíces. Restando el lado derecho de la ecuación del lado
izquierdo, obtenemos el problema de resolución de raíces que estamos
buscando. Es decir, hemos planteado nuestro problema de la siguiente
manera: encuentre la velocidad de flujo (q) de modo que el lado
izquierdo de la ecuación a continuación sea igual a cero. En este punto,
el cabezal de la bomba igualará las pérdidas de la cabeza del sistema.
El código
Para
evitar perder el panorama general de lo que estamos haciendo, no voy a
explicar cada pequeño detalle del código (supongo que ya tiene un fondo
razonable en Python). En cambio, voy a enfocar mis esfuerzos en
garantizar que la narración y la estructuración del código sea clara, y
entrar en más detalles según sea necesario. Como siempre, siéntase libre
de hacer cualquier pregunta si algo no está claro.
Configuración
Comenzaremos
importando todos los módulos necesarios. Se volverá evidente cómo cada
uno de los módulos se usa más adelante, pero vale la pena señalar que
las declaraciones de importación clave son las de SciPy. Estas son las
funciones que son específicas del problema en cuestión. El bloque de
código también establece algunas configuraciones de trazado
predeterminadas (al gusto personal), crea una carpeta para guardar las
figuras generadas y define algunas constantes de conversión de unidades
que facilitan nuestra vida más adelante en el código.
>>from dataclasses import dataclass
>>from pathlib import Path
>>import numpy as np
>>import matplotlib.pyplot as plt
>>import pandas as pd
>>from scipy.interpolate import interp1d
>>from scipy.optimize import root_scalar
>>plt.style.use('seaborn-v0_8-darkgrid')
>>plt.rcParams['font.family'] = 'Times New Roman'
>>plt.rcParams['font.size'] = 12
>>figsize = (6.4,4)
>>plots_folder = Path('plots')
>>plots_folder.mkdir(exist_ok=True)
>>INCHES_TO_METERS = 25.4/1000
>>FEET_TO_METERS = 12*INCHES_TO_METERS
>>GALLONS_TO_M3 = 0.0037854118
A continuación haremos una pitón dataclass
que esencialmente actúa como un contenedor para almacenar propiedades
de fluidos (densidad, viscosidad y gravedad). Por defecto, las
propiedades se establecen en las de agua. Tenga en cuenta que, aunque no
es completamente necesario, las clases de datos de Python son
increíblemente convenientes. Si son nuevos para usted, le recomiendo que
consulte este video .
>>@dataclass
>>class Fluid():
>>
>> rho: float = 997
>> mu: float = 0.0007972
>> g: float = 9.81
Modelo de tubería
El
siguiente paso es modelar las pérdidas de cabeza de la tubería (el
término HL (Q) en la ecuación de Bernoulli extendida anterior). Esto
generalmente se hace utilizando la ecuación de Darcy-Weisbach que se
muestra a continuación, donde F es un factor de fricción (más sobre esto
en breve), V es la velocidad de flujo, G es la gravedad, y L y D son la
longitud y el diámetro de la tubería, respectivamente.
Desafortunadamente, el factor de
fricción (F) no es constante, pero también depende de la velocidad del
flujo, las propiedades del fluido y las dimensiones de la tubería.
Existen varios modelos para calcular F, pero usaremos la ecuación de
Haaland, que se muestra a continuación.
En esta ecuación, Epsilon es la
rugosidad de la superficie de la tubería (que se puede encontrar en las
tablas de libros de texto de ingeniería), y RE es el famoso número de
Reynolds, calculado como a continuación.
Por último, podemos tener en
cuenta que el volumen barrido por unidad de tiempo o la velocidad de
flujo volumétrico (Q) es igual al área de sección transversal (a) de la
tubería de la velocidad del flujo (V). Por lo tanto, dada una velocidad
de flujo en la tubería, podemos calcular la velocidad del flujo
correspondiente en la tubería como:
Esperemos
que todas estas ecuaciones no estén rentando el panorama general, solo
estamos buscando un modelo particular de pérdida de la cabeza en una
tubería. Dada una tasa de flujo y las dimensiones de la tubería, primero
calcule la velocidad del flujo correspondiente, luego conéctelo a
través de las ecuaciones anteriores para calcular la pérdida de cabeza
de la tubería. Esto es exactamente lo que el Pipe La clase (que se muestra a continuación) implementa.
El
método de inicialización almacena las dimensiones de las tuberías
(todas que se supone que están en metros) y propiedades de fluido. El A El método calcula el área de sección transversal de la tubería (para aquellos que no están familiarizados con el @property decorador, este artículo lo explica muy bien). El Q_to_v El método convierte el caudal en galones por minuto (GPM) a una velocidad de flujo en m/s. El friction_factor El método evalúa la ecuación de Haaland como se describió anteriormente, y el head_lossy head_loss_feet Evalúe la pérdida de cabeza de la tubería en metros y pies respectivamente (usando la ecuación de Darcy-Weisbach).
>>class Pipe():
>> def __init__(self, L, D, epsilon, fluid: Fluid):
>>
>> self.L = L
>> self.D = D
>> self.epsilon= epsilon
>>
>> self.fluid = fluid
>> @property
>> def A(self):
>> """computes cross-sectional area of pipe in m^2"""
>> return np.pi*(self.D/2)**2
>> def Q_to_v(self, gpm):
>> """Converts gpm to fluid speed in pipe in m/s"""
>> Q = gpm*GALLONS_TO_M3/60
>> return Q/self.A
>> def friction_factor(self, gpm):
>> """computes Darcy friction factor, given flow rate in gpm
>> This method uses Haaland's equation, wich is an explicit approximation
>> of the well-known, but implicit Colebrook equation
"""
>>
>> v = self.Q_to_v(gpm)
>>
>> Re = self.fluid.rho*v*self.D/self.fluid.mu
>>
>> e_over_d = self.epsilon/self.D
>>
>> f = (-1.8*np.log10((e_over_d/3.7)**1.11 + 6.9/Re))**-2
>> return f
>> def head_loss(self, gpm):
"""computes head loss in meters, given flow rate in gpm"""
>>
>> v = self.Q_to_v(gpm)
>>
>> f = self.friction_factor(gpm)
>>
>> hl = 0.5*f*(self.L/self.D)*v**2/self.fluid.g
>> return hl
>> def head_loss_feet(self, gpm):
"""computes head loss in feet, given flow rate in gpm"""
>> hl_meters = self.head_loss(gpm)
>> return hl_meters/FEET_TO_METERS
Veamos la clase de tubería en acción. Primero, podemos crear un Fluid objeto para agua y un Pipe Objeto que tiene 100 pies de largo y 1.25 pulgadas de diámetro.
>>water = Fluid()
>>pipe = Pipe(L=100*FEET_TO_METERS,
>> D=1.25*INCHES_TO_METERS,
>> epsilon=0.00006*INCHES_TO_METERS,
>> fluid=water)
A continuación, trazaremos la
curva de pérdida de cabeza en función del caudal en GPM. ¿No es
increíble lo limpio y fácil que se vuelve el siguiente código cuando
explotamos la programación orientada a objetos?
>>gpm_arr = np.linspace(1,30,100)
>>hl = [pipe.head_loss_feet(gpm) for gpm in gpm_arr]
>>fig, ax = plt.subplots(figsize=figsize)
>>ax.plot(gpm_arr, hl)
>>ax.set_xlabel('Flow Rate [gpm]')
>>ax.set_ylabel('Head Loss [ft]')
>>fig.tight_layout()
>>fig.savefig(plots_folder/'pipe_loss_curve.png')
Modelo de bomba
Tenemos
un modelo de trabajo de las pérdidas de cabeza de tubería, ahora
necesitamos un modelo de la cabeza producido por una bomba, HP (Q).
Estoy seguro de que hay modelos analíticos que se pueden usar para
determinar un comportamiento de las bombas, pero asumiremos que ya
tenemos una bomba específica, a saber, la siguiente bomba de caballos
fraccionaria aleatoria que encontré en línea :
La
mayoría de las bombas tendrán una hoja de datos que incluye una curva
de bomba correspondiente que caracteriza el comportamiento de las
bombas. Para nuestra bomba, obtenemos la curva de la bomba a
continuación (tenga en cuenta que por razones de derechos de autor, esta
es una recreación de la figura proporcionada por el fabricante: el
original se puede encontrar aquí ).
En
este punto, tenemos una imagen que representa el comportamiento de las
bombas, pero no un modelo matemático que realmente pueda usarse para
determinar cómo funcionará en el sistema. Este problema surge todo el
tiempo, y la forma en que lo sigo es 1) digitalizar los datos, y luego
2) envolver los datos discretos con un esquema de interpolación para
producir una función continua. Déjame ilustrar.
Paso 1) Hay muchas herramientas que existen para digitalizar los datos de imágenes: mi favorito personal es el Utilitiony WebLotDigitizer .
Se carga en la imagen de interés de la parcela, alinee los ejes y luego
extrae las curvas de datos deseadas (ya sea manualmente o con la
herramienta de extracción automática). Los datos se pueden exportar a un
archivo .csv.
Paso
2) Ahora que tenemos datos digitalizados, solo necesitamos envolverlo
con algún tipo de interpolador, esto es exactamente lo que el Pipe
La clase hace a continuación. El método de inicialización toma el
nombre del archivo .csv, almacena el nombre del archivo, carga los datos
en un Pandas DataFrame, almacenándolo en un data atribuir y luego pasar los datos al interp1d función de scipy. El interp1d
La función luego genera una nueva función que, por defecto, utiliza la
interpolación lineal para convertir los puntos de datos discretos en una
función continua (documentación completa para interp1d La función se puede encontrar aquí ). La función de interpolación recién generada se almacena en un _interp atributo para el acceso posterior. El Pipe La clase también contiene un bounds
método que devuelve una lista que contenga los valores min/max de la
velocidad de flujo en los datos de la curva de la bomba (esto se
utilizará en el algoritmo de búsqueda de raíz) y un head_gain_feet método que toma el caudal en GPM y llama a la función de interpolación subyacente que fue generada por interp1d.
>>class Pump():
>> def __init__(self, file):
>>
>> self.file = file
>>
>> self.data = pd.read_csv(file, names=['gpm', 'head [ft]']).set_index('gpm')
>>
>> self._interp = interp1d(self.data.index.to_numpy(), self.data['head [ft]'].to_numpy())
>> @property
>> def bounds(self):
>> """returns min and max flow rates in pump curve data"""
>> return [self.data.index.min(), self.data.index.max()]
>> def head_gain_feet(self, gpm):
>> """return head (in feet) produced by the pump at a given flow rate"""
>> return self._interp(gpm)
Podemos crear un Pump Objeto y mire los datos sin procesar en los que hemos leído.
>>pump = Pump('pump_data.csv')
>>pump.data.head()
También podemos trazar los datos
de la curva de la bomba con la curva de pérdida de tubería para ver
visualmente dónde funcionará el sistema.
>>head_loss = [pipe.head_loss_feet(gpm) for gpm in pump.data.index]
>>fig, ax = plt.subplots(figsize=figsize)
>>ax.plot(pump.data, label='Pump Curve')
>>ax.plot(pump.data.index, head_loss, label='Pipe Head Loss')
>>ax.set_xlabel("Flow Rate [gpm]")
>>ax.set_ylabel("Head [ft]")
>>ax.legend(frameon=True, facecolor='w', framealpha=1, loc=6)
>>fig.tight_layout()
>>fig.savefig(plots_folder/'pump_curve_with_losses.png')
Modelo de sistema
Finalmente
tenemos la infraestructura en su lugar para resolver el punto de
funcionamiento del sistema de bomba/tubería. El último paso es crear un System clase que toma un Pipey Pump objeto y realiza la operación de resolución de raíces. Como se puede ver en el código a continuación, el System La clase toma y almacena un Pipey Pump objeto. Luego usa los dos objetos para crear un residual Método que calcula la diferencia entre el cabezal de la bomba y la pérdida de la cabeza de la tubería. Este residual el método se usa luego en el get_operating_point Método para resolver realmente el punto de funcionamiento del sistema. El método envuelve el root_scalar Funciona desde SciPy, que actúa como una interfaz para varios algoritmos de resolución de raíces. Dejaremos el root_scalar
Función Elija el algoritmo que considere mejor, pero para ayudarlo,
especificaremos un intervalo de soporte en el que sabemos que la raíz se
encuentra entre. En nuestro caso, este intervalo de corchete son los
límites de flujo superior e inferior de los datos de la curva de la
bomba. Documentación completa en el root_scalar La función se puede encontrar aquí .
Consejo profesional: el proceso de inyección del Pipe y Pump objetos en el System clase (en lugar de que la clase del sistema cree un Pipe y Pump
El objeto en la instancia) se llama "inyección de dependencia". Esto
generalmente se considera una buena práctica de codificación, ya que
hace que el código sea más modular, extensible y más fácil de
depurar/probar.
>>class System():
>> def __init__(self, pipe: Pipe, pump: Pump):
>> self.pipe = pipe
>> self.pump = pump
>> def residual(self, gpm):
"""
>> Computes the difference between the head produced by the pump
>> and the head loss in the pipe. At steady state, the pump head and
>> head loss will be equal and thus the residual function will go to zero
"""
>> return self.pump.head_gain_feet(gpm) - self.pipe.head_loss_feet(gpm)
>> def get_operating_point(self):
>> """solve for the flow rate where the residual function equals zero.
>> i.e. the pump head equals the pipe head loss"""
>> return root_scalar(self.residual, bracket=self.pump.bounds)
Creemos un System y corre el get_operating_point Método para observar los frutos de nuestro trabajo. Como puede ver la salida del código, el get_operating point El método simplemente desmonta la salida del root_scalar función, que es un RootResults objeto. Este objeto es esencialmente solo un contenedor que almacena varios atributos, el más importante es el root atributo tal como contiene la solución a nuestro problema.
>>sys = System(pipe, pump)
>>res = sys.get_operating_point()
>>res
Podemos trazar la misma bomba y
curvas de pérdida de cabeza nuevamente, esta vez agregando una línea
vertical en el punto de operación calculado de estado estacionario.
>>head_loss = [pipe.head_loss_feet(gpm) for gpm in pump.data.index]
>>fig, ax = plt.subplots(figsize=figsize)
>>ax.plot(pump.data, label='Pump Curve')
>>ax.plot(pump.data.index, head_loss, label='Pipe Head Loss')
>>ax.axvline(res.root, color='k', ls='--', lw=1)
>>ax.legend(frameon=True, facecolor='white', framealpha=1, loc=6)
>>ax.set_xlabel("Flow Rate [gpm]")
>>ax.set_ylabel("Head [ft]")
>>ax.set_title(f'Operating Point = {res.root:.1f} gpm')
>>fig.tight_layout()
>>fig.savefig(plots_folder/'intersection_solution.png')
¡Voila!
Hemos determinado programáticamente el punto de operación de nuestro
sistema. Y debido a que lo hemos hecho usando un marco de codificación
algo genérico, ¡podemos probar fácilmente el mismo análisis usando
diferentes bombas o tuberías! Incluso podríamos extender nuestro código
para incluir múltiples bombas o varios accesorios de tubería/ramas de
tubería.
Exploración de diseño
Como
un pequeño ejemplo final, destacando los beneficios de configurar el
código como lo hicimos, realizaremos una exploración de diseño. Usando
la misma bomba, nos gustaría comprender los impactos que tiene la
longitud de la tubería en el caudal volumétrico del sistema. Para hacer
esto, simplemente recorre una matriz de longitudes de tubería (que van
desde 100 a 1000 pies), actualicamos el atributo de longitud del Pipe objeto almacenado en el System,
y luego vuelva a competir el punto de funcionamiento del sistema,
agregando el resultado a una lista. Finalmente, trazamos la tasa de
flujo de agua en función de la longitud de la tubería.
>>lengths_feet = np.linspace(100, 1000, 1000)
>>lengths_meters = lengths_feet*FEET_TO_METERS
>>flow_rates = []
>>for l in lengths_meters:
>>
>> sys.pipe.L = l
>>
>> res = sys.get_operating_point()
>>
>> flow_rates.append(res.root)
>>fig, ax = plt.subplots(figsize=figsize)
>>ax.plot(lengths_feet, flow_rates)
>>ax.set_xlabel("Pipe Length [ft]")
>>ax.set_ylabel("Flow Rate [gpm]")
>>fig.tight_layout()
>>fig.savefig(plots_folder/'flow_vs_pipe_length.png')
En
unas pocas líneas de código, pudimos extraer una visión significativa
del comportamiento de nuestro sistema. Si este fuera un problema de
diseño, estas ideas podrían impulsar las opciones de diseño clave.
Conclusión
Este
artículo, aunque se centró fuertemente en un problema de ejemplo
específico de dominio, destaca algunos aspectos de un flujo de trabajo
que termino usando mucho. El problema del análisis de puntos operativos
surge constantemente en ingeniería y ciencia, y aunque hay muchas
maneras de abordar el problema, algunos métodos son más robustos,
extensibles y flexibles que otros. La metodología (principios de
formulación de problemas y estructuración de código) utilizadas en este
artículo me ha servido increíblemente bien, ¡y espero que otros se
inspiraran para adoptar un flujo de trabajo similar!
Siéntase
libre de dejar cualquier comentario o pregunta que pueda tener o
conectarse conmigo en LinkedIn: estaría más que feliz de aclarar
cualquier punto de incertidumbre. Finalmente, le animo a que juegue con
el código usted mismo (o incluso lo use como plantilla inicial para sus
propios flujos de trabajo): el cuaderno Jupyter para este artículo se
puede encontrar en mi Github .
Nicholas Hemenway
