Creación de un marco genérico para el análisis del punto de operación

Ciertas clases de problemas surgen con frecuencia en ingeniería. El enfoque de este artículo está en un tipo específico de problema que surge con tanta frecuencia en mi trabajo diario que pensé en compartir cómo lo soluciono usando Python. ¿De qué tipo de problema estamos hablando? ¡El problema de resolver el punto de funcionamiento de un sistema! Ilustremos lo que quiero decir con un ejemplo simple antes de sumergirnos en algo un poco más complejo con el código.

Nos gustaría resolver el punto de operación del circuito simple que se muestra a continuación. Esto se puede hacer reorganizando la ley de Ohm (V=IR) para aislar la corriente en términos del voltaje de entrada y la resistencia conocidos.

Sencillo, ¿verdad? Desafortunadamente, la mayoría de los problemas del mundo real nunca son tan fáciles. Por ejemplo, ¿qué pasaría si le dijera que a medida que la resistencia se calienta, su valor de resistencia cambia, esencialmente haciendo que la resistencia sea una función de la corriente? Terminamos con una ecuación de la siguiente forma:

Sin conocer la forma funcional real de la resistencia, no podemos resolver la corriente aislándola algebraicamente. Además, ¿qué pasa si la ecuación es complicada y no es posible aislar la corriente por sí misma? O, tal vez, la resistencia se da en términos de corriente como datos discretos tabulados, entonces ni siquiera tendríamos una expresión algebraica para manipular para tratar de resolver la corriente. Entonces, ¿cómo haríamos para determinar la corriente en el circuito? Necesitamos un enfoque más general para resolver este problema.

La solución general a un problema como este es plantearlo como un problema de búsqueda de raíces. En realidad, esto es increíblemente fácil de hacer: literalmente, solo tenemos que restar el lado derecho de la ecuación del lado izquierdo, de modo que obtengamos una ecuación que sea igual a cero. Hacerlo produce lo siguiente:

Al hacer esto, hemos vuelto a plantear nuestro problema. En lugar de resolver la corriente directamente en términos de todas las demás variables, podemos intentar encontrar el valor de la corriente que se puede ingresar en el lado izquierdo de la ecuación para que se evalúe como cero. ¿Por qué formulamos el problema de esta manera? ¡Porque hay un montón de algoritmos numéricos que existen (método de bisección, método de Newton, etc.) para resolver este tipo exacto de problema! Y a la mayoría de los algoritmos no les importa cuán complicado sea el lado izquierdo de la ecuación; ni siquiera tiene que tener una forma algebraica cerrada (es decir, podría estar compuesto de datos discretos interpolados, integrales evaluadas numéricamente o, literalmente, cualquier tipo de función de complejidad arbitraria a evaluar). Siempre que podamos plantear nuestro problema en la forma de f(x)=0,

El resto de este artículo le mostrará un ejemplo de cómo aplicar la metodología de búsqueda de raíces a un problema del mundo real un poco más complicado, con énfasis en técnicas sólidas de estructuración y organización de código en Python. Aunque el problema (determinar el caudal de agua en un sistema de tubería/bomba) es algo específico de un dominio, la metodología y las técnicas de codificación que se utilizan son completamente generales y aplicables a todos los dominios de la ingeniería. Con esto en mente, trataré de mantener los aspectos de modelado físico del problema a un alto nivel, de modo que, independientemente de los antecedentes técnicos, los objetivos de aprendizaje principales del artículo aún se vean claramente.

Como nota al margen, mi "especialidad" de dominio en estos días se encuentra en el ámbito de los controles de motores y la electrónica de potencia, y estoy muy alejado de las aplicaciones de bombeo/tuberías. No he tocado el tema en años, pero pensé que sería un ejemplo interesante del tema en cuestión. Estoy seguro de que hay muchas personas que están mucho más calificadas que yo para hablar sobre los detalles específicos del modelado de bombas/tuberías, pero mi intención con este artículo es la metodología, no cómo resolver problemas de tuberías/bombas. De todos modos, ¡acepto abiertamente los comentarios o sugerencias de mejora de aquellos que tienen más conocimientos en el campo!

El problema

Nos gustaría transferir agua de un tanque a otro. Ya tenemos una bomba y algunas tuberías que se pueden usar para conectar los dos tanques y queremos obtener una estimación de cuánto tiempo llevará transferir toda el agua. Se conoce el volumen de cada tanque, por lo que si podemos estimar la tasa de flujo de agua entre los tanques, podemos estimar cuánto tiempo tomará el proceso de transferencia. El aparato completo se muestra a continuación.

