Rotary Twin-Screw Compressor Sizer (Vi & Oil Cooling)

program rotary_twin_screw_compressor_sizing
    implicit none
    integer :: iostat_val, gas_select
    double precision :: P1_bar, T1_C, P2_bar, V_disp_m3h, N_male_rpm, Vi_built_in, oil_flow_lpm
    double precision :: gamma_k, MW_gas, R_gas, cp_gas, rho1_kgm3, m_dot_gas_kgs
    double precision :: Pi_internal, P_internal_bar, W_ideal_kW, W_ind_kW, P_shaft_kW
    double precision :: Q_oil_cooling_kW, T2_oil_C, T2_dry_C, eta_vol_pct, eta_ad_pct
    character(len=32) :: comp_regime
    double precision, parameter :: RHO_OIL = 870.0d0
    double precision, parameter :: CP_OIL = 2.05d0 ! kJ/kg.K

    ! Read inputs
    read(*,*,iostat=iostat_val) gas_select    ! 1=Air, 2=Ammonia NH3, 3=R134a, 4=Propane R290, 5=Biogas/CH4
    read(*,*,iostat=iostat_val) P1_bar        ! Suction Pressure P1 [bar a] (e.g. 1.013)
    read(*,*,iostat=iostat_val) T1_C          ! Suction Temperature T1 [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) P2_bar        ! Discharge Pressure P2 [bar a] (e.g. 8.5)
    read(*,*,iostat=iostat_val) V_disp_m3h    ! Swept Displacement V_disp [m3/h] (e.g. 650.0)
    read(*,*,iostat=iostat_val) N_male_rpm    ! Male Rotor Speed [rpm] (e.g. 2950.0)
    read(*,*,iostat=iostat_val) Vi_built_in   ! Built-in Volume Ratio Vi (e.g. 3.2)
    read(*,*,iostat=iostat_val) oil_flow_lpm  ! Injected Oil Flow Rate [L/min] (e.g. 45.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for twin-screw compressor calculation.'
        stop
    end if

    if (P1_bar <= 0.0d0 .or. P2_bar <= P1_bar .or. V_disp_m3h <= 0.0d0 .or. Vi_built_in <= 1.0d0) then
        write(*,*) 'ERROR: Pressures, flow, and volume ratio must be valid with P2 > P1 and Vi > 1.'
        stop
    end if

    ! Gas Properties
    if (gas_select == 1) then ! Air
        MW_gas = 28.97d0; gamma_k = 1.40d0; cp_gas = 1.005d0
    else if (gas_select == 2) then ! NH3 Ammonia
        MW_gas = 17.03d0; gamma_k = 1.31d0; cp_gas = 2.150d0
    else if (gas_select == 3) then ! R134a
        MW_gas = 102.03d0; gamma_k = 1.13d0; cp_gas = 0.850d0
    else if (gas_select == 4) then ! Propane R290
        MW_gas = 44.10d0; gamma_k = 1.14d0; cp_gas = 1.670d0
    else ! Biogas
        MW_gas = 22.50d0; gamma_k = 1.30d0; cp_gas = 1.450d0
    end if

    R_gas = 8314.46d0 / MW_gas
    rho1_kgm3 = (P1_bar * 1.0d5) / (R_gas * (T1_C + 273.15d0))

    ! Volumetric efficiency estimate based on pressure ratio PR
    eta_vol_pct = 94.0d0 - 1.2d0 * (P2_bar / P1_bar)
    eta_vol_pct = max(75.0d0, min(95.0d0, eta_vol_pct))

    m_dot_gas_kgs = (V_disp_m3h / 3600.0d0) * (eta_vol_pct / 100.0d0) * rho1_kgm3

    ! Built-in pressure ratio Pi
    Pi_internal = Vi_built_in**gamma_k
    P_internal_bar = P1_bar * Pi_internal

    if (P2_bar > P_internal_bar * 1.05d0) then
        comp_regime = 'UNDER-COMPRESSION (P2 > Pi*P1)'
    else if (P2_bar < P_internal_bar * 0.95d0) then
        comp_regime = 'OVER-COMPRESSION (P2 < Pi*P1)'
    else
        comp_regime = 'PERFECT MATCHED EXPANSION'
    end if

    ! Indicator Work [kW]
    W_ind_kW = (m_dot_gas_kgs * R_gas * (T1_C + 273.15d0) / 1000.0d0) * &
               ((gamma_k / (gamma_k - 1.0d0)) * (Vi_built_in**(gamma_k - 1.0d0) - 1.0d0) + &
                ((P2_bar - P_internal_bar) / P1_bar) / Vi_built_in)

    P_shaft_kW = W_ind_kW / 0.92d0 ! 92% mechanical & friction efficiency

    ! Adiabatic ideal power & efficiency
    W_ideal_kW = (m_dot_gas_kgs * (gamma_k / (gamma_k - 1.0d0)) * R_gas * (T1_C + 273.15d0) / 1000.0d0) * &
                 ((P2_bar / P1_bar)**((gamma_k - 1.0d0) / gamma_k) - 1.0d0)
    eta_ad_pct = (W_ideal_kW / max(1.0d0, P_shaft_kW)) * 100.0d0
    eta_ad_pct = min(86.0d0, max(60.0d0, eta_ad_pct))

    ! Discharge Temperatures (Dry vs Oil Injected)
    T2_dry_C = (T1_C + 273.15d0) * ((P2_bar / P1_bar)**((gamma_k - 1.0d0) / gamma_k)) - 273.15d0
    
    ! Oil cooling balance: m_oil * cp_oil * dT_oil = W_friction + compression heat
    Q_oil_cooling_kW = P_shaft_kW * 0.82d0
    T2_oil_C = T1_C + (Q_oil_cooling_kW / max(0.1d0, (oil_flow_lpm * RHO_OIL / 60000.0d0) * CP_OIL + m_dot_gas_kgs * cp_gas))
    T2_oil_C = min(115.0d0, max(50.0d0, T2_oil_C))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ROTARY TWIN-SCREW COMPRESSOR SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Pressure Ratio / Built-in Pi = ', (P2_bar/P1_bar), ' / ', Pi_internal, ''
    write(*,'(A,A)')               'Compression Matching Regime  = ', trim(comp_regime)
    write(*,'(A,F10.1,A,F10.1,A)') 'Discharge Temp (Oil / Dry)   = ', T2_oil_C, ' C / ', T2_dry_C, ' C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'SHAFT POWER CONSUMPTION      = ', P_shaft_kW, ' kW'
    write(*,'(A,F10.2,A,F10.2,A)') 'Adiabatic / Volumetric Eff.  = ', eta_ad_pct, ' % / ', eta_vol_pct, ' %'
    write(*,'(A,F10.3,A,F10.1,A)') 'Mass Flow / Oil Cooling Duty = ', m_dot_gas_kgs, ' kg/s / ', Q_oil_cooling_kW, ' kWth'
    write(*,'(A,F10.2,A,F10.1,A)') 'Built-in Vi / Swept Volume   = ', Vi_built_in, ' / ', V_disp_m3h, ' m3/h'
    write(*,'(A)') '============================================================'

end program rotary_twin_screw_compressor_sizing