Making a generic framework for working level evaluation
Sure lessons of issues come up incessantly in engineering. The main target of this text is on a selected kind of downside that comes up so usually in my every day work that I figured I‘d share how I am going about fixing it utilizing Python. What kind of downside are we speaking about? The issue of fixing the working level of a system! Let’s illustrate what I imply with a easy instance earlier than diving into one thing a little bit extra advanced with code.
We want to resolve the working level of the easy circuit proven under. This may be performed by rearranging Ohm’s regulation (V=IR) to isolate the present by way of the recognized enter voltage and resistance.
Easy, proper? Sadly, most actual world issues are by no means this straightforward. For instance, what if I advised you that because the resistor heats up, it’s resistance worth adjustments, primarily making the resistance a operate of present. We find yourself with an equation of the next type:
With out realizing the precise practical type of the resistance, we are able to’t simply resolve for present by isolating it algebraically. Moreover, what if the equation is difficult and it’s not attainable to isolate the present by itself? Or, maybe the resistance is given by way of present as tabulated discrete information — then we wouldn’t even have an algebraic expression to govern to attempt to resolve for present. How would we go about figuring out the present within the circuit then? We want a extra basic strategy to fixing this downside.
The overall answer to an issue like that is to pose it as a root discovering downside. That is really extremely straightforward to do — we actually simply need to subtract the precise hand facet of the equation from the left hand facet, such that we get an equation that equals zero. Doing so yields the next:
By doing this, we’ve re-posed our downside. As a substitute of fixing for the present immediately by way of all different variables, we are able to attempt to discover the worth of present that may be enter into the left hand facet of the equation to make it consider to zero. Why can we formulate the issue like this? As a result of there are a ton of numerical algorithms that exist (bisection technique, Newton’s technique, and many others.) to resolve this precise sort of downside! And many of the algorithms don’t care how difficult the left hand facet of the equation is — it doesn’t even need to have a closed algebraic type (i.e. it could possibly be composed of interpolated discrete information, numerically evaluated integrals, or actually any kind of operate of arbitrary complexity to guage). So long as we are able to pose our downside within the type of f(x)=0, we are able to (virtually) at all times discover a answer to the issue (and the code might be modified/prolonged simply if the issue assertion adjustments — no have to re-do algebra).
The remainder of this text goes to stroll by means of an instance of learn how to apply the foundation discovering methodology to a barely extra difficult, real-world downside, with emphasis on sound code structuring and group methods in Python. Though the issue (figuring out the stream price of water in a pipe/pump system) is considerably area particular, the methodology and coding methods which are used are utterly basic, and relevant to all engineering domains. With this in thoughts, I’ll attempt to preserve the bodily modeling features of the issue to a excessive degree, such that no matter ones technical background, the first studying goals of the article nonetheless come by means of clearly.
As a facet notice, my area “specialty” as of late lies within the realm of motor controls and energy electronics, and I’m very far faraway from pumping/piping functions. I haven’t touched on the topic in years, however thought it will make for an fascinating instance of the subject at hand. I’m certain there are many people on the market which are way more certified to talk on the specifics of pump/pipe modeling than myself, however my intent with this text is on the methodology — not on learn how to resolve pipe/pump issues. Regardless, I overtly welcome feedback or options for enchancment from these which are extra educated of the sector!
We want to switch water from one tank to a different. We have already got a pump and a few piping that can be utilized to attach the 2 tanks and wish to get an estimate of how lengthy it is going to take to switch all the water. The quantity of every tank is thought, so if we are able to estimate the stream price of water between the tanks, we are able to estimate how lengthy the switch course of will take. The complete equipment is proven under.
This particular downside (which might be categorised as an “inner stream” downside) may be very effectively understood throughout the subject of mechanical engineering. For these much less acquainted although, or in want of a fast evaluation, the best way we sometimes go about fixing these issues is with the Bernoulli equation (proven under).
The Bernoulli equation is basically an power conservation assertion that tells us how the power of a fluid particle is reworked between totally different power mechanisms because the fluid traverses alongside a streamline (the stream path that an imaginary particle would observe if dropped within the fluid). The left hand facet of the equation represents the entire per-weight power of a fluid particle at any arbitrary first location (location 1) throughout the fluid, and is the sum of a gravitational potential time period, kinetic time period, and strain time period. Because the fluid traverses the system, power should be conserved, and thus the entire power at any arbitrary second level (location 2) alongside the streamline (represented by the precise hand facet of the equation) should be equal to the entire power at location 1.
The above type of the Bernoulli equation is named the “head” type of the equation as a result of every time period has items of size/top. That is handy for our instinct as a result of we’re primarily equating the power of every time period to the equal gravitational potential power of a column of fluid with a top of the given head. One main limitation of the Bernoulli equation, nevertheless, is that it assumes there are not any losses within the system (which isn’t an incredible assumption). To beat this limitation, we are able to complement the equation with two further phrases as follows:
The Hp(Q) and Hl(Q) phrases characterize the pinnacle added to the system by a pump, and the pinnacle misplaced within the system from actual world results (like friction, viscosity, and many others.) respectively. Be aware that each phrases are capabilities of the system’s fluid stream price, Q. (As an fascinating consequence of the above paragraph describing the interpretation of head, the pinnacle of pump tells you the way excessive a pump may theoretically push a fluid). We’ll look at the pump and loss phrases extra totally in a bit, however earlier than we do, let’s simplify the above equation for our particular downside.
Trying on the system above once more, we are going to conveniently select our two areas for the Bernoulli equation, such that many of the phrases cancel. We will do that by selecting areas 1 and a pair of to be on the free floor of the water of every tank respectively, the place the strain is fixed and equal to atmospheric strain (P1=P2), and the speed is roughly fixed and 0 (V1=V2=0). We may even assume that the peak of the water within the two tanks is similar on the instantaneous we’re analyzing the system such that Z1=Z2. After simplifying the algebra, we see that just about all the phrases cancel, and we’re left with the truth that the pinnacle produced by the pump should equal the pinnacle misplaced within the system on account of non-idealities. Acknowledged one other means, the pump is making up for any power losses within the system.
This case might be seen qualitatively within the determine under. The pinnacle produced by a pump decreases with growing stream price, whereas, the losses in a piping system improve with growing stream price. The purpose the place the 2 curves intersect (pump head = head loss) determines the working level (stream price) of the system.
The final step earlier than we are able to bounce into code is to pose the issue as a root discovering downside. Subtracting the precise hand facet of the equation from the left hand facet, we acquire the foundation fixing downside that we’re searching for. That’s, we’ve posed our downside as follows: discover the stream price (Q) such that the left hand facet of the equation under equals zero. At this level, the pump head will equal the system head losses.
To keep away from dropping the large image of what we’re doing, I’m not going to clarify each little element of the code (I’m assuming you may have an affordable background in Python already). As a substitute, I’m going to focus my efforts on guaranteeing the narrative and structuring of the code is evident, and go into extra element as wanted. As at all times, be at liberty to ask any questions if one thing isn’t clear.
Setup
We are going to start by importing all the needed modules. It should grow to be obvious how every of the modules is used in a while, however it’s value noting that the important thing import statements are those from scipy. These are the capabilities which are particular to the issue at hand. The code block additionally set’s some default plotting settings (to private style), creates a folder to avoid wasting the generated figures in, and defines some unit conversion constants that make our life simpler later within the code.
from dataclasses import dataclass
from pathlib import Pathimport numpy as np
import matplotlib.pyplot as plt
import pandas as pd
#these are the important thing libraries for fixing the issue
from scipy.interpolate import interp1d
from scipy.optimize import root_scalar
#set plotting defaults
plt.fashion.use('seaborn-v0_8-darkgrid')
plt.rcParams['font.family'] = 'Instances New Roman'
plt.rcParams['font.size'] = 12
figsize = (6.4,4)
#make folder to avoid wasting plots to
plots_folder = Path('plots')
plots_folder.mkdir(exist_ok=True)
#outline 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
Subsequent we are going to make a Python dataclass that primarily acts as a container for storing fluid properties (density, viscosity, and gravity). By default the properties are set to these of water. Be aware that though not utterly needed, Python information lessons are extremely handy. In the event that they’re new to you, I extremely suggest you try this video.
@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
Pipe Mannequin
The following step is to mannequin the pinnacle losses of the pipe (the Hl(Q) time period within the prolonged Bernoulli equation above). That is sometimes performed utilizing the Darcy-Weisbach equation proven under, the place f is a friction issue (extra on this shortly), v is the stream velocity, g is gravity, and L and D are the pipe size and diameter respectively.
Sadly, the friction issue (f) isn’t fixed, but additionally is determined by the stream velocity, fluid properties, and pipe dimensions. Varied fashions exist to compute f, however we are going to use Haaland’s equation, proven under.
On this equation, epsilon is the floor roughness of the pipe (which might be present in engineering textbook tables), and Re is the well-known Reynolds quantity, computed as per under.
Lastly, we are able to notice that the swept quantity per unit time, or volumetric stream price (Q) equals the cross-sectional space (A) of the pipe instances the stream velocity (v). Thus, given a stream price within the pipe, we are able to compute the corresponding stream velocity within the pipe as:
Hopefully all of those equations aren’t detracting from the larger image— we’re simply taking a look at one explicit mannequin of computing head loss in a pipe. Given a stream price and pipe dimensions, first compute the corresponding stream velocity, then plug by means of the equations above to compute the pipe’s head loss. That is precisely what the Pipe class (proven under) implements.
The initialization technique shops the pipes dimensions (all assumed to be in meters) and fluid properties. The A technique computes the cross-sectional space of the pipe (for these unfamiliar with the@property decorator, this article explains it very effectively). The Q_to_v technique converts the stream price in gallons per minute (gpm) to a stream velocity in m/s. The friction_factor technique evaluates Haaland’s equation as described above, and the head_loss and head_loss_feet consider the pipe’s head loss in meters and toes respectively (utilizing the Darcy-Weisbach equation).
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 space of pipe in m^2"""
return np.pi*(self.D/2)**2 #space in m^2
def Q_to_v(self, gpm):
"""Converts gpm to fluid velocity in pipe in m/s"""
Q = gpm*GALLONS_TO_M3/60 #stream price in m^3/s
return Q/self.A #stream velocity in m/s
def friction_factor(self, gpm):
"""computes Darcy friction issue, given stream price in gpm
This technique makes use of Haaland's equation, wich is an express approximation
of the well-known, however implicit Colebrook equation
"""
#first get stream velocity from stream price and pipe dimensions
v = self.Q_to_v(gpm)
#compute Reynold's quantity
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 stream price in gpm"""
#get stream velocity
v = self.Q_to_v(gpm)
#get Darcy friction issue
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 toes, given stream price in gpm"""
hl_meters = self.head_loss(gpm)
return hl_meters/FEET_TO_METERS
Let’s see the pipe class in motion. First, we are able to create a Fluid object for water and a Pipe object that’s 100 toes lengthy and 1.25 inches in diameter.
#create fluid object for water
water = Fluid()#create pipe section 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)
Subsequent we’ll plot the pinnacle loss curve as a operate of the stream price in gpm. Isn’t it unimaginable how clear and simple the code under turns into after we exploit object oriented programming?
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('Movement Price [gpm]')
ax.set_ylabel('Head Loss [ft]')
fig.tight_layout()
fig.savefig(plots_folder/'pipe_loss_curve.png')
Pump Mannequin
We’ve acquired a working mannequin of the pipe head losses — now we’d like a mannequin of the pinnacle produced by a pump, Hp(Q). I’m certain there are analytic fashions that can be utilized to find out a pumps conduct, however we are going to assume that we have already got a selected pump — specifically the next random, fractional horspower pump that I discovered on-line:
Most pumps may have an information sheet that features a corresponding pump curve that characterizes the pumps conduct. For our pump, we get the pump curve under (notice that for copyright causes, it is a recreation of the determine supplied by the producer — the unique might be discovered right here).
At this level we’ve a picture depicting the pumps conduct, however not a mathematical mannequin that may really be used to find out the way it will carry out within the system. This challenge comes up on a regular basis, and the best way I get round it’s to 1) digitize the information, after which 2) wrap the discrete information with an interpolation scheme to supply a continous operate. Let me illustrate.
Step 1) There are various instruments that exist to digitize picture information — my private favourite is the free on-line utilitiy WebPlotDigitizer. You load within the plot picture of curiosity, align the axes, after which extract the specified information curves (both manually, or with the automated extraction software). The information can then be exported to a .csv file.
Step 2) Now that we’ve acquired digitized information, we simply have to wrap it with some kind of interpolator — that is precisely what the Pipe class does under. The initialization technique takes within the .csv file identify, shops the file identify, masses the information right into a pandas DataFrame, storing it in a information attribute, after which passes the information to the interp1d operate from scipy. The interp1d operate then generates a brand new operate that, by default, makes use of linear interpolation to show the discrete datapoints right into a continous operate (full documentation for interp1d operate might be discovered right here). The newly generated interpolation operate is then saved in a _interp attribute for later entry. The Pipe class additionally incorporates a bounds technique that returns a listing containing the min/max values of stream price within the pump curve information (this shall be used within the root discovering algorithm), and a head_gain_feet technique that takes within the stream price in gpm, and calls the underlying interpolation operate that was generated by interp1d.
class Pump():def __init__(self, file):
#retailer file identify
self.file = file
#learn information into pandas dataframe and assign column names
self.information = pd.read_csv(file, names=['gpm', 'head [ft]']).set_index('gpm')
#create steady interpolation operate
self._interp = interp1d(self.information.index.to_numpy(), self.information['head [ft]'].to_numpy())
@property
def bounds(self):
"""returns min and max stream charges in pump curve information"""
return [self.data.index.min(), self.data.index.max()]
def head_gain_feet(self, gpm):
"""return head (in toes) produced by the pump at a given stream price"""
return self._interp(gpm)
We will create a Pump object and take a look at the uncooked information that we’ve learn in.
pump = Pump('pump_data.csv')
pump.information.head()
We will additionally plot the pump curve information with the pipe loss curve to visually see the place the system will function.
head_loss = [pipe.head_loss_feet(gpm) for gpm in pump.data.index]fig, ax = plt.subplots(figsize=figsize)
ax.plot(pump.information, label='Pump Curve')
ax.plot(pump.information.index, head_loss, label='Pipe Head Loss')
ax.set_xlabel("Movement Price [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')
System Mannequin
We lastly have the infrastructure in place to resolve the working level of the pump/pipe system. The final step is to create a System class that takes in a Pipe and Pump object and performs the foundation fixing operation. As might be seen within the code under, the System class takes in and shops a Pipe and Pump object. It then makes use of the 2 objects to create a residual technique that computes the distinction between the pump head and pipe head loss. This residual technique is then used within theget_operating_point technique to really resolve the working level of the system. The tactic wraps the root_scalar operate from scipy, which acts as an interface for varied root fixing algorithms. We are going to let the root_scalar operate select whichever algorithm it sees greatest match, however to assist it, we are going to specify a bracketing interval that we all know the foundation lies between. In our case, this bracketing interval is the higher and decrease flow-rate bounds of the pump curve information. Full documentation on the root_scalar operate might be discovered right here.
Professional tip: the method of injecting the
PipeandPumpobjects into theSystemclass (versus having the system class create aPipeandPumpobject at instantiation) is known as “dependency injection”. That is sometimes thought of a very good coding follow because it makes code extra modular, extensible, and simpler to debug/take a look at.
class System():def __init__(self, pipe: Pipe, pump: Pump):
self.pipe = pipe
self.pump = pump
def residual(self, gpm):
"""
Computes the distinction between the pinnacle produced by the pump
and the pinnacle loss within the pipe. At regular state, the pump head and
head loss shall be equal and thus the residual operate will go to zero
"""
return self.pump.head_gain_feet(gpm) - self.pipe.head_loss_feet(gpm)
def get_operating_point(self):
"""resolve for the stream price the place the residual operate equals zero.
i.e. the pump head equals the pipe head loss"""
return root_scalar(sys.residual, bracket=pump.bounds)
Let’s create a System and run the get_operating_point technique to look at the fruits of our labor. As might be seen by the code output, the get_operating level technique merely passes out the output of the root_scalar operate, which is a RootResults object. This object is basically only a container that shops varied attributes, an important of which, is the root attribute because it incorporates the answer to our downside.
sys = System(pipe, pump)res = sys.get_operating_point()
res
We will plot the identical pump and head loss curves once more, this time including in a vertical line on the computed steady-state working level.
head_loss = [pipe.head_loss_feet(gpm) for gpm in pump.data.index]fig, ax = plt.subplots(figsize=figsize)
ax.plot(pump.information, label='Pump Curve')
ax.plot(pump.information.index, head_loss, label='Pipe Head Loss')
#plot vertical line at working level
ax.axvline(res.root, colour='okay', ls='--', lw=1)
ax.legend(frameon=True, facecolor='white', framealpha=1, loc=6)
ax.set_xlabel("Movement Price [gpm]")
ax.set_ylabel("Head [ft]")
ax.set_title(f'Working Level = {res.root:.1f} gpm')
fig.tight_layout()
fig.savefig(plots_folder/'intersection_solution.png')
Voila! We’ve programmatically decided the working level of our system. And since we’ve performed it utilizing a considerably generic coding framework, we are able to simply attempt the identical evaluation utilizing totally different pumps or piping! We may even prolong our code to incorporate a number of pumps, or varied pipe fittings/pipe branches.
Design Exploration
As a last little instance, highlighting the advantages of organising the code the best way we did, we are going to carry out a design exploration. Utilizing the identical pump, we want to perceive the impacts that pipe size has on the volumetric stream price within the system. To do that, we merely loop over an array of pipe lengths (starting from 100 to 1000 toes), replace the size attribute of the Pipe object saved within the System, after which re-compute the working level of the system, appending the consequence to a listing. Lastly we plot water stream price as a operate of pipe size.
#sweep pipe size from 100 to 1000 toes
lengths_feet = np.linspace(100, 1000, 1000)
lengths_meters = lengths_feet*FEET_TO_METERSflow_rates = []
for l in lengths_meters:
#replace pipe size
sys.pipe.L = l
#compute new stream price answer
res = sys.get_operating_point()
#append answer to stream charges checklist
flow_rates.append(res.root)
#plot outcomes
fig, ax = plt.subplots(figsize=figsize)
ax.plot(lengths_feet, flow_rates)
ax.set_xlabel("Pipe Size [ft]")
ax.set_ylabel("Movement Price [gpm]")
# ax.set_ylim(backside=0)
fig.tight_layout()
fig.savefig(plots_folder/'flow_vs_pipe_length.png')
In a number of strains of code, we have been capable of extract significant perception into the conduct of our system. If this was a design downside, these insights may drive key design decisions.
Conclusion
This text, though centered strongly on a website particular instance downside, highlights a number of features of a piece stream that I find yourself utilizing lots. The issue of working level evaluation comes up consistently in engineering and science, and though there are various methods to strategy the problem, some strategies are extra sturdy, extensible, and versatile than others. The methodology (downside formulation and code structuring rules) used on this article have served me extremely effectively, and I hope that others are impressed to undertake the same work stream!
Be happy to depart any feedback or questions you might have or join with me on Linkedin — I’d be very happy to make clear any factors of uncertainty. Lastly, I encourage you to mess around with the code your self (and even use it as a beginning template on your personal workflows ) — the Jupyter Pocket book for this text might be discovered on my Github.