Este problema específico (que se puede clasificar como un problema de "flujo interno") se entiende muy bien dentro del campo de la ingeniería mecánica. Sin embargo, para aquellos menos familiarizados, o que necesitan una revisión rápida, la forma en que generalmente resolvemos estos problemas es con la ecuación de Bernoulli (que se muestra a continuación).

La ecuación de Bernoulli es esencialmente una declaración de conservación de energía que nos dice cómo la energía de una partícula de fluido se transforma entre diferentes mecanismos de energía a medida que el fluido atraviesa una línea de corriente (la ruta de flujo que seguiría una partícula imaginaria si se dejara caer en el fluido). El lado izquierdo de la ecuación representa la energía total por peso de una partícula de fluido en cualquier primera ubicación arbitraria (ubicación 1) dentro del fluido, y es la suma de un término de potencial gravitacional, un término cinético y un término de presión. A medida que el fluido atraviesa el sistema, la energía debe conservarse y, por lo tanto, la energía total en cualquier segundo punto arbitrario (ubicación 2) a lo largo de la línea de corriente (representada por el lado derecho de la ecuación) debe ser igual a la energía total en la ubicación 1 .

La forma anterior de la ecuación de Bernoulli se conoce como la forma de "cabeza" de la ecuación porque cada término tiene unidades de longitud/altura. Esto es conveniente para nuestra intuición porque esencialmente estamos equiparando la energía de cada término con la energía potencial gravitatoria equivalente de una columna de fluido con una altura de la cabeza dada. Sin embargo, una limitación importante de la ecuación de Bernoulli es que asume que no hay pérdidas en el sistema (lo cual no es una gran suposición). Para superar esta limitación, podemos complementar la ecuación con dos términos adicionales de la siguiente manera:

Los términos Hp(Q) y Hl(Q) representan la carga añadida al sistema por una bomba y la carga perdida en el sistema por efectos del mundo real (como fricción, viscosidad, etc.) respectivamente. Tenga en cuenta que ambos términos son funciones de la velocidad de flujo del fluido del sistema, Q. (Como una consecuencia interesante del párrafo anterior que describe la interpretación de la cabeza, la cabeza de la bomba le dice qué tan alto teóricamente podría empujar una bomba un fluido). Examinaremos los términos de bombeo y pérdida más a fondo en un momento, pero antes de hacerlo, 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 que estén 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 supondremos que la altura del agua en los dos tanques es la misma en el momento en que analizamos el sistema, de modo que Z1=Z2. Después de simplificar el álgebra, vemos que casi todos los términos se cancelan y nos quedamos con el hecho de que la carga producida por la bomba debe ser igual a la carga 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 siguiente figura. La cabeza producida por una bomba disminuye al aumentar el caudal, mientras que las pérdidas en un sistema de tuberías aumentan al aumentar el caudal. El punto donde se cruzan las dos curvas (cabeza de la bomba = pérdida de carga) determina el punto de funcionamiento (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íz que estamos buscando. Es decir, hemos planteado nuestro problema de la siguiente manera: encuentre el caudal (Q) tal que el lado izquierdo de la siguiente ecuación sea igual a cero. En este punto, la carga de la bomba será igual a las pérdidas de carga 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 una experiencia razonable en Python). En su lugar, voy a centrar mis esfuerzos en garantizar que la narrativa y la estructuración del código sean claras, 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. Será evidente cómo se usa cada uno de los módulos 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 para el problema en cuestión. El bloque de código también establece algunas configuraciones de trazado predeterminadas (a gusto personal), crea una carpeta para guardar las cifras generadas y define algunas constantes de conversión de unidades que nos hacen la vida más fácil 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

#these are the key libraries for solving the problem
from scipy.interpolate import interp1d
from scipy.optimize import root_scalar

#set plotting defaults
plt.style.use('seaborn-v0_8-darkgrid')
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 12
figsize = (6.4,4)

#make folder to save plots to
plots_folder = Path('plots')
plots_folder.mkdir(exist_ok=True)

#define conversion constants for ease of use later
INCHES_TO_METERS = 25.4/1000
FEET_TO_METERS = 12*INCHES_TO_METERS
GALLONS_TO_M3 = 0.0037854118 #convert gallons to m^3

@dataclass
class Fluid():
    #fluid defaults to water properties
    rho: float = 997 #kg/m^3
    mu: float = 0.0007972 #N-s/m^2 = kg/m-s
    g: float = 9.81 #m/s^2

El siguiente paso es modelar las pérdidas de carga de la tubería (el término Hl(Q) en la ecuación de Bernoulli extendida anterior). Esto generalmente se hace usando 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 del 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, sino que 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, épsilon es la rugosidad de la superficie de la tubería (que se puede encontrar en las tablas de los libros de texto de ingeniería) y Re es el famoso número de Reynolds, calculado como se muestra a continuación.

