Usando Python para resolver um dos problemas mais comuns em engenharia

Jan 04 2023

Criando uma Estrutura Genérica para Análise de Ponto de Operação Certas classes de problemas surgem com frequência na engenharia. O foco deste artigo é um tipo específico de problema que surge com tanta frequência em meu trabalho diário que decidi compartilhar como resolvo isso usando Python.

Criando uma Estrutura Genérica para Análise de Ponto Operacional

Foto da SpaceX no Unsplash

Certas classes de problemas surgem com frequência na engenharia. O foco deste artigo é um tipo específico de problema que surge com tanta frequência em meu trabalho diário que decidi compartilhar como resolvo isso usando Python. De que tipo de problema estamos falando? O problema de resolver o ponto de operação de um sistema! Vamos ilustrar o que quero dizer com um exemplo simples antes de mergulhar em algo um pouco mais complexo com código.

Gostaríamos de resolver o ponto de operação do circuito simples mostrado abaixo. Isso pode ser feito rearranjando a lei de Ohm (V=IR) para isolar a corrente em termos da tensão de entrada e resistência conhecidas.

Simples, certo? Infelizmente, a maioria dos problemas do mundo real nunca é tão fácil. Por exemplo, e se eu dissesse que conforme o resistor esquenta, seu valor de resistência muda, essencialmente tornando a resistência uma função da corrente. Chegamos a uma equação da seguinte forma:

Sem conhecer a forma funcional real da resistência, não podemos simplesmente calcular a corrente isolando-a algebricamente. Além disso, e se a equação for complicada e não for possível isolar a corrente sozinha? Ou talvez a resistência seja dada em termos de corrente como dados discretos tabulados - então não teríamos nem mesmo uma expressão algébrica para manipular para tentar resolver a corrente. Como faríamos para determinar a corrente no circuito então? Precisamos de uma abordagem mais geral para resolver esse problema.

A solução geral para um problema como esse é colocá-lo como um problema de descoberta de raiz. Na verdade, isso é incrivelmente fácil de fazer - literalmente, só temos que subtrair o lado direito da equação do lado esquerdo, de modo que obtenhamos uma equação igual a zero. Fazer isso produz o seguinte:

Ao fazer isso, redefinimos nosso problema. Em vez de resolver a corrente diretamente em termos de todas as outras variáveis, podemos tentar encontrar o valor da corrente que pode ser inserido no lado esquerdo da equação para torná-la avaliada como zero. Por que formulamos o problema dessa maneira? Porque existe uma tonelada de algoritmos numéricos (método da bisseção, método de Newton, etc.) para resolver esse tipo exato de problema! E a maioria dos algoritmos não se importa com o quão complicado é o lado esquerdo da equação - ele nem precisa ter uma forma algébrica fechada (ou seja, pode ser composto de dados discretos interpolados, integrais avaliados numericamente ou literalmente qualquer tipo de função de complexidade arbitrária a avaliar). Desde que possamos colocar nosso problema na forma de f(x)=0,

O restante deste artigo mostrará um exemplo de como aplicar a metodologia de localização de raiz a um problema do mundo real um pouco mais complicado, com ênfase na estruturação de código de som e técnicas de organização em Python. Embora o problema (determinar a taxa de fluxo de água em um sistema de tubulação/bomba) seja um pouco específico do domínio, a metodologia e as técnicas de codificação usadas são completamente gerais e aplicáveis a todos os domínios da engenharia. Com isso em mente, tentarei manter os aspectos de modelagem física do problema em um nível elevado, de modo que, independentemente de sua formação técnica, os objetivos primários de aprendizado do artigo ainda apareçam claramente.

Como uma observação lateral, minha “especialidade” de domínio atualmente está no domínio dos controles de motores e eletrônica de potência, e estou muito distante das aplicações de bombeamento/tubulação. Não toco no assunto há anos, mas pensei que seria um exemplo interessante do tópico em questão. Tenho certeza de que existem muitas pessoas por aí que são muito mais qualificadas para falar sobre as especificidades da modelagem de bombas/tubos do que eu, mas minha intenção com este artigo é a metodologia - não como resolver problemas de tubos/bombas. Independentemente disso, aceito abertamente comentários ou sugestões de melhoria daqueles que são mais conhecedores do campo!

