Usare Python per risolvere uno dei problemi più comuni in ingegneria
Creazione di un framework generico per l'analisi del punto operativo
Alcune classi di problemi emergono frequentemente in ingegneria. Il focus di questo articolo è su un tipo specifico di problema che si presenta così spesso nel mio lavoro quotidiano che ho pensato di condividere come lo risolvo usando Python. Di che tipo di problema stiamo parlando? Il problema di risolvere il punto di funzionamento di un sistema! Illustriamo cosa intendo con un semplice esempio prima di immergerci in qualcosa di un po' più complesso con il codice.
Vorremmo risolvere il punto di funzionamento del semplice circuito mostrato sotto. Questo può essere fatto riorganizzando la legge di Ohm (V=IR) per isolare la corrente in termini di tensione e resistenza di ingresso note.
Semplice, vero? Sfortunatamente, la maggior parte dei problemi del mondo reale non sono mai così facili. Ad esempio, se ti dicessi che mentre il resistore si riscalda, il suo valore di resistenza cambia, essenzialmente rendendo la resistenza una funzione della corrente. Ci ritroviamo con un'equazione della seguente forma:
Senza conoscere l'effettiva forma funzionale della resistenza, non possiamo semplicemente risolvere la corrente isolandola algebricamente. Inoltre, cosa succede se l'equazione è complicata e non è possibile isolare la corrente da sola? Oppure, forse la resistenza è data in termini di corrente come dati discreti tabulati - quindi non avremmo nemmeno un'espressione algebrica da manipolare per cercare di risolvere la corrente. Come faremmo allora per determinare la corrente nel circuito? Abbiamo bisogno di un approccio più generale per risolvere questo problema.
La soluzione generale a un problema come questo è porlo come un problema di ricerca della radice. Questo è in realtà incredibilmente facile da fare: dobbiamo letteralmente sottrarre il lato destro dell'equazione dal lato sinistro, in modo tale da ottenere un'equazione uguale a zero. In questo modo si ottiene quanto segue:
In questo modo, abbiamo riproposto il nostro problema. Invece di risolvere la corrente direttamente in termini di tutte le altre variabili, possiamo provare a trovare il valore della corrente che può essere inserito nella parte sinistra dell'equazione per farla valutare a zero. Perché formuliamo il problema in questo modo? Perché esistono un sacco di algoritmi numerici (metodo della bisezione, metodo di Newton, ecc.) per risolvere esattamente questo tipo di problema! E alla maggior parte degli algoritmi non importa quanto sia complicato il lato sinistro dell'equazione — non deve nemmeno avere una forma algebrica chiusa (cioè potrebbe essere composto da dati discreti interpolati, integrali valutati numericamente o letteralmente qualsiasi tipo di funzione di complessità arbitraria da valutare). Finché possiamo porre il nostro problema nella forma di f(x)=0,
Il resto di questo articolo illustrerà un esempio di come applicare la metodologia di ricerca delle radici a un problema del mondo reale leggermente più complicato, con enfasi sulla strutturazione del codice del suono e sulle tecniche di organizzazione in Python. Sebbene il problema (determinare la portata dell'acqua in un sistema di tubi/pompe) sia in qualche modo specifico del dominio, la metodologia e le tecniche di codifica utilizzate sono del tutto generali e applicabili a tutti i domini dell'ingegneria. Con questo in mente, cercherò di mantenere gli aspetti di modellazione fisica del problema ad un livello elevato, in modo tale che, indipendentemente dal proprio background tecnico, gli obiettivi di apprendimento primari dell'articolo risultino ancora chiaramente.
Come nota a margine, la mia "specialità" di dominio in questi giorni risiede nel regno dei controlli motore e dell'elettronica di potenza, e sono molto lontano dalle applicazioni di pompaggio/tubazioni. Non tocco l'argomento da anni, ma ho pensato che sarebbe stato un esempio interessante dell'argomento in questione. Sono sicuro che ci sono molte persone là fuori che sono molto più qualificate di me per parlare delle specifiche della modellazione di pompe/tubazioni, ma il mio intento con questo articolo è sulla metodologia, non su come risolvere i problemi di tubazioni/pompe. Indipendentemente da ciò, accolgo apertamente commenti o suggerimenti per il miglioramento da parte di coloro che sono più esperti del settore!
Il problema
Vorremmo trasferire l'acqua da un serbatoio all'altro. Abbiamo già una pompa e delle tubazioni che possono essere utilizzate per collegare i due serbatoi e vogliamo avere una stima di quanto tempo ci vorrà per trasferire tutta l'acqua. Il volume di ciascun serbatoio è noto, quindi se possiamo stimare la portata dell'acqua tra i serbatoi, possiamo stimare quanto tempo impiegherà il processo di trasferimento. Di seguito è riportato l'apparato completo.
Questo problema specifico (che può essere classificato come un problema di "flusso interno") è molto ben compreso nel campo dell'ingegneria meccanica. Per quelli meno familiari, o che necessitano di una rapida revisione, il modo in cui in genere risolviamo questi problemi è con l'equazione di Bernoulli (mostrata sotto).
L'equazione di Bernoulli è essenzialmente un'affermazione di conservazione dell'energia che ci dice come l'energia di una particella fluida viene trasformata tra diversi meccanismi energetici mentre il fluido attraversa una linea di flusso (il percorso di flusso che una particella immaginaria seguirebbe se lasciata cadere nel fluido). Il lato sinistro dell'equazione rappresenta l'energia totale per peso di una particella fluida in qualsiasi prima posizione arbitraria (posizione 1) all'interno del fluido ed è la somma di un termine potenziale gravitazionale, termine cinetico e termine pressione. Mentre il fluido attraversa il sistema, l'energia deve essere conservata, e quindi l'energia totale in qualsiasi secondo punto arbitrario (posizione 2) lungo la linea di flusso (rappresentata dal lato destro dell'equazione) deve essere uguale all'energia totale nella posizione 1 .
La forma sopra dell'equazione di Bernoulli è conosciuta come la forma "testa" dell'equazione perché ogni termine ha unità di lunghezza/altezza. Questo è conveniente per la nostra intuizione perché stiamo essenzialmente equiparando l'energia di ogni termine all'energia potenziale gravitazionale equivalente di una colonna di fluido con un'altezza della testa data. Una delle principali limitazioni dell'equazione di Bernoulli, tuttavia, è che presuppone che non ci siano perdite nel sistema (il che non è un grande presupposto). Per superare questa limitazione, possiamo integrare l'equazione con due termini aggiuntivi come segue:
I termini Hp(Q) e Hl(Q) rappresentano rispettivamente la prevalenza aggiunta al sistema da una pompa e la prevalenza persa nel sistema a causa di effetti del mondo reale (come attrito, viscosità, ecc.). Si noti che entrambi i termini sono funzioni della portata del fluido del sistema, Q. (Come interessante conseguenza del paragrafo precedente che descrive l'interpretazione della prevalenza, la prevalenza della pompa indica quanto in alto una pompa potrebbe teoricamente spingere un fluido). Esamineremo i termini di pompaggio e perdita in modo più approfondito tra un po', ma prima di farlo, semplifichiamo l'equazione precedente per il nostro problema specifico.
Guardando di nuovo il sistema sopra, sceglieremo convenientemente le nostre due posizioni per l'equazione di Bernoulli, in modo tale che la maggior parte dei termini si annulli. Possiamo farlo scegliendo rispettivamente le posizioni 1 e 2 sulla superficie libera dell'acqua di ciascun serbatoio, dove la pressione è costante e uguale alla pressione atmosferica (P1=P2) e la velocità è approssimativamente costante e zero (V1 =V2=0). Supponiamo inoltre che l'altezza dell'acqua nei due serbatoi sia la stessa nell'istante in cui stiamo analizzando il sistema tale che Z1=Z2. Dopo aver semplificato l'algebra, vediamo che quasi tutti i termini si annullano e rimane il fatto che la prevalenza prodotta dalla pompa deve essere uguale alla prevalenza persa nel sistema a causa delle non idealità. Detto in altro modo, la pompa sta compensando eventuali perdite di energia nel sistema.
Questa situazione può essere vista qualitativamente nella figura sottostante. La prevalenza prodotta da una pompa diminuisce all'aumentare della portata, mentre le perdite in un sistema di tubazioni aumentano all'aumentare della portata. Il punto di intersezione delle due curve (prevalenza pompa = perdita di carico) determina il punto di lavoro (portata) dell'impianto.
L'ultimo passo prima di poter entrare nel codice è porre il problema come un problema di ricerca della radice. Sottraendo il membro destro dell'equazione dal membro sinistro, otteniamo il problema di root solving che stiamo cercando. Cioè, abbiamo posto il nostro problema come segue: trova la portata (Q) tale che il lato sinistro dell'equazione sottostante sia uguale a zero. A questo punto, la prevalenza della pompa sarà uguale alle perdite di carico dell'impianto.
Il codice
Per evitare di perdere il quadro generale di ciò che stiamo facendo, non spiegherò ogni piccolo dettaglio del codice (presumo che tu abbia già un discreto background in Python). Invece, concentrerò i miei sforzi sull'assicurare che la narrativa e la strutturazione del codice siano chiare e, se necessario, approfondirò i dettagli. Come sempre, sentiti libero di fare qualsiasi domanda se qualcosa non è chiaro.
Impostare
Inizieremo importando tutti i moduli necessari. Diventerà evidente come ciascuno dei moduli viene utilizzato in seguito, ma vale la pena notare che le dichiarazioni di importazione chiave sono quelle di scipy. Queste sono le funzioni specifiche per il problema in questione. Il blocco di codice imposta anche alcune impostazioni di stampa predefinite (secondo il gusto personale), crea una cartella in cui salvare le cifre generate e definisce alcune costanti di conversione delle unità che ci semplificano la vita più avanti nel codice.
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
Il passo successivo è modellare le perdite di carico del tubo (il termine Hl(Q) nell'equazione di Bernoulli estesa sopra). Questo viene in genere fatto utilizzando l'equazione di Darcy-Weisbach mostrata di seguito, dove f è un fattore di attrito (ne parleremo tra poco), v è la velocità del flusso, g è la gravità e L e D sono rispettivamente la lunghezza e il diametro del tubo.
Sfortunatamente, il fattore di attrito (f) non è costante, ma dipende anche dalla velocità del flusso, dalle proprietà del fluido e dalle dimensioni del tubo. Esistono vari modelli per calcolare f, ma useremo l'equazione di Haaland, mostrata sotto.
In questa equazione, epsilon è la rugosità superficiale del tubo (che può essere trovata nelle tabelle dei manuali di ingegneria) e Re è il famoso numero di Reynolds, calcolato come di seguito.
Infine, possiamo notare che il volume spostato per unità di tempo, o portata volumetrica (Q), è uguale all'area della sezione trasversale (A) del tubo moltiplicata per la velocità del flusso (v). Pertanto, data una portata nel tubo, possiamo calcolare la corrispondente velocità del flusso nel tubo come:
Si spera che tutte queste equazioni non stiano sminuendo il quadro più ampio: stiamo solo osservando un particolare modello di calcolo della perdita di carico in un tubo. Data una portata e le dimensioni del tubo, calcolare prima la velocità del flusso corrispondente, quindi collegare le equazioni precedenti per calcolare la perdita di carico del tubo. Questo è esattamente ciò che Pipeimplementa la classe (mostrata sotto).
Il metodo di inizializzazione memorizza le dimensioni dei tubi (tutti assunti in metri) e le proprietà del fluido. Il Ametodo calcola l'area della sezione trasversale del tubo (per chi non ha familiarità con il @propertydecoratore, questo articolo lo spiega molto bene). Il Q_to_vmetodo converte la portata in galloni al minuto (gpm) in una velocità di flusso in m/s. Il friction_factormetodo valuta l'equazione di Haaland come descritto sopra, head_losse head_loss_feetvaluta la perdita di carico del tubo rispettivamente in metri e piedi (usando l'equazione di 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')
Abbiamo un modello funzionante delle perdite di carico del tubo — ora abbiamo bisogno di un modello del carico prodotto da una pompa, Hp(Q). Sono sicuro che ci sono modelli analitici che possono essere utilizzati per determinare il comportamento di una pompa, ma supponiamo di avere già una pompa specifica, vale a dire la seguente pompa di potenza frazionaria casuale che ho trovato online :
La maggior parte delle pompe avrà una scheda tecnica che include una curva della pompa corrispondente che caratterizza il comportamento delle pompe. Per la pompa sopra, otteniamo la seguente curva della pompa:
A questo punto abbiamo un'immagine che descrive il comportamento delle pompe, ma non un modello matematico che può effettivamente essere utilizzato per determinare come si comporterà nel sistema. Questo problema si presenta continuamente e il modo in cui lo risolvo è 1) digitalizzare i dati, quindi 2) avvolgere i dati discreti con uno schema di interpolazione per produrre una funzione continua. Lasciatemi illustrare.
Passaggio 1) Esistono molti strumenti per digitalizzare i dati delle immagini: il mio preferito è l'utilità online gratuita WebPlotDigitizer . Si carica l'immagine del grafico di interesse, si allineano gli assi e quindi si estraggono le curve di dati desiderate (manualmente o con lo strumento di estrazione automatica). I dati possono quindi essere esportati in un file .csv.
Passaggio 2) Ora che abbiamo i dati digitalizzati, dobbiamo solo avvolgerli con una sorta di interpolatore: questo è esattamente ciò che Pipefa la classe di seguito. Il metodo di inizializzazione accetta il nome del file .csv, memorizza il nome del file, carica i dati in un DataFrame panda, memorizzandoli in un dataattributo e quindi passa i dati alla interp1dfunzione da scipy. La interp1dfunzione genera quindi una nuova funzione che, per impostazione predefinita, utilizza l'interpolazione lineare per trasformare i punti dati discreti in una funzione continua (la documentazione completa per la interp1dfunzione è disponibile qui ). La funzione di interpolazione appena generata viene quindi memorizzata in un _interpattributo per un successivo accesso. La Pipeclasse contiene anche aboundsmetodo che restituisce un elenco contenente i valori min/max della portata nei dati della curva della pompa (questo verrà utilizzato nell'algoritmo di ricerca della radice) e un head_gain_feetmetodo che prende la portata in gpm e chiama la funzione di interpolazione sottostante che è stato generato da 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 disponiamo dell'infrastruttura per risolvere il punto di funzionamento del sistema pompa/tubazione. L'ultimo passaggio consiste nel creare una Systemclasse che accetti un oggetto Pipee Pumped esegua l'operazione di risoluzione delle radici. Come si può vedere nel codice seguente, la Systemclasse accetta e memorizza un oggetto Pipeand Pump. Quindi utilizza i due oggetti per creare un residualmetodo che calcola la differenza tra la prevalenza della pompa e la perdita di carico del tubo. Questo residualmetodo viene quindi utilizzato nel get_operating_pointmetodo per risolvere effettivamente il punto operativo del sistema. Il metodo esegue il wrapping della root_scalarfunzione da scipy, che funge da interfaccia per vari algoritmi di risoluzione delle radici. Lasceremo ilroot_scalarfunzione sceglie l'algoritmo che ritiene più adatto, ma per aiutarlo, specificheremo un intervallo di parentesi tra cui sappiamo che si trova la radice. Nel nostro caso, questo intervallo di parentesi rappresenta i limiti di portata superiore e inferiore dei dati della curva della pompa. La documentazione completa sulla root_scalarfunzione può essere trovata qui .
Suggerimento pro: il processo di inserimento degli oggetti Pipee nella classe (invece di fare in modo che la classe di sistema crei un oggetto e all'istanza) è chiamato "iniezione di dipendenza". Questa è in genere considerata una buona pratica di codifica in quanto rende il codice più modulare, estensibile e più facile da eseguire il debug/test.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')
Esplorazione del design
Come ultimo piccolo esempio, evidenziando i vantaggi dell'impostazione del codice nel modo in cui lo abbiamo fatto, eseguiremo un'esplorazione del progetto. Utilizzando la stessa pompa, vorremmo comprendere gli impatti che la lunghezza del tubo ha sulla portata volumetrica nel sistema. Per fare ciò, eseguiamo semplicemente il loop su una serie di lunghezze di tubo (che vanno da 100 a 1000 piedi), aggiorniamo l'attributo di lunghezza Pipedell'oggetto memorizzato in System, quindi ricalcoliamo il punto operativo del sistema, aggiungendo il risultato a un elenco. Infine tracciamo la portata dell'acqua in funzione della lunghezza del 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')
Conclusione
Questo articolo, sebbene incentrato fortemente su un problema di esempio specifico del dominio, evidenzia alcuni aspetti di un flusso di lavoro che finisco per utilizzare molto. Il problema dell'analisi del punto operativo si presenta costantemente in ingegneria e scienza, e sebbene ci siano molti modi per affrontare il problema, alcuni metodi sono più robusti, estensibili e flessibili di altri. La metodologia (formulazione del problema e principi di strutturazione del codice) utilizzata in questo articolo mi è stata incredibilmente utile e spero che altri siano ispirati ad adottare un flusso di lavoro simile!
Sentiti libero di lasciare qualsiasi commento o domanda che potresti avere o di connetterti con me su Linkedin: sarei più che felice di chiarire eventuali punti di incertezza. Infine, ti incoraggio a giocare tu stesso con il codice (o persino a usarlo come modello di partenza per i tuoi flussi di lavoro): il Jupyter Notebook per questo articolo può essere trovato sul mio Github .
Nicola Hemenway
- Se ti è piaciuto, seguimi su Medium
- Prendi in considerazione la possibilità di iscriverti agli aggiornamenti via e-mail
- Interessato a collaborare? Connettiamoci su LinkedIn
![Che cos'è un elenco collegato, comunque? [Parte 1]](https://post.nghiatu.com/assets/images/m/max/724/1*Xokk6XOjWyIGCBujkJsCzQ.jpeg)
