Centrifugal Compressor Stage Sizer (Schulz / ASME PTC 10)

program centrifugal_compressor_stage_sizing
    implicit none
    integer :: iostat_val, gas_type
    double precision :: P1_bar, T1_C, m_dot_kgs, P2_bar, D2_mm, N_rpm, eta_poly_target
    double precision :: MW_gas, cp_gas, gamma_k, R_gas, Z1, Z2, Z_avg
    double precision :: P1_Pa, P2_Pa, T1_K, PR, n_poly, n_poly_exp, Hp_Jkg, Hp_kJkg
    double precision :: U2_ms, a1_ms, Mu2, head_coeff_mu, T2_K, T2_C, P_shaft_kW, P_shaft_MW
    double precision :: surge_margin_pct
    character(len=32) :: mach_status
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: R_UNIV = 8314.462618d0

    ! Read inputs
    read(*,*,iostat=iostat_val) gas_type         ! 1=Methane/Natural Gas, 2=Air, 3=Nitrogen, 4=Carbon Dioxide, 5=Refrigerant R134a
    read(*,*,iostat=iostat_val) P1_bar           ! Inlet Pressure P1 [bar a] (e.g. 25.0)
    read(*,*,iostat=iostat_val) T1_C             ! Inlet Temperature T1 [deg C] (e.g. 35.0)
    read(*,*,iostat=iostat_val) m_dot_kgs        ! Gas Mass Flow Rate [kg/s] (e.g. 18.5)
    read(*,*,iostat=iostat_val) P2_bar           ! Target Discharge Pressure P2 [bar a] (e.g. 55.0)
    read(*,*,iostat=iostat_val) D2_mm            ! Impeller Outer Tip Diameter D2 [mm] (e.g. 480.0)
    read(*,*,iostat=iostat_val) N_rpm            ! Rotational Speed N [rpm] (e.g. 11500.0)
    read(*,*,iostat=iostat_val) eta_poly_target  ! Polytropic Efficiency [0.70 to 0.88] (e.g. 0.82)

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

    if (P1_bar <= 0.0d0 .or. P2_bar <= P1_bar .or. m_dot_kgs <= 0.0d0 .or. D2_mm <= 0.0d0 .or. N_rpm <= 0.0d0) then
        write(*,*) 'ERROR: Pressures, flow rate, impeller diameter, and speed must be valid and positive.'
        stop
    end if

    ! Gas Properties
    if (gas_type == 1) then ! Natural Gas / Methane
        MW_gas = 16.04d0; gamma_k = 1.31d0; Z1 = 0.96d0; Z2 = 0.93d0
    else if (gas_type == 2) then ! Air
        MW_gas = 28.97d0; gamma_k = 1.40d0; Z1 = 1.00d0; Z2 = 1.00d0
    else if (gas_type == 3) then ! Nitrogen
        MW_gas = 28.01d0; gamma_k = 1.40d0; Z1 = 0.99d0; Z2 = 0.99d0
    else if (gas_type == 4) then ! CO2
        MW_gas = 44.01d0; gamma_k = 1.28d0; Z1 = 0.95d0; Z2 = 0.90d0
    else ! R134a
        MW_gas = 102.03d0; gamma_k = 1.13d0; Z1 = 0.92d0; Z2 = 0.85d0
    end if

    R_gas = R_UNIV / MW_gas
    Z_avg = 0.5d0 * (Z1 + Z2)
    T1_K = T1_C + 273.15d0
    P1_Pa = P1_bar * 1.0d5
    P2_Pa = P2_bar * 1.0d5
    PR = P2_bar / P1_bar

    ! Polytropic Exponent (Schulz / ASME PTC 10)
    n_poly_exp = (gamma_k - 1.0d0) / (gamma_k * max(0.50d0, eta_poly_target))
    n_poly = 1.0d0 / (1.0d0 - n_poly_exp)

    ! Polytropic Head Hp [J/kg]
    Hp_Jkg = Z_avg * R_gas * T1_K * (1.0d0 / n_poly_exp) * ((PR**n_poly_exp) - 1.0d0)
    Hp_kJkg = Hp_Jkg / 1000.0d0

    ! Impeller Kinematics
    U2_ms = PI * (D2_mm / 1000.0d0) * (N_rpm / 60.0d0)
    a1_ms = sqrt(gamma_k * Z1 * R_gas * T1_K)
    Mu2 = U2_ms / a1_ms

    if (Mu2 > 1.25d0) then
        mach_status = 'TRANSONIC / SUPERSONIC TIP'
    else if (Mu2 > 0.85d0) then
        mach_status = 'HIGH SUBSONIC TIP'
    else
        mach_status = 'MODERATE SUBSONIC TIP'
    end if

    ! Head Coefficient (Work coefficient mu_p)
    head_coeff_mu = Hp_Jkg / max(1.0d0, U2_ms**2)

    ! Discharge Temperature T2
    T2_K = T1_K * (PR**n_poly_exp)
    T2_C = T2_K - 273.15d0

    ! Shaft Power Consumption [kW / MW] (assuming 98% mechanical efficiency)
    P_shaft_kW = (m_dot_kgs * Hp_kJkg) / (eta_poly_target * 0.98d0)
    P_shaft_MW = P_shaft_kW / 1000.0d0

    ! Surge Margin Estimation
    surge_margin_pct = 22.0d0 + 5.0d0 * (head_coeff_mu - 0.50d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CENTRIFUGAL COMPRESSOR STAGE SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Pressure Ratio (PR) / Tip Speed = ', PR, ' / ', U2_ms, ' m/s'
    write(*,'(A,F10.2,A,A)')       'Impeller Tip Mach (Mu2)         = ', Mu2, ' | ', trim(mach_status)
    write(*,'(A,F10.3,A,F10.1,A)') 'Work Coeff (μ) / Surge Margin   = ', head_coeff_mu, ' / ', surge_margin_pct, ' %'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'POLYTROPIC HEAD (Hp)         = ', Hp_kJkg, ' kJ/kg'
    write(*,'(A,F10.2,A,F10.3,A)') 'SHAFT POWER (P_shaft)        = ', P_shaft_kW, ' kW (', P_shaft_MW, ' MW)'
    write(*,'(A,F10.1,A)')  'Discharge Gas Temperature T2 = ', T2_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Polytropic Exponent (n)      = ', n_poly, ''
    write(*,'(A,F10.1,A)')  'Speed of Sound at Inlet (a1) = ', a1_ms, ' m/s'
    write(*,'(A)') '============================================================'

end program centrifugal_compressor_stage_sizing