O problema

Gostaríamos de transferir água de um tanque para outro. Já temos uma bomba e algumas tubulações que podem ser usadas para conectar os dois tanques e queremos ter uma estimativa de quanto tempo levará para transferir toda a água. O volume de cada tanque é conhecido, portanto, se pudermos estimar a vazão de água entre os tanques, podemos estimar quanto tempo levará o processo de transferência. O aparato completo é mostrado abaixo.

Este problema específico (que pode ser classificado como um problema de “fluxo interno”) é muito bem compreendido no campo da engenharia mecânica. No entanto, para aqueles menos familiarizados ou que precisam de uma revisão rápida, a maneira como normalmente resolvemos esses problemas é com a equação de Bernoulli (mostrada abaixo).

A equação de Bernoulli é essencialmente uma declaração de conservação de energia que nos diz como a energia de uma partícula de fluido é transformada entre diferentes mecanismos de energia conforme o fluido atravessa uma linha de corrente (o caminho do fluxo que uma partícula imaginária seguiria se caísse no fluido). O lado esquerdo da equação representa a energia total por peso de uma partícula de fluido em qualquer primeira localização arbitrária (localização 1) dentro do fluido e é a soma de um termo de potencial gravitacional, termo cinético e termo de pressão. À medida que o fluido atravessa o sistema, a energia deve ser conservada e, portanto, a energia total em qualquer segundo ponto arbitrário (localização 2) ao longo da linha de corrente (representada pelo lado direito da equação) deve ser igual à energia total na localização 1 .

A forma acima da equação de Bernoulli é conhecida como a forma “cabeça” da equação porque cada termo tem unidades de comprimento/altura. Isso é conveniente para a nossa intuição porque estamos essencialmente igualando a energia de cada termo à energia potencial gravitacional equivalente de uma coluna de fluido com a altura da carga dada. Uma grande limitação da equação de Bernoulli, no entanto, é que ela assume que não há perdas no sistema (o que não é uma boa suposição). Para superar essa limitação, podemos complementar a equação com dois termos adicionais, como segue:

Os termos Hp(Q) e Hl(Q) representam a carga adicionada ao sistema por uma bomba e a carga perdida no sistema devido a efeitos do mundo real (como fricção, viscosidade etc.), respectivamente. Observe que ambos os termos são funções da taxa de fluxo de fluido do sistema, Q. (Como uma consequência interessante do parágrafo acima, que descreve a interpretação da altura manométrica, a altura manométrica da bomba informa a que altura uma bomba poderia, teoricamente, empurrar um fluido). Examinaremos os termos de bomba e perda mais detalhadamente daqui a pouco, mas antes disso, vamos simplificar a equação acima para nosso problema específico.

Olhando novamente para o sistema acima, escolheremos convenientemente nossas duas localizações para a equação de Bernoulli, de modo que a maioria dos termos se cancelem. Podemos fazer isso escolhendo os locais 1 e 2 para estar na superfície livre da água de cada tanque respectivamente, onde a pressão é constante e igual à pressão atmosférica (P1=P2), e a velocidade é aproximadamente constante e zero (V1 =V2=0). Também assumiremos que a altura da água nos dois tanques é a mesma no instante em que estamos analisando o sistema de forma que Z1=Z2. Depois de simplificar a álgebra, vemos que quase todos os termos se cancelam e ficamos com o fato de que a altura manométrica produzida pela bomba deve ser igual à carga perdida no sistema devido a não idealidades. Dito de outra forma, a bomba está compensando quaisquer perdas de energia no sistema.

Esta situação pode ser vista qualitativamente na figura abaixo. A carga produzida por uma bomba diminui com o aumento da vazão, enquanto as perdas em um sistema de tubulação aumentam com o aumento da vazão. O ponto de interseção das duas curvas (caudal da bomba = perda de carga) determina o ponto de operação (vazão) do sistema.

A última etapa antes de entrarmos no código é colocar o problema como um problema de localização de raiz. Subtraindo o lado direito da equação do lado esquerdo, obtemos a raiz do problema que procuramos. Ou seja, colocamos nosso problema da seguinte forma: encontre a taxa de fluxo (Q) tal que o lado esquerdo da equação abaixo seja igual a zero. Neste ponto, a carga da bomba será igual às perdas de carga do sistema.

