transient tank filling

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.

Input data

First of all, let’s analyze the transient vessel pump system.

Pump vessel system

Therefore, we use the following data.

  • System fluid: Water. It is assumed to be incompressible. Therefore, the density (\rho = 1000 kg/m³) in the system is constant.
  • We assume a total of the pump efficiency \eta = 0.7.
  • The maximum discharge head of the pump h_d = 7 m, while the power input is 8 kW.
  • The diameter of the vessel is D_v = 2 m, while the initial fill level is 2 m.
  • Furthermore, the maximum water level in the vessel is h_{max} = 4.75 m, while the minimum water level h_{min} = 2.75 m.
  • Diameter of the outlet d_{out} = 0.1 m

Equation system

The volume flow of the pump is calculated from the power input P, total efficiency \eta, density \rho, gravitational acceleration g, discharge head (pump) h_d:

(1)   \begin{equation*} \dot V_{in} = \frac{P \cdot \eta}{\rho \cdot g \cdot h_d}. \end{equation*}

Furthermore, we can calculate the outflow velocity c_{out} from the vessel using the Torricelli’s equation:

(2)   \begin{equation*} c_{out}= \sqrt{2 \cdot g \cdot h}. \end{equation*}

Afterwards, the outlet volume flow \dot V_{out} can be calculated using the outlet area A_{out}

(3)   \begin{equation*} \dot V_{out}= c_{out} \cdot A_{out}. \end{equation*}

Finally, the time dependent change in height h is calculated using the vessel area A_v and the inlet \dot V_{in} and outlet volume flow \dot V_{out}:

(4)   \begin{equation*} \frac{dh}{dt} = \frac{V_{in} - V_{out}}{A_v}. \end{equation*}

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 h_{max} 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 h_{min} 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 h 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
#Author: johannes <>, 22.12.17
#License: MIT License (

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.integrate import odeint

# gravitational acceleration
g = 9.81

def area(d):
    return math.pi*d**2/4.

def toricelliVelocity(h):
    return (2.*g*h)**0.5

def volumeFlowPump(P, rho, eta, h):
    return (P*eta)/(rho*g*h)

def vesselPumpSystem(h, t):
    global hflag
    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)
    return dhdtVessel
# Initial height in the vessel
h0 = 2 #m

# Sampling
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 (\dot V = 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.

Transient fill level of the vessel


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.