pythonnumpyscipypde

Solve waterhammer PDE in numpy/scipy

I’ve been trying to solve the water hammer PDE’s from the Maple example linked below in python (numpy/scipy). I’m getting very unstable results. Can anyone see my mistake? Guessing something is wrong with the boundary conditions.

https://www.maplesoft.com/support/help/view.aspx?path=applications/WaterHammer

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

## Parameters
Dia = 0.1
V = 14.19058741 # Stead state
p = 1000 # Liquid density
u = 0.001 # Viscosity
L = 25 
e = 0.0001 # Roughness
Psource = 0.510E6
thick = 0.001
E= 7010*10**9
K=20010E6
Vsteady= 14.19058741
Ks = 1/((1/K)+(Dia/E*thick))

# Darcy-Weisbach
def Friction(V):
    Rey = ((Dia*V*p)/u)
    fL = 64/Rey
    fT = 1/((1.8*np.log10((6.9/Rey) + (e/(3.7*Dia))**1.11))**2)
    
    if Rey >= 0 and Rey < 2000:
        return fL
    
    if Rey >= 2000 and Rey<4000:
        return fL + ((fT-fL)*(Rey-2000))/(4000-2000)
     
    if Rey >= 4000:
        return fT
    
    return 0

def model(D, t):
    V = D[:N]
    P = D[N:]
    
    dVdt = np.zeros(N)
    for i in range(1, len(dVdt)-1):
        dVdt[i] = -(1/p)*((P[i+1]-P[i-1])/2*dx)-((Friction(np.abs(V[i]))*(np.abs(V[i])**2))/(2*Dia))
    
    dPdt = np.zeros(N)
    for i in range(1, len(dPdt)-1):
        dPdt[i] = -((V[i+1]-V[i-1])/(2*dx))*Ks
        
    if t < 2:
        dVdt[29] = 0
    else:
        dVdt[29] = -1
     
    dPdt[29] = 0
    dVdt[0] = dVdt[1]
           
    return np.append(dVdt,dPdt)


N = 30
x = np.linspace(0, L, N)
dx = x[1] - x[0]

## Initial conditions
Vi_0 = np.ones(N)*Vsteady
Pi_0 = np.arange(N)
for i in Pi_0:
    Pi_0[i] = Psource - (i*dx/L)*Psource

# initial condition
y0 = np.append(Vi_0, Pi_0)

# time points
t = np.linspace(0,3,10000)

# solve ODE
y = odeint(model,y0,t)

Vr = y[:,0:N]
Pr = y[:,N:]

plt.plot(t,Pr[:,5])

Pressure plot

Solution

  • I have just come across this posting and have modified your code to obtain exactly the same results in python as from the Maple code you cite. The following points are worth making:

    The python code:

    import numpy as np
import scipy as sc
from scipy import integrate
# from scipy.integrate import odeint
import matplotlib.pyplot as plt
import time

t_start = time.time()

## Parameters
p = 1000                       # Liquid density
K = 200*1E6                    # Bulk Modulus of Liquid
u = 0.001                      # Viscosity
Dia = 0.1                      # Pipe Diameter
thick = 0.001                  # Pipe Thickness
e = 0.0001                     # Roughness
L = 25                         # Length of Pipe
E = 70*1E9                     # Young's Modulus
A = (np.pi*Dia**2)/4           # Cross-sectional Area
Psource = 0.5*1E6              # Pressure at Start of Pipeline
Ks = 1/((1/K)+(Dia/(E*thick))) # Effective Modulus of system

Vsteady = 14.19058741          # Steady State (from Maple)
a0 = np.sqrt(K/p)              # Speed of Sound

# Darcy-Weisbach
def Friction(V):
    Rey = ((Dia*V*p)/u)
    fL = 64/Rey
    fT = 1/((1.8*np.log10((6.9/Rey) + (e/(3.7*Dia))**1.11))**2)
    
    if Rey > 0 and Rey < 2000:
        return fL    
    elif Rey >= 2000 and Rey<4000:
        return fL + ((fT-fL)*(Rey-2000))/(4000-2000)
    elif Rey >= 4000:
        return fT    
    else: 
      return 0

def model(t, D):
    V = np.zeros(N+2)    # Allocate Nodes with ghost cells   
    P = np.zeros(N+2)
    
    V[1:N+1] = D[0:N]    # Populate Node Data
    P[1:N+1] = D[N:2*N]
    
    dVdt = np.zeros(N+2) # Allocate Derivative vars with ghost cells
    dPdt = np.zeros(N+2)
    if t >= t_trip:
      # Set ghost cells equal to BCs, as Maple
      V[0] = V[1]
      P[0] = Psource
      V[N+1] = Vsteady*np.exp(-70*(t-t_trip))
      P[N+1] = 0
      
      for i in range(1, N+1):
          dVdt[i] = -(1/p)*(P[i+1]-P[i-1])/(2*dx)- \
                     ((Friction(np.abs(V[i]))*(V[i]*np.abs(V[i])))/(2*Dia))   
      for i in range(1, N+1):
          dPdt[i] = -((V[i+1]-V[i-1])/(2*dx))*Ks

    return np.append(dVdt[1:N+1],dPdt[1:N+1])

N = 30
x = np.linspace(0, L, N)
dx = x[1] - x[0]

## Initial conditions
Vi_0 = np.ones(N)*Vsteady
Pi_0 = np.zeros(N)
for i in range(N):
    Pi_0[i] = Psource - (i*dx/L)*Psource
y0 = np.append(Vi_0, Pi_0)

# time points
t0 = 0.0
tf = 3.0
t_trip = 2

# solve ODE
method = "DOP853" #"LSODA","DOP853","BDF","RK23","RK45,Radau""
tol = 1.e-12
out = sc.integrate.solve_ivp(fun=model,
      t_span=[t0,tf], y0=y0,
      method=method, first_step=1e-6, max_step=0.0001,
      rtol=tol, atol=tol)

print(r"Integration Success =", out.success)
print(out.message)

dP_max = p*a0*Vsteady
print(f"Maximum Theoretical Pressure = {dP_max/10**6:.2f} [MPa]")

t_end = time.time()
total_time = t_end-t_start
print(f"Elapsed Time = {total_time:.2f}")

t = out.t
sol = out.y

Vr = sol[0:N,:]
Pr = sol[N:2*N,:]

plotNode = N-1
fig = plt.figure(figsize=(10,6))
ax1, ax2 = fig.subplots(2)
fig.suptitle(f'Water Hammer Plots at Node: {plotNode}',
             fontsize=20)
ax1.plot(t,Pr[plotNode,:])
ax2.plot(t,Vr[plotNode,:])

ax1.set_xlabel('time [s]', fontsize=16)
ax1.set_ylabel('Pr [Pa]', fontsize=16)
ax2.set_xlabel('time [s]', fontsize=16) 
ax2.set_ylabel('Vr [m/s]',fontsize=16)
ax1.grid()
ax2.grid()
fig.tight_layout()

    The output plots generated by the python code are included. enter image description here