Por último, podemos notar que el volumen barrido por unidad de tiempo, o tasa de flujo volumétrico (Q) es igual al área de la sección transversal (A) de la tubería por la velocidad del flujo (v). Por lo tanto, dada una tasa de flujo en la tubería, podemos calcular la velocidad de flujo correspondiente en la tubería como:

Con suerte, todas estas ecuaciones no restan valor al panorama general: solo estamos viendo un modelo particular de cálculo de la pérdida de carga en una tubería. Dada una tasa de flujo y las dimensiones de la tubería, primero calcule la velocidad de flujo correspondiente, luego complete las ecuaciones anteriores para calcular la pérdida de carga de la tubería. Esto es exactamente lo que Pipeimplementa la clase (que se muestra a continuación).

El método de inicialización almacena las dimensiones de las tuberías (se supone que todas están en metros) y las propiedades del fluido. El Amétodo calcula el área de la sección transversal de la tubería (para aquellos que no están familiarizados con el @propertydecorador, este artículo lo explica muy bien). El Q_to_vmétodo convierte la tasa de flujo en galones por minuto (gpm) a una velocidad de flujo en m/s. El friction_factormétodo evalúa la ecuación de Haaland como se describe arriba, head_lossy head_loss_feetevalúa la pérdida de carga 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):
        
        #pipe dimensions are all assumed to be in meters
        self.L = L
        self.D = D
        self.epsilon= epsilon
        
        #fluid properties
        self.fluid = fluid
        
    @property
    def A(self):
        """computes cross-sectional area of pipe in m^2"""
        
        return np.pi*(self.D/2)**2 #area in m^2
    
    def Q_to_v(self, gpm):
        """Converts gpm to fluid speed in pipe in m/s"""
        
        Q = gpm*GALLONS_TO_M3/60 #flow rate in m^3/s
        return  Q/self.A #flow velocity in m/s
    
    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
        """
        #first get flow velocity from flow rate and pipe dimensions
        v = self.Q_to_v(gpm)
        #compute Reynold's number
        Re = self.fluid.rho*v*self.D/self.fluid.mu
        #compute relative roughness
        e_over_d = self.epsilon/self.D
        #use Haaland's equation
        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"""
        
        #get flow velocity
        v = self.Q_to_v(gpm)
        #get Darcy friction factor
        f = self.friction_factor(gpm)
        #compute head loss in meters
        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

#create fluid object for water
water = Fluid()

#create pipe segment with water flowing in it
pipe = Pipe(L=100*FEET_TO_METERS, 
            D=1.25*INCHES_TO_METERS, 
            epsilon=0.00006*INCHES_TO_METERS,
            fluid=water)

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')

      
                

                    
                    
                
            

Tenemos un modelo funcional de las pérdidas de carga de la tubería; ahora necesitamos un modelo de la carga producida por una bomba, Hp(Q). Estoy seguro de que hay modelos analíticos que se pueden usar para determinar el comportamiento de una bomba, pero supondremos que ya tenemos una bomba específica, es decir, la siguiente bomba de potencia 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 la bomba anterior, obtenemos la siguiente curva de bomba:

En este punto, tenemos una imagen que representa el comportamiento de las bombas, pero no un modelo matemático que pueda usarse para determinar cómo funcionará en el sistema. Este problema surge todo el tiempo, y la forma en que lo soluciono 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) Existen muchas herramientas para digitalizar datos de imágenes: mi favorita es la utilidad en línea gratuita WebPlotDigitizer . Cargue la imagen de la trama de interés, alinee los ejes y luego extraiga las curvas de datos deseadas (ya sea manualmente o con la herramienta de extracción automática). Luego, los datos se pueden exportar a un archivo .csv.

Paso 2) Ahora que tenemos datos digitalizados, solo necesitamos envolverlos con algún tipo de interpolador; esto es exactamente lo que Pipehace la clase 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 DataFrame de pandas, los almacena en un dataatributo y luego pasa los datos a la interp1dfunción desde scipy. Luego, la interp1dfunción genera una nueva función que, de forma predeterminada, utiliza la interpolación lineal para convertir los puntos de datos discretos en una función continua (la documentación completa de la interp1dfunción se puede encontrar aquí ). La función de interpolación recién generada se almacena luego en un _interpatributo para un acceso posterior. La Pipeclase también contiene unboundsmétodo que devuelve una lista que contiene los valores mínimos/máximos de la tasa de flujo en los datos de la curva de la bomba (esto se usará en el algoritmo de búsqueda de raíz) y un head_gain_feetmétodo que toma la tasa de flujo en gpm y llama a la función de interpolación subyacente que fue generado por interp1d.

