💻 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.

Adiabatic Flame Temperature (Tools)

Core Numerical Engine in Fortran 90 • 28 total downloads

combustion_aft.f90
! =========================================================================
! Source File: combustion_aft.f90
! =========================================================================

program combustion_aft
    implicit none
    integer :: fuel_type, iostat_val, i, n_sweep, iter
    double precision :: excess_air, T_inlet, P_kPa
    double precision :: C_at, H_at, O_at, M_fuel, LHV
    double precision :: a_stoich, air_stoich, AFR_stoich, AFR_actual
    double precision :: actual_O2, n_CO2, n_H2O, n_N2, n_O2_ex, n_total
    double precision :: y_CO2, y_H2O, y_N2, y_O2
    double precision :: T_ad, T_low, T_high, T_mid
    double precision :: H_react, H_prod, hf_fuel
    double precision :: cp_CO2, cp_H2O, cp_N2, cp_O2
    double precision :: hf_CO2, hf_H2O
    double precision :: P_H2O, T_dew
    double precision :: exc_sw, T_ad_sw, AFR_sw

    ! Average cp in J/(mol.K)
    cp_CO2 = 54.0d0
    cp_H2O = 41.0d0
    cp_N2  = 32.0d0
    cp_O2  = 34.0d0

    ! Enthalpies of formation in kJ/mol
    hf_CO2 = -393.5d0
    hf_H2O = -241.8d0

    ! ── Read inputs ─────────────────────────────────────────────
    read(*,*,iostat=iostat_val) fuel_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid fuel type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) excess_air
    read(*,*,iostat=iostat_val) T_inlet
    read(*,*,iostat=iostat_val) P_kPa
    read(*,*,iostat=iostat_val) C_at
    read(*,*,iostat=iostat_val) H_at
    read(*,*,iostat=iostat_val) O_at
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all inputs.'
        stop
    end if

    ! ── Input validation ────────────────────────────────────────
    if (excess_air < 0.0d0) excess_air = 0.0d0
    if (T_inlet <= 0.0d0) T_inlet = 298.15d0
    if (P_kPa <= 0.0d0) P_kPa = 101.325d0

    ! ── Fuel properties ─────────────────────────────────────────
    if (fuel_type == 1) then
        C_at = 1.0d0; H_at = 4.0d0; O_at = 0.0d0
        M_fuel = 16.04d0; LHV = 50050.0d0
    else if (fuel_type == 2) then
        C_at = 3.0d0; H_at = 8.0d0; O_at = 0.0d0
        M_fuel = 44.10d0; LHV = 46357.0d0
    else if (fuel_type == 3) then
        C_at = 8.0d0; H_at = 18.0d0; O_at = 0.0d0
        M_fuel = 114.23d0; LHV = 44427.0d0
    else if (fuel_type == 4) then
        C_at = 0.0d0; H_at = 2.0d0; O_at = 0.0d0
        M_fuel = 2.016d0; LHV = 119960.0d0
    else
        ! Custom fuel
        if (C_at < 0.0d0) C_at = 1.0d0
        if (H_at < 0.0d0) H_at = 4.0d0
        if (O_at < 0.0d0) O_at = 0.0d0
        M_fuel = C_at*12.011d0 + H_at*1.008d0 + O_at*16.0d0
        if (M_fuel < 1.0d0) M_fuel = 16.04d0
        ! Estimate LHV from Dulong-like formula (kJ/kg)
        LHV = 33900.0d0*C_at*12.011d0/M_fuel + &
              142300.0d0*H_at*1.008d0/M_fuel - &
              9090.0d0*O_at*16.0d0/M_fuel
        if (LHV < 1000.0d0) LHV = 50000.0d0
    end if

    ! ── Stoichiometry ───────────────────────────────────────────
    ! CxHyOz + a(O2 + 3.76 N2) -> C CO2 + H/2 H2O + 3.76a N2
    ! a_stoich = C + H/4 - O_fuel/2
    a_stoich = C_at + H_at / 4.0d0 - O_at / 2.0d0
    if (a_stoich < 0.0d0) a_stoich = 0.001d0

    air_stoich = a_stoich / 0.21d0    ! moles air per mole fuel
    AFR_stoich = air_stoich * 28.97d0 / M_fuel   ! mass basis
    AFR_actual = AFR_stoich * (1.0d0 + excess_air / 100.0d0)

    ! Actual O2 supplied
    actual_O2 = a_stoich * (1.0d0 + excess_air / 100.0d0)

    ! Product moles per mole fuel
    n_CO2 = C_at
    n_H2O = H_at / 2.0d0
    n_N2  = 3.76d0 * actual_O2
    n_O2_ex = a_stoich * excess_air / 100.0d0
    n_total = n_CO2 + n_H2O + n_N2 + n_O2_ex

    ! Mole fractions
    y_CO2 = n_CO2 / max(n_total, 1.0d-10)
    y_H2O = n_H2O / max(n_total, 1.0d-10)
    y_N2  = n_N2  / max(n_total, 1.0d-10)
    y_O2  = n_O2_ex / max(n_total, 1.0d-10)

    ! ── Adiabatic flame temperature (bisection) ────────────────
    ! Enthalpy of formation of fuel from LHV:
    ! LHV = -hf_fuel + C*(-hf_CO2) + H/2*(-hf_H2O_g) [per mole fuel, roughly]
    ! Actually: hf_fuel = C*hf_CO2 + H/2*hf_H2O - (-LHV*M_fuel/1000)
    hf_fuel = n_CO2 * hf_CO2 + n_H2O * hf_H2O + LHV * M_fuel / 1000.0d0

    ! Energy balance: H_react = H_prod(T_ad)
    ! H_react = hf_fuel + actual_O2*(0) + n_N2_react*(0) + cp contributions at T_inlet
    ! H_react preheat = [actual_O2*cp_O2 + 3.76*actual_O2*cp_N2] * (T_inlet-298.15) / 1000
    H_react = hf_fuel + (actual_O2 * cp_O2 + 3.76d0*actual_O2 * cp_N2) * &
              (T_inlet - 298.15d0) / 1000.0d0

    ! H_prod(T) = n_CO2*hf_CO2 + n_H2O*hf_H2O + [n_CO2*cp_CO2 + n_H2O*cp_H2O + n_N2*cp_N2 + n_O2_ex*cp_O2]*(T-298.15)/1000
    ! Solve: H_react = H_prod(T_ad)
    ! => hf_fuel + preheat = n_CO2*hf_CO2 + n_H2O*hf_H2O + cp_prod_total*(T_ad-298.15)/1000
    ! => T_ad = 298.15 + 1000 * (hf_fuel + preheat - n_CO2*hf_CO2 - n_H2O*hf_H2O) / cp_prod_total

    T_low = 300.0d0
    T_high = 5000.0d0
    do iter = 1, 100
        T_mid = 0.5d0 * (T_low + T_high)
        H_prod = n_CO2*hf_CO2 + n_H2O*hf_H2O + &
                 (n_CO2*cp_CO2 + n_H2O*cp_H2O + n_N2*cp_N2 + n_O2_ex*cp_O2) * &
                 (T_mid - 298.15d0) / 1000.0d0
        if (H_prod < H_react) then
            T_low = T_mid
        else
            T_high = T_mid
        end if
        if (abs(T_high - T_low) < 0.1d0) exit
    end do
    T_ad = 0.5d0 * (T_low + T_high)

    ! ── Dew point ───────────────────────────────────────────────
    P_H2O = y_H2O * P_kPa   ! partial pressure of H2O in kPa
    ! Simplified Antoine for water: T_dew(K) ~ 273.15 + 100 * (P_H2O/101.325)^0.25
    ! More accurate simplified: T_dew = 273.15 + 45.0*log(P_H2O) + 1.0 (rough fit for low P)
    if (P_H2O > 0.0d0) then
        T_dew = 273.15d0 + 45.0d0 * log(max(P_H2O, 0.001d0)) + 7.0d0
    else
        T_dew = 273.15d0
    end if
    if (T_dew < 273.15d0) T_dew = 273.15d0

    ! ── Output ──────────────────────────────────────────────────
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   COMBUSTION & ADIABATIC FLAME TEMPERATURE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Fuel Type                 = ', fuel_type
    write(*,'(A,F10.2,A)')  '  Fuel Formula              : C', C_at, ''
    write(*,'(A,F10.2)')    '                              H', H_at
    if (O_at > 0.0d0) &
    write(*,'(A,F10.2)')    '                              O', O_at
    write(*,'(A,F10.2,A)')  '  Molar Mass M_fuel         = ', M_fuel, ' g/mol'
    write(*,'(A,F12.2,A)')  '  LHV                       = ', LHV, ' kJ/kg'
    write(*,'(A,F10.2,A)')  '  Excess Air                = ', excess_air, ' percent'
    write(*,'(A,F12.2,A)')  '  Inlet Temperature         = ', T_inlet, ' K'
    write(*,'(A,F12.4,A)')  '  Pressure                  = ', P_kPa, ' kPa'
    write(*,*)
    write(*,'(A)') '--- STOICHIOMETRY -------------------------------------------'
    write(*,'(A,F10.4)')    '  O2 stoichiometric (mol)   = ', a_stoich
    write(*,'(A,F10.4)')    '  Air stoich (mol/mol fuel) = ', air_stoich
    write(*,'(A,F10.4)')    '  AFR stoichiometric (mass) = ', AFR_stoich
    write(*,'(A,F10.4)')    '  AFR actual (mass)         = ', AFR_actual
    write(*,'(A,F10.4)')    '  Equivalence Ratio phi     = ', 1.0d0/(1.0d0+excess_air/100.0d0)
    write(*,*)
    write(*,'(A)') '--- FLUE GAS COMPOSITION (per mole fuel) -------------------'
    write(*,'(A,F10.4,A)')  '  CO2                       = ', n_CO2, ' mol'
    write(*,'(A,F10.4,A)')  '  H2O                       = ', n_H2O, ' mol'
    write(*,'(A,F10.4,A)')  '  N2                        = ', n_N2, ' mol'
    write(*,'(A,F10.4,A)')  '  O2 (excess)               = ', n_O2_ex, ' mol'
    write(*,'(A,F10.4,A)')  '  Total products            = ', n_total, ' mol'
    write(*,*)
    write(*,'(A)') '--- FLUE GAS MOLE FRACTIONS ---------------------------------'
    write(*,'(A,F10.4,A)')  '  y_CO2                     = ', y_CO2*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_H2O                     = ', y_H2O*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_N2                      = ', y_N2*100.0d0, ' percent'
    write(*,'(A,F10.4,A)')  '  y_O2                      = ', y_O2*100.0d0, ' percent'
    write(*,*)
    write(*,'(A)') '--- ADIABATIC FLAME TEMPERATURE -----------------------------'
    write(*,'(A,F12.2,A)')  '  T_ad                      = ', T_ad, ' K'
    write(*,'(A,F12.2,A)')  '  T_ad                      = ', T_ad - 273.15d0, ' C'
    write(*,*)
    write(*,'(A)') '--- DEW POINT -----------------------------------------------'
    write(*,'(A,F10.4,A)')  '  P_H2O (partial)           = ', P_H2O, ' kPa'
    write(*,'(A,F12.2,A)')  '  T_dew                     = ', T_dew, ' K'
    write(*,'(A,F12.2,A)')  '  T_dew                     = ', T_dew - 273.15d0, ' C'
    write(*,*)

    ! ── Sensitivity sweep: excess air 0% to 200% ──────────────
    n_sweep = 40
    write(*,'(A)') '--- SENSITIVITY: T_AD VS EXCESS AIR -------------------------'
    write(*,'(A)') '  ExcessAir[%]  T_ad[K]       AFR_actual'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, n_sweep
        exc_sw = dble(i-1) * 200.0d0 / dble(n_sweep - 1)
        ! Recompute for this excess air
        actual_O2 = a_stoich * (1.0d0 + exc_sw / 100.0d0)
        n_CO2 = C_at
        n_H2O = H_at / 2.0d0
        n_N2  = 3.76d0 * actual_O2
        n_O2_ex = a_stoich * exc_sw / 100.0d0
        n_total = n_CO2 + n_H2O + n_N2 + n_O2_ex

        H_react = hf_fuel + (actual_O2 * cp_O2 + 3.76d0*actual_O2 * cp_N2) * &
                  (T_inlet - 298.15d0) / 1000.0d0

        T_low = 300.0d0; T_high = 5000.0d0
        do iter = 1, 80
            T_mid = 0.5d0*(T_low+T_high)
            H_prod = n_CO2*hf_CO2 + n_H2O*hf_H2O + &
                     (n_CO2*cp_CO2+n_H2O*cp_H2O+n_N2*cp_N2+n_O2_ex*cp_O2) * &
                     (T_mid-298.15d0)/1000.0d0
            if (H_prod < H_react) then; T_low=T_mid; else; T_high=T_mid; end if
            if (abs(T_high-T_low)<0.1d0) exit
        end do
        T_ad_sw = 0.5d0*(T_low+T_high)
        AFR_sw = AFR_stoich * (1.0d0 + exc_sw/100.0d0)
        write(*,'(F10.2,4X,F12.2,4X,F10.4)') exc_sw, T_ad_sw, AFR_sw
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Stoichiometric: CxHyOz + a(O2+3.76N2) -> CO2+H2O+N2'
    write(*,'(A)') '  a = C + H/4 - O_fuel/2'
    write(*,'(A)') '  AFR = a/0.21 * 28.97 / M_fuel (mass basis)'
    write(*,'(A)') '  T_ad by bisection on energy balance (constant cp)'
    write(*,'(A)') '  Dew point from partial pressure of H2O (simplified)'

end program combustion_aft


Solver Description

Calculates adiabatic flame temperature under constant pressure combustion conditions using thermodynamic enthalpy balances.

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 combustion_aft.f90 -o combustion_aft

Execution Command:

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

combustion_aft < input.txt

📥 Downloads & Local Files

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

! Fuel type (1=Methane, 2=Ethane, 3=Propane, 4=Butane, 5=Hydrogen, 6=Methanol, 7=Ethanol)
1
! Inlet temperature [K]
298.15
! Excess air percentage [%]
0.0
! Pressure P [atm]
1.0