๐Ÿ” Search & Filter All 178+ Engineering Calculators โ†’

๐Ÿ’ป Fortran Source Code Library

We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.

Pump & System Curve Calculator

Core Numerical Engine in Fortran 90 โ€ข 34 total downloads

pump_system.f90 
! =========================================================================
! Source File: 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

Solver Description

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.

Key Numerical Methods & Architecture

  • Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
  • Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
  • Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

๐Ÿ› ๏ธ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 pump_system.f90 -o pump_system_calc

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

pump_system_calc < input.txt

๐Ÿ“ฅ Downloads & Local Files

Preview of the required input file (input.txt):

! Static head [m]
10.0
! Pipe diameter [mm]
100.0
! Pipe length [m]
50.0
! Pipe roughness [mm]
0.05
! Minor loss sum
2.0
! Fluid density [kg/m3]
1000.0
! Fluid viscosity [Pa-s]
0.001
! Curve option (1=3pts, 2=coeffs, 3=preset)
3
! Q1
0.0 0.0
! H1
0.0 0.0
! Q2
0.0 0.0
! H2
0.0 0.0
! Q3
2
! H3
2.0
! a
3.0
! b
1.5
! Preset ID
2.34
! Suction static head [m]
3.0
! Suction pipe length [m]
75.0
๐Ÿ“ฆ Download Offline Bundle (ZIP) 9 dl ๐Ÿ“„ Download input.txt 7 dl ๐Ÿ“ฅ Download pump_system.f90 18 dl