fork download
copy 
  1. import numpy as np
  2. from scipy.interpolate import interp1d
  3.  
  4. def calculate_reynolds(Q, D):
  5.     """Calculate Reynolds number"""
  6.     A = np.pi * (D/2)**2
  7.     v = Q / A
  8.     return (v * D) / (1.005e-6)
  9.  
  10. def friction_factor(Re):
  11.     """Calculate friction factor for smooth pipe using Blasius equation"""
  12.     if Re < 2300:
  13.         return 64/Re
  14.     else:
  15.         return 0.316 * Re**(-0.25)  # Blasius equation for turbulent flow
  16.  
  17. def head_loss(Q):
  18.     """Calculate total head loss for given flow rate"""
  19.     # Convert Q from l/s to m³/s
  20.     Q = Q / 1000
  21.  
  22.     # Pipe dimensions
  23.     D_large = 0.028  # m
  24.     D_small = 0.0136  # m
  25.  
  26.     # Lengths
  27.     L_main = 7.5  # m
  28.     L_large_fig2 = (0.56 + 0.24 + 0.09 + 0.095 + 0.8)  # m
  29.     L_small = (0.4 + 0.4)  # m
  30.  
  31.     # Calculate velocities
  32.     A_large = np.pi * (D_large/2)**2
  33.     A_small = np.pi * (D_small/2)**2
  34.     v_large = Q / A_large
  35.     v_small = Q / A_small
  36.  
  37.     # Reynolds numbers
  38.     Re_large = calculate_reynolds(Q, D_large)
  39.     Re_small = calculate_reynolds(Q, D_small)
  40.  
  41.     # Friction factors
  42.     f_large = friction_factor(Re_large)
  43.     f_small = friction_factor(Re_small)
  44.  
  45.     # Major losses
  46.     h_major_main = f_large * (L_main + L_large_fig2) * v_large**2 / (2 * 9.81 * D_large)
  47.     h_major_small = f_small * L_small * v_small**2 / (2 * 9.81 * D_small)
  48.  
  49.     # Minor losses
  50.     # Sum of K values from Fig 2.1 = 2
  51.     # Additional K values from Fig 2.2: bends and area changes
  52.     K_bends = 4 * 0.5  # Assuming K ≈ 0.5 for 90° bends with r/D = 1
  53.     K_contractions = 0.5  # Approximate value for sudden contraction
  54.     K_expansions = 1.0  # Approximate value for sudden expansion
  55.  
  56.     h_minor_large = (2 + K_bends) * v_large**2 / (2 * 9.81)
  57.     h_minor_small = (K_contractions + K_expansions) * v_small**2 / (2 * 9.81)
  58.  
  59.     return h_major_main + h_major_small + h_minor_large + h_minor_small
  60.  
  61. # Pump curve data (Q in l/s, H in m)
  62. Q_pump = np.array([0, 1, 2, 3, 4])
  63. H_pump = np.array([6.8, 6.2, 5.4, 4.2, 2.8])
  64.  
  65. # Create interpolation function for pump curve
  66. pump_curve = interp1d(Q_pump, H_pump, kind='linear')
  67.  
  68. # Find intersection point
  69. Q_range = np.linspace(0, 4, 1000)
  70. H_system = [head_loss(q) for q in Q_range]
  71. H_pump_range = pump_curve(Q_range)
  72.  
  73. # Find where difference is minimal
  74. intersection_idx = np.argmin(np.abs(H_system - H_pump_range))
  75. Q_operating = Q_range[intersection_idx]
  76. H_operating = H_system[intersection_idx]
  77.  
  78. print(f"Operating point:")
  79. print(f"Q = {Q_operating:.2f} l/s")
  80. print(f"H = {H_operating:.2f} m")
Success #stdin #stdout #stderr 0.25s 46100KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Operating point:
Q = 0.00 l/s
H = nan m
stderr
copy 
./prog.py:13: RuntimeWarning: divide by zero encountered in double_scalars
./prog.py:46: RuntimeWarning: invalid value encountered in double_scalars
./prog.py:47: RuntimeWarning: invalid value encountered in double_scalars
https://ideone.com/fQOv7w
language:
Python 3 (python  3.12)
created:
10 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!