O código

Para evitar perder a visão geral do que estamos fazendo, não vou explicar cada pequeno detalhe do código (presumo que você já tenha um conhecimento razoável em Python). Em vez disso, vou concentrar meus esforços em garantir que a narrativa e a estruturação do código sejam claras e entrar em mais detalhes conforme necessário. Como sempre, sinta-se à vontade para fazer qualquer pergunta se algo não estiver claro.

Configurar

Começaremos importando todos os módulos necessários. Ficará aparente como cada um dos módulos é usado mais tarde, mas vale a pena notar que as principais declarações de importação são as do scipy. Estas são as funções que são específicas para o problema em questão. O bloco de código também define algumas configurações de plotagem padrão (para gosto pessoal), cria uma pasta para salvar as figuras geradas e define algumas constantes de conversão de unidades que facilitam nossa vida posteriormente no 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

O próximo passo é modelar as perdas de carga do tubo (o termo Hl(Q) na equação estendida de Bernoulli acima). Isso geralmente é feito usando a equação de Darcy-Weisbach mostrada abaixo, onde f é um fator de atrito (mais sobre isso em breve), v é a velocidade do fluxo, g é a gravidade e L e D são o comprimento e o diâmetro do tubo, respectivamente.

Infelizmente, o fator de atrito (f) não é constante, mas também depende da velocidade do fluxo, das propriedades do fluido e das dimensões do tubo. Existem vários modelos para calcular f, mas usaremos a equação de Haaland, mostrada abaixo.

Nesta equação, epsilon é a rugosidade da superfície do tubo (que pode ser encontrada nas tabelas dos livros de engenharia) e Re é o famoso número de Reynolds, calculado conforme abaixo.

Por fim, podemos observar que o volume varrido por unidade de tempo, ou vazão volumétrica (Q) é igual à área da seção transversal (A) do tubo vezes a velocidade do fluxo (v). Assim, dada uma taxa de fluxo no tubo, podemos calcular a velocidade de fluxo correspondente no tubo como:

Esperançosamente, todas essas equações não prejudicam o quadro geral - estamos apenas analisando um modelo específico de cálculo da perda de carga em um tubo. Dada uma taxa de fluxo e as dimensões do tubo, primeiro calcule a velocidade do fluxo correspondente e, em seguida, conecte as equações acima para calcular a perda de carga do tubo. Isso é exatamente o que a Pipeclasse (mostrada abaixo) implementa.

O método de inicialização armazena as dimensões dos tubos (todas consideradas em metros) e as propriedades do fluido. O Amétodo calcula a área da seção transversal do tubo (para aqueles que não estão familiarizados com o @propertydecorador, este artigo explica muito bem). O Q_to_vmétodo converte a taxa de fluxo em galões por minuto (gpm) em uma velocidade de fluxo em m/s. O friction_factormétodo avalia a equação de Haaland conforme descrito acima, head_losse head_loss_feetavalia a perda de carga do tubo em metros e pés, respectivamente (usando a equação 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')

Temos um modelo funcional das perdas de carga do tubo — agora precisamos de um modelo da carga produzida por uma bomba, Hp(Q). Tenho certeza de que existem modelos analíticos que podem ser usados para determinar o comportamento de uma bomba, mas vamos assumir que já temos uma bomba específica - ou seja, a seguinte bomba de potência fracionada aleatória que encontrei online :

A maioria das bombas terá uma folha de dados que inclui uma curva de bomba correspondente que caracteriza o comportamento das bombas. Para a bomba acima, obtemos a seguinte curva da bomba:

Neste ponto, temos uma imagem que descreve o comportamento das bombas, mas não um modelo matemático que possa realmente ser usado para determinar como ela funcionará no sistema. Esse problema surge o tempo todo, e a maneira como contorná-lo é 1) digitalizar os dados e, em seguida, 2) agrupar os dados discretos com um esquema de interpolação para produzir uma função contínua. Deixe-me ilustrar.

