Gas Turbine Combustor Sizer (Lefebvre Model)

program gas_turbine_combustor_sizing_lefebvre
    implicit none
    integer :: iostat_val, fuel_type, num_cans
    double precision :: m_dot_air_kgs, P3_bar, T3_K, T4_TIT_K, dp_loss_pct, length_diameter_ratio
    double precision :: LHV_MJkg, FAR_stoich, FAR_actual, m_dot_fuel_kgs, Q_thermal_MW
    double precision :: cp_air, cp_gas, T_pz_K, V_comb_m3, D_liner_m, L_comb_m
    double precision :: m_primary_kgs, m_secondary_kgs, m_dilution_kgs
    double precision :: theta_loading, eta_comb_pct, pattern_factor_otdf
    character(len=32) :: flame_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fuel_type             ! 1=Natural Gas / Methane, 2=Kerosene Jet-A, 3=Diesel #2, 4=Hydrogen H2
    read(*,*,iostat=iostat_val) m_dot_air_kgs         ! Combustor Total Air Mass Flow [kg/s] (e.g. 60.0)
    read(*,*,iostat=iostat_val) P3_bar                ! Compressor Exit Pressure P3 [bar a] (e.g. 16.0)
    read(*,*,iostat=iostat_val) T3_K                  ! Compressor Exit Temperature T3 [K] (e.g. 680.0)
    read(*,*,iostat=iostat_val) T4_TIT_K              ! Turbine Inlet Temperature TIT / T4 [K] (e.g. 1450.0)
    read(*,*,iostat=iostat_val) num_cans              ! Number of Combustor Cans (1=Annular, 8, 12, 14=Can-Annular) (e.g. 1)
    read(*,*,iostat=iostat_val) dp_loss_pct           ! Combustor Pressure Loss ΔP/P3 [%] (e.g. 4.0)
    read(*,*,iostat=iostat_val) length_diameter_ratio ! Combustor L/D Ratio (e.g. 2.8)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for gas turbine combustor calculation.'
        stop
    end if

    if (m_dot_air_kgs <= 0.0d0 .or. P3_bar <= 0.0d0 .or. T3_K <= 0.0d0 .or. T4_TIT_K <= T3_K) then
        write(*,*) 'ERROR: Flow rate, pressures, and temperatures must be positive with T4 > T3.'
        stop
    end if

    ! Fuel Properties
    if (fuel_type == 1) then ! Natural Gas
        LHV_MJkg = 50.0d0; FAR_stoich = 0.058d0
    else if (fuel_type == 2) then ! Kerosene Jet-A
        LHV_MJkg = 43.1d0; FAR_stoich = 0.068d0
    else if (fuel_type == 3) then ! Diesel #2
        LHV_MJkg = 42.5d0; FAR_stoich = 0.069d0
    else ! Hydrogen H2
        LHV_MJkg = 120.0d0; FAR_stoich = 0.029d0
    end if

    cp_air = 1060.0d0
    cp_gas = 1150.0d0

    ! Fuel-Air Ratio & Fuel Flow
    FAR_actual = (cp_gas * T4_TIT_K - cp_air * T3_K) / (LHV_MJkg * 1.0d6 * 0.995d0 - cp_gas * T4_TIT_K)
    m_dot_fuel_kgs = m_dot_air_kgs * FAR_actual
    Q_thermal_MW = (m_dot_fuel_kgs * LHV_MJkg)

    ! Zone Air Splits
    m_primary_kgs = m_dot_air_kgs * 0.28d0    ! 28% in primary zone
    m_secondary_kgs = m_dot_air_kgs * 0.22d0  ! 22% in secondary zone
    m_dilution_kgs = m_dot_air_kgs * 0.50d0   ! 50% in dilution zone

    ! Primary flame zone temperature
    T_pz_K = min(2450.0d0, T3_K + 1650.0d0)

    ! Lefebvre Volume & Dimensions
    V_comb_m3 = (m_dot_air_kgs * 287.0d0 * T3_K) / (P3_bar * 1.0d5 * 25.0d0) ! Residence time ~ 5-8 ms
    D_liner_m = ((4.0d0 * V_comb_m3) / (PI * length_diameter_ratio * dble(max(1, num_cans))))**(1.0d0/3.0d0)
    L_comb_m = D_liner_m * length_diameter_ratio

    ! Lefebvre Loading Parameter & Efficiency
    theta_loading = ((P3_bar * 1.0d5 / 1.0d5)**1.75d0 * V_comb_m3 * exp(T3_K / 300.0d0)) / m_dot_air_kgs
    eta_comb_pct = (1.0d0 - exp(-0.00015d0 * theta_loading)) * 100.0d0
    eta_comb_pct = min(99.95d0, max(95.0d0, eta_comb_pct))

    ! Pattern Factor (OTDF)
    pattern_factor_otdf = 0.22d0 - 0.05d0 * (L_comb_m / max(0.2d0, D_liner_m) - 2.0d0)
    pattern_factor_otdf = max(0.12d0, min(0.30d0, pattern_factor_otdf))

    if (fuel_type == 4) then
        flame_status = 'ZERO-CARBON HYDROGEN FLAME'
    else if (T_pz_K > 2200.0d0) then
        flame_status = 'HIGH THERMAL NOX ZONE'
    else
        flame_status = 'LEAN PREMIXED LOW-NOX'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — GAS TURBINE COMBUSTOR (LEFEBVRE SIZER)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Thermal Duty (Q_th) / Fuel   = ', Q_thermal_MW, ' MWth / ', m_dot_fuel_kgs, ' kg/s'
    write(*,'(A,F10.4,A,F10.2,A)') 'Fuel-to-Air Ratio (FAR)      = ', FAR_actual, ' (Equivalence phi = ', FAR_actual/FAR_stoich, ')'
    write(*,'(A,F10.1,A,A)')       'Primary Flame Temperature    = ', T_pz_K, ' K | ', trim(flame_status)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Combustor Volume (V_comb)    = ', V_comb_m3, ' m3'
    write(*,'(A,F10.1,A,F10.1,A)') 'Liner Diameter / Length      = ', D_liner_m*1000.0d0, ' mm / ', L_comb_m*1000.0d0, ' mm'
    write(*,'(A,F10.2,A,F10.2,A)') 'Combustion Efficiency / OTDF = ', eta_comb_pct, ' % / ', pattern_factor_otdf, ''
    write(*,'(A,F10.1,A,F10.1,A)') 'Primary / Dilution Air Flow  = ', m_primary_kgs, ' kg/s / ', m_dilution_kgs, ' kg/s'
    write(*,'(A,F10.2,A,F10.1,A)') 'Combustor Pressure Drop (ΔP) = ', (P3_bar * dp_loss_pct / 100.0d0), ' bar (', dp_loss_pct, '%)'
    write(*,'(A)') '============================================================'

end program gas_turbine_combustor_sizing_lefebvre