Pressure drop in a circular section pipe simulation. Read the full article at https://firsttimeprogrammer.blogspot.com/2016/07/fluid-dynamics-pressure-drop-modelling.html

pressure_drop_revised.py

# Imports
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from numpy import vectorize

# Let's define the variables
Q = np.linspace(0.01,1,50) # Flow rate m^3/s
d = np.linspace(0.5,1,50)  # Pipe's diameter m
mu = 0.001                 # Water viscosity Pa s @25 C
rho = 1000                 # Water density kg/m^3
L = 1                      # Pipe's length m
e = 0.000005               # Pipe's rugosity

# This function calculates the velocity of the fluid
def vel(Q, d):
    v = 4*Q/(np.pi * d**2)
    return v
    
# This function calculates the friction coefficient Lambda
def lbd(v, d, mu=mu, rho=rho,e=e):
    reynolds = rho*v*d/mu
    # If RE < 2100 flow is laminar, otherwise it is turbolent
    if reynolds <= 2100:
        lbda = 64/reynolds
    else:
        lbda = 1.325/(np.log(e/(3.7*d)+5.74/(reynolds**0.9)))**2
    return lbda

# This function calculates the pressure drop
def pressureDrop(Q, d, L=L):
    v = vel(Q, d)
    pDrop = 0.5 * lbdv(v,d) * v**2 * d / L
    return pDrop

# This function is used to plot the contour plot
def contour(X, Y, Z, N=20):
    plt.figure()
    CS = plt.contour(X, Y, Z, N)
    plt.clabel(CS, inline=1, fontsize=10)
    plt.title('Pressure drop(m) in 1m of tube')
    plt.xlabel('Q (m^3/h)')
    plt.ylabel('d (cm)')
    plt.show()

# This function is used to make a 3d plot
def dplot(X, Y, Z):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(Q, d, pdrop, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    ax.set_xlabel('Q (m^3/h)')
    ax.set_ylabel('d (cm)')
    ax.set_zlabel('P (m)')
    ax.set_title('Pressure drop (m) in 1m of tube')
    fig.colorbar(surf, shrink=0.5, aspect=5)
    plt.show()

# Vectorize function for calculating lambda. Sooo useful vectorize!!!
lbdv = vectorize(lbd)

# Grid
Q, d = np.meshgrid(Q, d)

# Calculate pressure drop
pdrop = pressureDrop(Q, d)

# Convert pressure drop in meters
pdrop = pdrop * rho / 9.81
# Convert Q in m^3/h and d in cm (for better visualization)
Q = Q * 3600
d = d * 100

# Plot
dplot(Q,d,pdrop)

# Contour plot
contour(Q,d,pdrop)