Centrifugal Pump Impeller & Slip Factor (Wiesner / Stodola)

program centrifugal_pump_impeller_vane_slip
    implicit none
    integer :: iostat_val, num_vanes
    double precision :: Q_m3h, H_target_m, N_rpm, rho_liquid_kgm3
    double precision :: D1_eye_mm, D2_tip_mm, b2_width_mm, beta2b_deg, eta_h_target
    double precision :: Q_m3s, D2_m, b2_m, U2_ms, Vr2_ms, Vu2_inf, He_inf_m
    double precision :: sigma_wiesner, sigma_stodola, V_slip_ms, Vu2_actual, He_actual_m
    double precision :: H_actual_m, Ns_metric, P_hyd_kW, P_shaft_kW, eta_overall
    character(len=32) :: pump_regime
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: G_ACC = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Q_m3h             ! Flow Rate Q [m3/h] (e.g. 350.0)
    read(*,*,iostat=iostat_val) H_target_m        ! Head Target H [m] (e.g. 65.0)
    read(*,*,iostat=iostat_val) N_rpm             ! Shaft Speed N [rpm] (e.g. 2950.0)
    read(*,*,iostat=iostat_val) rho_liquid_kgm3   ! Liquid Density [kg/m3] (e.g. 1000.0)
    read(*,*,iostat=iostat_val) D1_eye_mm         ! Eye Suction Diameter D1 [mm] (e.g. 140.0)
    read(*,*,iostat=iostat_val) D2_tip_mm         ! Tip Outer Diameter D2 [mm] (e.g. 260.0)
    read(*,*,iostat=iostat_val) b2_width_mm       ! Exit Width b2 [mm] (e.g. 22.0)
    read(*,*,iostat=iostat_val) beta2b_deg        ! Blade Exit Angle β2b [deg] (e.g. 27.5)
    read(*,*,iostat=iostat_val) num_vanes         ! Number of Vanes Z (e.g. 6)
    read(*,*,iostat=iostat_val) eta_h_target      ! Hydraulic Efficiency [0.75 to 0.90] (e.g. 0.84)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for centrifugal pump impeller sizing.'
        stop
    end if

    if (Q_m3h <= 0.0d0 .or. N_rpm <= 0.0d0 .or. D2_tip_mm <= 0.0d0 .or. num_vanes < 3) then
        write(*,*) 'ERROR: Flow rate, speed, diameter, and vanes count must be positive and valid.'
        stop
    end if

    Q_m3s = Q_m3h / 3600.0d0
    D2_m = D2_tip_mm / 1000.0d0
    b2_m = b2_width_mm / 1000.0d0

    ! Specific speed Metric Ns = N * sqrt(Q_m3s) / (H^0.75)
    Ns_metric = N_rpm * sqrt(Q_m3s) / (max(1.0d0, H_target_m)**0.75d0)

    if (Ns_metric < 30.0d0) then
        pump_regime = 'RADIAL IMPELLER (LOW Ns)'
    else if (Ns_metric < 80.0d0) then
        pump_regime = 'MIXED-FLOW IMPELLER (MEDIUM Ns)'
    else
        pump_regime = 'AXIAL-FLOW PROPELLER (HIGH Ns)'
    end if

    ! Impeller Kinematics
    U2_ms = PI * D2_m * (N_rpm / 60.0d0)
    Vr2_ms = Q_m3s / (PI * D2_m * b2_m * 0.92d0) ! 92% blockage allowance

    ! Infinite Vane Euler Head
    Vu2_inf = U2_ms - (Vr2_ms / tan(beta2b_deg * PI / 180.0d0))
    He_inf_m = (U2_ms * Vu2_inf) / G_ACC

    ! Slip Factor Correlations (Wiesner & Stodola)
    sigma_wiesner = 1.0d0 - (sqrt(sin(beta2b_deg * PI / 180.0d0)) / (dble(num_vanes)**0.70d0))
    V_slip_ms = (PI * U2_ms * sin(beta2b_deg * PI / 180.0d0)) / dble(num_vanes)
    sigma_stodola = 1.0d0 - (V_slip_ms / max(1.0d0, U2_ms))

    ! Real Euler & Actual Head
    Vu2_actual = sigma_wiesner * Vu2_inf
    He_actual_m = (U2_ms * Vu2_actual) / G_ACC
    H_actual_m = eta_h_target * He_actual_m

    ! Power and Efficiencies
    P_hyd_kW = (rho_liquid_kgm3 * G_ACC * Q_m3s * H_actual_m) / 1000.0d0
    eta_overall = eta_h_target * 0.96d0 ! 96% mechanical & volumetric
    P_shaft_kW = P_hyd_kW / eta_overall

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CENTRIFUGAL PUMP IMPELLER & SLIP SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Tip Speed U2 / Exit Radial Vr2= ', U2_ms, ' m/s / ', Vr2_ms, ' m/s'
    write(*,'(A,F10.1,A,A)')       'Metric Specific Speed (Ns)   = ', Ns_metric, ' | ', trim(pump_regime)
    write(*,'(A,F10.3,A,F10.3,A)') 'Wiesner / Stodola Slip Factor= ', sigma_wiesner, ' / ', sigma_stodola, ''
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'ACTUAL TOTAL HEAD DELIVERED  = ', H_actual_m, ' m (Target: ', H_target_m, ' m)'
    write(*,'(A,F10.2,A,F10.2,A)') 'Euler Head (Infinite / Real) = ', He_inf_m, ' m / ', He_actual_m, ' m'
    write(*,'(A,F10.2,A,F10.1,A)') 'SHAFT POWER DEMAND           = ', P_shaft_kW, ' kW (Overall Eff: ', eta_overall*100.0d0, '%)'
    write(*,'(A,F10.2,A,I3,A)')    'Vane Exit Angle / Vanes Count= ', beta2b_deg, ' deg / ', num_vanes, ' vanes'
    write(*,'(A)') '============================================================'

end program centrifugal_pump_impeller_vane_slip