🔬
Solver Purpose & Physical Scope
Fits quadratic centrifugal pump curves from coordinates or coefficients, models pipeline system resistance curves using Colebrook-White friction factors, solves for active operating points across single, series, or parallel configurations, and evaluates suction line cavitation safety margins via dynamic NPSH checks.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 382 times
📄 Source File:
pump_system.f90
📁 calcul/FluidMechanics /
pump_system.f90
! ============================================================================
! ThermoFluidCalc — Pump & System Curve Solver
! Reference: Karassik, Pump Handbook; Cengel & Cimbala, Fluid Mechanics Ch.14
! ============================================================================
program pump_system
implicit none
! Inputs: System parameters
double precision :: delta_z ! Static head [m]
double precision :: D ! Pipe diameter [mm]
double precision :: L ! Pipe length [m]
double precision :: epsilon ! Pipe roughness [mm]
double precision :: sum_K ! Minor loss coefficient sum
double precision :: rho ! Fluid density [kg/m³]
double precision :: mu ! Fluid dynamic viscosity [Pa·s]
! Inputs: Pump parameters
integer :: curve_option ! 1 = 3 Points, 2 = Coefficients, 3 = Preset
double precision :: Q1, H1, Q2, H2, Q3, H3 ! Option 1: 3 Points
double precision :: coeff_a, coeff_b ! Option 2: a - bQ^2
integer :: preset_id ! Option 3: Preset (1=Low, 2=Med, 3=High)
! Inputs: Suction and NPSH parameters
double precision :: delta_z_s ! Suction static head [m]
double precision :: L_s ! Suction pipe length [m]
double precision :: sum_Ks ! Suction minor loss sum
double precision :: Pv ! Vapor pressure [kPa]
double precision :: NPSHr ! NPSH required [m]
double precision :: efficiency_peak ! Peak pump efficiency [%]
! Constants
double precision, parameter :: g = 9.81d0
double precision, parameter :: pi = 3.141592653589793d0
double precision, parameter :: Patm = 101.325d0 ! atmospheric pressure [kPa]
! Solver intermediate variables
double precision :: A, B, C ! Pump curve coefficients H = A*Q^2 + B*Q + C (Q in m³/h)
double precision :: det, det_A, det_B, det_C
double precision :: Q_op1, H_op1 ! Operating point: Single
double precision :: Q_op2, H_op2 ! Operating point: Parallel
double precision :: Q_op3, H_op3 ! Operating point: Series
double precision :: eff_op1, eff_op2, eff_op3 ! Efficiencies at operating points
double precision :: power_op1, power_op2, power_op3 ! Required shaft powers [kW]
double precision :: NPSHa_op1, NPSHa_op2, NPSHa_op3 ! NPSH available
double precision :: runout_flow ! Runout flow (zero static head)
double precision :: h_Ls_op1, h_Ls_op2, h_Ls_op3 ! Suction head losses
! Read input from stdin
read(*,*) delta_z
read(*,*) D
read(*,*) L
read(*,*) epsilon
read(*,*) sum_K
read(*,*) rho
read(*,*) mu
read(*,*) curve_option
read(*,*) Q1, H1
read(*,*) Q2, H2
read(*,*) Q3, H3
read(*,*) coeff_a, coeff_b
read(*,*) preset_id
read(*,*) delta_z_s
read(*,*) L_s
read(*,*) sum_Ks
read(*,*) Pv
read(*,*) NPSHr
read(*,*) efficiency_peak
! Set defaults and safety bounds
if (rho <= 0.0d0) rho = 1000.0d0
if (mu <= 0.0d0) mu = 0.001d0
if (D <= 0.0d0) D = 50.0d0
if (L < 0.0d0) L = 0.0d0
if (epsilon < 0.0d0) epsilon = 0.0d0
if (sum_K < 0.0d0) sum_K = 0.0d0
if (L_s < 0.0d0) L_s = 2.0d0
if (sum_Ks < 0.0d0) sum_Ks = 1.5d0
if (Pv < 0.0d0) Pv = 2.34d0 ! vapor pressure of water at 20°C
if (NPSHr < 0.0d0) NPSHr = 3.0d0
if (efficiency_peak <= 0.0d0) efficiency_peak = 75.0d0
! ── FIT PUMP CURVE H = A*Q^2 + B*Q + C (Q in m³/h) ───────
if (curve_option == 1) then
! 3 points: solve system using Cramer's rule
det = Q1**2 * (Q2 - Q3) + Q2**2 * (Q3 - Q1) + Q3**2 * (Q1 - Q2)
if (abs(det) > 1.0d-6) then
det_A = H1 * (Q2 - Q3) + H2 * (Q3 - Q1) + H3 * (Q1 - Q2)
det_B = Q1**2 * (H2 - H3) + Q2**2 * (H3 - H1) + Q3**2 * (H1 - H2)
det_C = Q1**2 * (Q2*H3 - Q3*H2) + Q2**2 * (Q3*H1 - Q1*H3) + Q3**2 * (Q1*H2 - Q2*H1)
A = det_A / det
B = det_B / det
C = det_C / det
else
! Fallback if points are collinear
A = -0.0005d0
B = -0.05d0
C = H1
end if
elseif (curve_option == 2) then
! Coefficients: H = a - bQ^2
A = -coeff_b
B = 0.0d0
C = coeff_a
else
! Preset options
select case (preset_id)
case (1) ! Low Head
A = -0.001d0
B = -0.05d0
C = 15.0d0
case (2) ! Medium Head
A = -0.0005d0
B = -0.05d0
C = 30.0d0
case (3) ! High Head
A = -0.0004444d0
B = -0.066667d0
C = 60.0d0
case default
A = -0.0005d0
B = -0.05d0
C = 30.0d0
end select
end if
! ── SOLVE OPERATING POINTS (Q in m³/h) ───────────────────
call find_op_bisection(1, delta_z, D, L, epsilon, sum_K, rho, mu, A, B, C, Q_op1)
call find_op_bisection(2, delta_z, D, L, epsilon, sum_K, rho, mu, A, B, C, Q_op2)
call find_op_bisection(3, delta_z, D, L, epsilon, sum_K, rho, mu, A, B, C, Q_op3)
! Compute Heads at operating points
H_op1 = system_head(Q_op1, delta_z, D, L, epsilon, sum_K, rho, mu)
H_op2 = system_head(Q_op2, delta_z, D, L, epsilon, sum_K, rho, mu)
H_op3 = system_head(Q_op3, delta_z, D, L, epsilon, sum_K, rho, mu)
! Fit a parabolic efficiency curve peaking at BEP
! BEP is assumed at 60% of Q_max (where H_pump = 0)
! Q_max is solved from A*Q^2 + B*Q + C = 0
! eta = efficiency_peak * [1 - ( (Q - Q_bep)/Q_bep )^2]
call compute_efficiency_and_power(Q_op1, H_op1, A, B, C, efficiency_peak, rho, eff_op1, power_op1)
! For parallel configuration, each pump operates at Q_op2 / 2
call compute_efficiency_and_power(Q_op2 / 2.0d0, H_op2, A, B, C, efficiency_peak, rho, eff_op2, power_op2)
! Total power is 2 * individual pump power
power_op2 = 2.0d0 * power_op2
! For series configuration, each pump operates at Q_op3 and generates H_op3 / 2
call compute_efficiency_and_power(Q_op3, H_op3 / 2.0d0, A, B, C, efficiency_peak, rho, eff_op3, power_op3)
power_op3 = 2.0d0 * power_op3
! ── NPSH AVAILABLE CHECKS ────────────────────────────────
! h_Ls = dynamic friction loss in suction pipe at flow rate
h_Ls_op1 = suction_head_loss(Q_op1, D, L_s, epsilon, sum_Ks, rho, mu)
h_Ls_op2 = suction_head_loss(Q_op2 / 2.0d0, D, L_s, epsilon, sum_Ks, rho, mu)
h_Ls_op3 = suction_head_loss(Q_op3, D, L_s, epsilon, sum_Ks, rho, mu)
NPSHa_op1 = ((Patm - Pv) * 1000.0d0) / (rho * g) + delta_z_s - h_Ls_op1
NPSHa_op2 = ((Patm - Pv) * 1000.0d0) / (rho * g) + delta_z_s - h_Ls_op2
NPSHa_op3 = ((Patm - Pv) * 1000.0d0) / (rho * g) + delta_z_s - h_Ls_op3
! Runout flow rate (operating point with zero static head delta_z = 0)
call find_op_bisection(1, 0.0d0, D, L, epsilon, sum_K, rho, mu, A, B, C, runout_flow)
! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
write(*, '(A, F14.6)') "Pump A = ", A
write(*, '(A, F14.6)') "Pump B = ", B
write(*, '(A, F14.6)') "Pump C = ", C
write(*, '(A, F14.4)') "Single Q = ", Q_op1
write(*, '(A, F14.4)') "Single H = ", H_op1
write(*, '(A, F14.2)') "Single Eff = ", eff_op1
write(*, '(A, F14.4)') "Single Power = ", power_op1
write(*, '(A, F14.4)') "Single NPSHa = ", NPSHa_op1
write(*, '(A, F14.4)') "Parallel Q = ", Q_op2
write(*, '(A, F14.4)') "Parallel H = ", H_op2
write(*, '(A, F14.2)') "Parallel Eff = ", eff_op2
write(*, '(A, F14.4)') "Parallel Power = ", power_op2
write(*, '(A, F14.4)') "Parallel NPSHa = ", NPSHa_op2
write(*, '(A, F14.4)') "Series Q = ", Q_op3
write(*, '(A, F14.4)') "Series H = ", H_op3
write(*, '(A, F14.2)') "Series Eff = ", eff_op3
write(*, '(A, F14.4)') "Series Power = ", power_op3
write(*, '(A, F14.4)') "Series NPSHa = ", NPSHa_op3
write(*, '(A, F14.4)') "Runout Flow = ", runout_flow
contains
! Colebrook-White friction factor solver
double precision function colebrook_f(Re, eps_over_D)
double precision, intent(in) :: Re, eps_over_D
double precision :: f, f_old, term
integer :: i
if (Re < 2300.0d0) then
if (Re > 0.0d0) then
colebrook_f = 64.0d0 / Re
else
colebrook_f = 0.0d0
end if
return
end if
! Haaland approximation as initial guess
term = (eps_over_D / 3.7d0)**1.11d0 + 6.9d0 / Re
f = 1.0d0 / (-1.8d0 * log10(term))**2
! Successive approximation loops
do i = 1, 5
f_old = f
term = eps_over_D / 3.7d0 + 2.51d0 / (Re * sqrt(f_old))
f = 1.0d0 / (-2.0d0 * log10(term))**2
if (abs(f - f_old) < 1.0d-8) exit
end do
colebrook_f = f
end function colebrook_f
! System head loss curve calculation
double precision function system_head(Q, dz, diam, length, eps, K_sum, r, viscosity)
double precision, intent(in) :: Q, dz, diam, length, eps, K_sum, r, viscosity
double precision :: Q_m3s, diam_m, eps_over_D, V, Re, f
if (Q <= 0.0d0) then
system_head = dz
return
end if
Q_m3s = Q / 3600.0d0
diam_m = diam / 1000.0d0
eps_over_D = (eps / 1000.0d0) / diam_m
V = 4.0d0 * Q_m3s / (pi * diam_m**2)
Re = r * V * diam_m / viscosity
f = colebrook_f(Re, eps_over_D)
system_head = dz + (f * length / diam_m + K_sum) * (V**2 / (2.0d0 * g))
end function system_head
! Suction pipe head loss calculation for NPSH
double precision function suction_head_loss(Q, diam, length, eps, K_sum, r, viscosity)
double precision, intent(in) :: Q, diam, length, eps, K_sum, r, viscosity
double precision :: Q_m3s, diam_m, eps_over_D, V, Re, f
if (Q <= 0.0d0) then
suction_head_loss = 0.0d0
return
end if
Q_m3s = Q / 3600.0d0
diam_m = diam / 1000.0d0
eps_over_D = (eps / 1000.0d0) / diam_m
V = 4.0d0 * Q_m3s / (pi * diam_m**2)
Re = r * V * diam_m / viscosity
f = colebrook_f(Re, eps_over_D)
suction_head_loss = (f * length / diam_m + K_sum) * (V**2 / (2.0d0 * g))
end function suction_head_loss
! Operating point solver using bisection
subroutine find_op_bisection(config, dz, diam, length, eps, K_sum, r, viscosity, c_A, c_B, c_C, Q_sol)
integer, intent(in) :: config
double precision, intent(in) :: dz, diam, length, eps, K_sum, r, viscosity, c_A, c_B, c_C
double precision, intent(out) :: Q_sol
double precision :: low, high, mid, f_low, f_mid, H_p, H_s
double precision :: discriminant, Q_max_pump
integer :: iter
! Determine upper bound of flow sweep Q_max_pump where H_pump = 0
discriminant = c_B**2 - 4.0d0 * c_A * c_C
if (discriminant >= 0.0d0 .and. c_A /= 0.0d0) then
Q_max_pump = (-c_B - sqrt(discriminant)) / (2.0d0 * c_A)
else
Q_max_pump = 1000.0d0
end if
if (Q_max_pump <= 0.0d0) Q_max_pump = 1000.0d0
! Check if pump can overcome static head at Q = 0
H_p = c_C
if (config == 3) H_p = 2.0d0 * c_C
if (H_p < dz) then
Q_sol = 0.0d0
return
end if
low = 0.0d0
high = Q_max_pump * 1.5d0
do iter = 1, 100
mid = (low + high) / 2.0d0
! Pump Head
select case (config)
case (1) ! Single
H_p = c_A * mid**2 + c_B * mid + c_C
case (2) ! Parallel
H_p = c_A * (mid / 2.0d0)**2 + c_B * (mid / 2.0d0) + c_C
case (3) ! Series
H_p = 2.0d0 * (c_A * mid**2 + c_B * mid + c_C)
end select
! System Head
H_s = system_head(mid, dz, diam, length, eps, K_sum, r, viscosity)
f_mid = H_p - H_s
! Check low bound function value
select case (config)
case (1)
H_p = c_C
case (2)
H_p = c_C
case (3)
H_p = 2.0d0 * c_C
end select
f_low = H_p - dz
if (f_mid * f_low < 0.0d0) then
high = mid
else
low = mid
end if
if (abs(high - low) < 1.0d-6) exit
end do
Q_sol = mid
end subroutine find_op_bisection
! Helper to compute efficiency and power at a flow rate
subroutine compute_efficiency_and_power(Q, H, c_A, c_B, c_C, eta_peak, r, eff, power)
double precision, intent(in) :: Q, H, c_A, c_B, c_C, eta_peak, r
double precision, intent(out) :: eff, power
double precision :: discriminant, Q_max, Q_bep, dev
! Determine Q_max of single pump
discriminant = c_B**2 - 4.0d0 * c_A * c_C
if (discriminant >= 0.0d0 .and. c_A /= 0.0d0) then
Q_max = (-c_B - sqrt(discriminant)) / (2.0d0 * c_A)
else
Q_max = 500.0d0
end if
! BEP is typically around 60% of Q_max
Q_bep = Q_max * 0.60d0
if (Q_bep <= 0.0d0) Q_bep = 100.0d0
! Compute efficiency using standard parabolic decay centered at Q_bep
dev = (Q - Q_bep) / Q_bep
eff = eta_peak * (1.0d0 - dev**2)
if (eff < 5.0d0) eff = 5.0d0
if (eff > eta_peak) eff = eta_peak
! Water power: W_w = rho * g * Q * H (converting Q in m³/h to m³/s)
! Shaft power: W_s = W_w / eff
power = (r * g * (Q / 3600.0d0) * H) / (eff / 100.0d0) / 1000.0d0 ! in kW
if (Q <= 0.0d0) then
eff = 0.0d0
power = 0.0d0
end if
end subroutine compute_efficiency_and_power
end program pump_system
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 pump_system.f90 -o pump_system_calc
2. Execution with input.txt redirection:
pump_system_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
10.0 100.0 50.0 0.05 2.0 1000.0 0.001 3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2 2.0 3.0 1.5 2.34 3.0 75.0
Parameter Description:
Static head [m] Pipe diameter [mm] Pipe length [m] Pipe roughness [mm] Minor loss sum Fluid density [kg/m3] Fluid viscosity [Pa-s] Curve option (1=3pts, 2=coeffs, 3=preset) Q1 H1 Q2 H2 Q3 H3 a b Preset ID Suction static head [m] Suction pipe length [m] Suction minor loss sum Vapor pressure [kPa] NPSH required [m] Peak efficiency [%]