class Pump():
    
    def __init__(self, file):    
        #store file name
        self.file = file 
        #read data into pandas dataframe and assign column names
        self.data = pd.read_csv(file, names=['gpm', 'head [ft]']).set_index('gpm')
        #create continuous interpolation function
        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)

pump = Pump('pump_data.csv')
pump.data.head()

      
                

                    
                    
                
            

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')

      
                

                    
                    
                
            

Finalmente contamos con la infraestructura para solucionar el punto de operación del sistema bomba/tubería. El último paso es crear una Systemclase que tome un objeto Pipey Pumprealice la operación de resolución raíz. Como se puede ver en el código a continuación, la Systemclase toma y almacena un objeto Pipey Pump. A continuación, utiliza los dos objetos para crear un residualmétodo que calcula la diferencia entre la pérdida de carga de la tubería y la altura de la bomba. Este residualmétodo se usa luego en el get_operating_pointmétodo para resolver realmente el punto de operación del sistema. El método envuelve la root_scalarfunción de scipy, que actúa como una interfaz para varios algoritmos de resolución de raíz. Dejaremos que elroot_scalarelija el algoritmo que mejor se ajuste, pero para ayudarlo, especificaremos un intervalo de horquillado entre el que sabemos que se encuentra la raíz. En nuestro caso, este intervalo de horquillado son los límites superior e inferior del caudal de los datos de la curva de la bomba. La documentación completa sobre la root_scalarfunción se puede encontrar aquí .

Consejo profesional: el proceso de inyectar los objetos Pipey en la clase (en lugar de que la clase del sistema cree un objeto y en la creación de instancias) 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.PumpSystemPipePump

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(sys.residual, bracket=pump.bounds)

sys = System(pipe, pump)

res = sys.get_operating_point()
res

      
                

                    
                    
                
            

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')
#plot vertical line at operating point
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')

      
                

                    
                    
                
            

Exploración de diseño

Como pequeño ejemplo final, destacando los beneficios de configurar el código de la forma en que lo hicimos, realizaremos una exploración del diseño. Usando la misma bomba, nos gustaría comprender los impactos que tiene la longitud de la tubería en la tasa de flujo volumétrico en el sistema. Para hacer esto, simplemente recorremos una serie de longitudes de tubería (que van desde 100 a 1000 pies), actualizamos el atributo de longitud del Pipeobjeto almacenado en el Systemarchivo y luego volvemos a calcular el punto operativo del sistema, agregando el resultado a un lista. Finalmente, representamos gráficamente el caudal de agua en función de la longitud de la tubería.

#sweep pipe length from 100 to 1000 feet
lengths_feet = np.linspace(100, 1000, 1000)
lengths_meters = lengths_feet*FEET_TO_METERS

flow_rates = []
for l in lengths_meters:
    #update pipe length
    sys.pipe.L = l 
    #compute new flow rate solution
    res = sys.get_operating_point()
    #append solution to flow rates list
    flow_rates.append(res.root) 
    
#plot results
fig, ax = plt.subplots(figsize=figsize)
ax.plot(lengths_feet, flow_rates)
ax.set_xlabel("Pipe Length [ft]")
ax.set_ylabel("Flow Rate [gpm]")
# ax.set_ylim(bottom=0)
fig.tight_layout()
fig.savefig(plots_folder/'flow_vs_pipe_length.png')

      
                

                    
                    
                
            

Conclusión

Este artículo, aunque se centró fuertemente en un problema de ejemplo específico de un dominio, destaca algunos aspectos de un flujo de trabajo que termino usando mucho. El problema del análisis del punto de operación surge constantemente en la ingeniería y la ciencia, y aunque hay muchas formas de abordar el problema, algunos métodos son más robustos, extensibles y flexibles que otros. La metodología (formulación de problemas y principios de estructuración de código) utilizada en este artículo me ha servido increíblemente bien, ¡y espero que otros se sientan inspirados 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é más que feliz de aclarar cualquier punto de incertidumbre. Finalmente, lo animo a que juegue con el código usted mismo (o incluso lo use como una plantilla de inicio para sus propios flujos de trabajo): el Jupyter Notebook para este artículo se puede encontrar en mi Github .

Nicolás Hemenway

• Si disfrutaste esto, sígueme en Medium
• Considere suscribirse a las actualizaciones por correo electrónico
• ¿Interesado en colaborar? Conectémonos en LinkedIn