🔬
Solver Purpose & Physical Scope
Calculates adiabatic flame temperature under constant pressure combustion conditions using thermodynamic enthalpy balances.
📂 Discipline: Tools
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 375 times
📄 Source File:
combustion_aft.f90
📁 calcul/Tools /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 combustion_aft.f90 -o combustion_aft
2. Execution with input.txt redirection:
combustion_aft < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 298.15 0.0 1.0
Parameter Description:
Fuel type (1=Methane, 2=Ethane, 3=Propane, 4=Butane, 5=Hydrogen, 6=Methanol, 7=Ethanol) Inlet temperature [K] Excess air percentage [%] Pressure P [atm]