Apply an Integrator to a Vessel Pump System using Python
In the following blog post we will have a look at a transient vessel pump system. It is meant to give an insight into transient systems including transition conditions. Furthermore, we will have a glance at integrators.
First of all, let’s analyze the transient vessel pump system.
Therefore, we use the following data.
- System fluid: Water. It is assumed to be incompressible. Therefore, the density ( = 1000 kg/m³) in the system is constant.
- We assume a total of the pump efficiency = 0.7.
- The maximum discharge head of the pump = 7 m, while the power input is 8 kW.
- The diameter of the vessel is = 2 m, while the initial fill level is 2 m.
- Furthermore, the maximum water level in the vessel is = 4.75 m, while the minimum water level = 2.75 m.
- Diameter of the outlet = 0.1 m
The volume flow of the pump is calculated from the power input , total efficiency , density , gravitational acceleration , discharge head (pump) :
Furthermore, we can calculate the outflow velocity from the vessel using the Torricelli’s equation:
Afterwards, the outlet volume flow can be calculated using the outlet area
Finally, the time dependent change in height is calculated using the vessel area and the inlet and outlet volume flow :
Behavior of the vessel pump system
The pump conveys water from an infinite reservoir into a vessel. Simultaneously, water is drained from the vessel via the outlet. Since the inlet volume flow is larger than the outlet volume flow the water level in the vessel rises.
When the water level reaches the maximum level a signal is sent to the control valve and the pump via the level switch high high (LSHH). Due to that, the control valve shuts and the pump stops. Hence, the level in the vessel drops.
When the minimum level is reached, the level switch high low (LSHL) sends a signal to the control valve and pump. Therefore, the control valve opens while the pump starts and the level in the vessel rises.
Solving the Differential Equation using an Integrator
We can describe the time dependent fill level using a differential equation. Hence, it can be solved, using an integration method. Also, a very basic integration method is the (explicit) Euler method.
In this example we use the integrator package odeint from the scipy library. It is a package, where an adequate integrator is chosen automatically. Basically, odeint is using various backward differentiation formula for stiff differential equations. In addition, for non-stiff differential equations multistep methods.
Writing the Python script
So let us go over the Python script.
- Import the dependencies e.g. numpy
- Define the some basic functions like area, Toricelli equation and the volume flow
- Specify the system function. Therefore, input variables are set. Furthermore, the transition condition for the set points are set. Also, the equations are called and defined.
- Define various sampling points (discretization)
- Integrate and therefore solve the equation system
- Finally, plot the results
#License: MIT License (http://opensource.org/licenses/MIT)
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# gravitational acceleration
g = 9.81
def volumeFlowPump(P, rho, eta, h):
def vesselPumpSystem(h, t):
eta = 0.7
rho = 1000 #kg/m^2
dOut = 0.1 #m
dVessel = 2 #m
powerPump = 8000 #W
hDischarge = 7 #m
hMaxVessel = 4.75 #m
hMinVessel = 2.75 #m
if h >= hMaxVessel:
P = 0 #W
hflag = True
elif h >= hMinVessel and hflag:
P = 0 #W
P = powerPump
hflag = False
volFlowIn = volumeFlowPump(P, rho, eta, hDischarge)
volFlowout = toricelliVelocity(h)*area(dOut)
dhdtVessel = (volFlowIn-volFlowout)/area(dVessel)
# Initial height in the vessel
h0 = 2 #m
t2 = np.linspace(0, 1500, 50)
t1 = np.linspace(0, 1500, 500)
t3 = np.linspace(0, 1500, 5000)
# Start integration
y1 = odeint(vesselPumpSystem, h0, t1)
y2 = odeint(vesselPumpSystem, h0, t2)
y3 = odeint(vesselPumpSystem, h0, t3)
# Plot results
plt.plot(t2, y2, label="50 samples", linewidth=1.5, color='b')
plt.plot(t1, y1, label="500 samples", linewidth=1.5, color='g')
plt.plot(t3, y3, label="5000 samples", linewidth=1.5, color='r')
plt.xlabel('Time in sec')
plt.ylabel('Fill level in m')
As we can see, the vessel is filled until the maximum fill level is reached. Afterwards, the inlet volume flow stops. Hence, the vessel is emptied until the minimum fill level is reached. Afterwards, the vessel fill level varies between the maximum and minimum fill level.
In addition, let us observe the results regarding the variegated sample steps (discretization) of the integrator. We can see, that the results differ slightly. This is based on the fact, that the width of the integration step (distance between two sample points) differs. As a result, the transition condition ( = 0) can occur at slightly different moments in time. Therefore, the results may differ. Furthermore, the gradient information is propagated through time. Hence, differences in the solutions accumulate over time.
The transient vessel pump system can be described using a mixture of algebraic equations and a differential equation. Therefore, it is solved numerically using integration methods. Furthermore, transition conditions are included to start and stop the feed volume flow.
In order to solve the differential equation system of the vessel pump system, a suitable integrator needs to be chosen. Therefore, integrator packages that switch automatically between different integrators are helpful. Especially, when there is a lack of in-depth knowledge. Furthermore, the amount of discretization steps especially relevant because they can impact the solution. As a result, I would recommend an automatic stepper, which adapts the step size automatically.
What did you think?
I’d like to hear what you think about this post.
Let me know by leaving a comment below.