Etapa 1) Existem muitas ferramentas para digitalizar dados de imagem - minha favorita é o utilitário on-line gratuito WebPlotDigitizer . Você carrega a imagem do gráfico de interesse, alinha os eixos e extrai as curvas de dados desejadas (manualmente ou com a ferramenta de extração automática). Os dados podem então ser exportados para um arquivo .csv.

Etapa 2) Agora que temos os dados digitalizados, só precisamos envolvê-los com algum tipo de interpolador — é exatamente isso que a Pipeclasse faz abaixo. O método de inicialização pega o nome do arquivo .csv, armazena o nome do arquivo, carrega os dados em um Pandas DataFrame, armazenando-o em um dataatributo e, em seguida, passa os dados para a interp1dfunção do scipy. A interp1dfunção então gera uma nova função que, por padrão, usa interpolação linear para transformar os pontos de dados discretos em uma função contínua (a documentação completa da interp1dfunção pode ser encontrada aqui ). A função de interpolação recém-gerada é armazenada em um _interpatributo para acesso posterior. A Pipeclasse também contém umboundsmétodo que retorna uma lista contendo os valores mínimo/máximo da taxa de fluxo nos dados da curva da bomba (isso será usado no algoritmo de localização da raiz) e um head_gain_feetmétodo que recebe a taxa de fluxo em gpm e chama a função de interpolação subjacente que foi gerado 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 temos a infraestrutura instalada para resolver o ponto de operação do sistema bomba/tubo. A última etapa é criar uma Systemclasse que receba um objeto Pipeand Pumpe execute a operação de solução de raiz. Como pode ser visto no código abaixo, a Systemclasse recebe e armazena um objeto Pipeand . PumpEm seguida, ele usa os dois objetos para criar um residualmétodo que calcula a diferença entre a carga da bomba e a perda de carga da tubulação. Este residualmétodo é então usado no get_operating_pointmétodo para realmente resolver o ponto de operação do sistema. O método envolve a root_scalarfunção do scipy, que atua como uma interface para vários algoritmos de solução de raiz. Nós vamos deixar oroot_scalarA função escolhe o algoritmo que achar melhor, mas para ajudá-la, especificaremos um intervalo entre colchetes que sabemos que a raiz está entre. Em nosso caso, esse intervalo de variação é os limites superior e inferior da vazão dos dados da curva da bomba. A documentação completa sobre a root_scalarfunção pode ser encontrada aqui .

Dica profissional: o processo de injetar os objetos Pipee na classe (ao contrário de fazer com que a classe do sistema crie um objeto and na instanciação) é chamado de “injeção de dependência”. Isso normalmente é considerado uma boa prática de codificação, pois torna o código mais modular, extensível e mais fácil de depurar/testar.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')

Exploração de projeto

Como um pequeno exemplo final, destacando os benefícios de configurar o código da maneira que fizemos, faremos uma exploração do design. Usando a mesma bomba, gostaríamos de entender os impactos que o comprimento da tubulação tem na vazão volumétrica do sistema. Para fazer isso, simplesmente percorremos uma matriz de comprimentos de tubos (variando de 100 a 1000 pés), atualizamos o atributo de comprimento do Pipeobjeto armazenado no Systeme, em seguida, recalculamos o ponto de operação do sistema, anexando o resultado a um Lista. Finalmente, plotamos a taxa de fluxo de água em função do comprimento do tubo.

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

Conclusão

Este artigo, embora focado fortemente em um problema de exemplo específico de domínio, destaca alguns aspectos de um fluxo de trabalho que acabo usando muito. O problema da análise do ponto operacional surge constantemente na engenharia e na ciência e, embora existam muitas maneiras de abordar o problema, alguns métodos são mais robustos, extensíveis e flexíveis do que outros. A metodologia (formulação de problemas e princípios de estruturação de código) usada neste artigo me serviu incrivelmente bem e espero que outros se inspirem a adotar um fluxo de trabalho semelhante!

Sinta-se à vontade para deixar quaisquer comentários ou perguntas que possa ter ou conecte-se comigo no Linkedin - ficarei mais do que feliz em esclarecer quaisquer pontos de incerteza. Por fim, encorajo você a brincar com o código por conta própria (ou até mesmo usá-lo como um modelo inicial para seus próprios fluxos de trabalho) — o Jupyter Notebook para este artigo pode ser encontrado no meu Github .