🔬
Solver Purpose & Physical Scope
Solves the electrochemical and thermodynamic equations of hydrogen-oxygen fuel cells (PEMFC, SOFC, AFC, DMFC). Computes Nernst reversible open-circuit potential, activation polarization (Tafel equation), ohmic polarization, concentration polarization losses, actual cell voltage, stack electrical power output, thermal efficiency, and rate of heat generation.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 241 times
📄 Source File:
fuel_cell_thermo.f90
📁 calcul/Thermodynamics /
fuel_cell_thermo.f90
program fuel_cell_thermo
implicit none
integer :: cell_type, n_cells, iostat_val, i, n_sweep
double precision :: T_cell, P_cell, i_dens, A_active
double precision :: P_H2, P_O2, alpha_a, i0_a, r_ohm, i_lim
double precision :: n_e, dG0, dH0, dS0, E0_rev, M_fuel
double precision :: F_const, R_gas
double precision :: E_rev_T, E_Nernst, V_act, V_ohm, V_conc, V_cell
double precision :: P_density, P_cell_W, P_stack
double precision :: eta_thermo, eta_voltage, eta_overall
double precision :: Q_heat, mdot_fuel
double precision :: i_sw, V_sw, P_sw
character(len=50) :: cell_name, fuel_name
F_const = 96485.0d0 ! C/mol
R_gas = 8.314462d0 ! J/(mol.K)
read(*,*,iostat=iostat_val) cell_type
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid cell type.'; stop; end if
read(*,*,iostat=iostat_val) T_cell
read(*,*,iostat=iostat_val) P_cell
read(*,*,iostat=iostat_val) i_dens
read(*,*,iostat=iostat_val) A_active
read(*,*,iostat=iostat_val) n_cells
read(*,*,iostat=iostat_val) P_H2
read(*,*,iostat=iostat_val) P_O2
read(*,*,iostat=iostat_val) alpha_a
read(*,*,iostat=iostat_val) i0_a
read(*,*,iostat=iostat_val) r_ohm
read(*,*,iostat=iostat_val) i_lim
if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if
if (T_cell <= 0.0d0) T_cell = 353.15d0
if (P_cell <= 0.0d0) P_cell = 1.0d0
if (i_dens <= 0.0d0) i_dens = 0.5d0
if (A_active <= 0.0d0) A_active = 100.0d0
if (n_cells < 1) n_cells = 1
if (P_H2 <= 0.0d0) P_H2 = 1.0d0
if (P_O2 <= 0.0d0) P_O2 = 0.21d0
if (alpha_a <= 0.0d0) alpha_a = 0.5d0
if (i0_a <= 0.0d0) i0_a = 1.0d-4
if (r_ohm <= 0.0d0) r_ohm = 0.15d0
if (i_lim <= 0.0d0) i_lim = 2.0d0
! ── Cell reaction data ─────────────────────────────────────
select case(cell_type)
case(1)
cell_name = 'PEM Fuel Cell (H2/O2)'
fuel_name = 'Hydrogen H2'
n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
M_fuel = 2.016d0
case(2)
cell_name = 'Solid Oxide Fuel Cell (SOFC, H2/O2)'
fuel_name = 'Hydrogen H2'
n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
M_fuel = 2.016d0
if (T_cell < 873.15d0) T_cell = 1073.15d0 ! SOFC operates at high T
case(3)
cell_name = 'Direct Methanol FC (DMFC)'
fuel_name = 'Methanol CH3OH'
n_e = 6.0d0; dG0 = -702.5d0; dH0 = -726.6d0; E0_rev = 1.213d0
M_fuel = 32.04d0
case default
cell_name = 'Alkaline Fuel Cell (AFC, H2/O2)'
fuel_name = 'Hydrogen H2'
n_e = 2.0d0; dG0 = -237.2d0; dH0 = -285.8d0; E0_rev = 1.229d0
M_fuel = 2.016d0
cell_type = 4
end select
! ── Thermodynamic calculations ─────────────────────────────
! Entropy change
dS0 = (dH0 - dG0) * 1000.0d0 / 298.15d0 ! J/(mol.K)
! Reversible voltage at temperature T
E_rev_T = E0_rev + dS0 / (n_e * F_const) * (T_cell - 298.15d0)
! Nernst equation (partial pressure correction)
! For H2 + 0.5O2 -> H2O:
! E_Nernst = E_rev(T) + (RT/(n_e*F))*ln(P_H2 * P_O2^0.5)
E_Nernst = E_rev_T + (R_gas * T_cell / (n_e * F_const)) * &
log(P_H2 * sqrt(P_O2))
! ── Polarization losses ────────────────────────────────────
! Activation overpotential (Butler-Volmer / Tafel)
if (i_dens > 0.0d0 .and. i0_a > 0.0d0) then
V_act = (R_gas * T_cell / (alpha_a * n_e * F_const)) * &
log(i_dens / i0_a + sqrt((i_dens/i0_a)**2 + 1.0d0))
else
V_act = 0.0d0
end if
! Ohmic losses
V_ohm = i_dens * r_ohm
! Concentration (mass transport) losses
if (i_dens < i_lim .and. i_dens > 0.0d0) then
V_conc = -(R_gas * T_cell / (n_e * F_const)) * log(1.0d0 - i_dens / i_lim)
else
V_conc = 0.5d0 ! saturated
end if
! Cell voltage
V_cell = E_Nernst - V_act - V_ohm - V_conc
if (V_cell < 0.0d0) V_cell = 0.0d0
! Power density
P_density = V_cell * i_dens ! W/cm2
! Cell and stack power
P_cell_W = V_cell * i_dens * A_active ! W
P_stack = P_cell_W * dble(n_cells) ! W
! Efficiencies
eta_thermo = abs(dG0 / dH0) ! max thermodynamic efficiency
if (E_Nernst > 1.0d-10) then
eta_voltage = V_cell / E_Nernst
else
eta_voltage = 0.0d0
end if
eta_overall = eta_voltage * eta_thermo
! Heat generation
Q_heat = (E_Nernst - V_cell) * i_dens * A_active * dble(n_cells) ! W
! Fuel consumption
mdot_fuel = i_dens * A_active * dble(n_cells) * M_fuel / (n_e * F_const) ! g/s
! ── Output ──────────────────────────────────────────────────
write(*,'(A)') '============================================================'
write(*,'(A)') ' FUEL CELL THERMODYNAMICS'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Cell Type = ', trim(cell_name)
write(*,'(A,A)') ' Fuel = ', trim(fuel_name)
write(*,'(A,F12.2,A)') ' Cell Temperature = ', T_cell, ' K'
write(*,'(A,F12.2,A)') ' Cell Temperature = ', T_cell-273.15d0, ' C'
write(*,'(A,F10.4,A)') ' Current Density = ', i_dens, ' A/cm2'
write(*,'(A,F12.2,A)') ' Active Area = ', A_active, ' cm2'
write(*,'(A,I8)') ' Number of Cells = ', n_cells
write(*,'(A,F10.4,A)') ' P_H2 = ', P_H2, ' atm'
write(*,'(A,F10.4,A)') ' P_O2 = ', P_O2, ' atm'
write(*,*)
write(*,'(A)') '--- ELECTROCHEMICAL PARAMETERS ------------------------------'
write(*,'(A,F10.4)') ' Transfer Coefficient a = ', alpha_a
write(*,'(A,ES12.4,A)') ' Exchange Current i0 = ', i0_a, ' A/cm2'
write(*,'(A,F10.4,A)') ' Ohmic Resistance r = ', r_ohm, ' ohm.cm2'
write(*,'(A,F10.4,A)') ' Limiting Current i_lim = ', i_lim, ' A/cm2'
write(*,*)
write(*,'(A)') '--- THERMODYNAMICS ------------------------------------------'
write(*,'(A,F12.2,A)') ' Delta_G0 = ', dG0, ' kJ/mol'
write(*,'(A,F12.2,A)') ' Delta_H0 = ', dH0, ' kJ/mol'
write(*,'(A,F12.4,A)') ' Delta_S0 = ', dS0, ' J/(mol.K)'
write(*,'(A,F10.6)') ' n_electrons = ', n_e
write(*,'(A,F10.6,A)') ' E0_rev (25C) = ', E0_rev, ' V'
write(*,'(A,F10.6,A)') ' E_rev(T) = ', E_rev_T, ' V'
write(*,'(A,F10.6,A)') ' E_Nernst = ', E_Nernst, ' V'
write(*,*)
write(*,'(A)') '--- POLARIZATION LOSSES -------------------------------------'
write(*,'(A,F10.6,A)') ' Activation V_act = ', V_act, ' V'
write(*,'(A,F10.6,A)') ' Ohmic V_ohm = ', V_ohm, ' V'
write(*,'(A,F10.6,A)') ' Concentration V_conc = ', V_conc, ' V'
write(*,'(A,F10.6,A)') ' Total Losses = ', V_act+V_ohm+V_conc, ' V'
write(*,*)
write(*,'(A)') '--- CELL PERFORMANCE ----------------------------------------'
write(*,'(A,F10.6,A)') ' Cell Voltage V_cell = ', V_cell, ' V'
write(*,'(A,F10.6,A)') ' Power Density = ', P_density, ' W/cm2'
write(*,'(A,F12.4,A)') ' Cell Power = ', P_cell_W, ' W'
write(*,'(A,F12.4,A)') ' Stack Power = ', P_stack, ' W'
write(*,'(A,F12.2,A)') ' Stack Power = ', P_stack/1000.0d0, ' kW'
write(*,*)
write(*,'(A)') '--- EFFICIENCY ----------------------------------------------'
write(*,'(A,F10.4)') ' Thermo Efficiency (max) = ', eta_thermo
write(*,'(A,F10.2,A)') ' Thermo Efficiency = ', eta_thermo*100.0d0, ' percent'
write(*,'(A,F10.4)') ' Voltage Efficiency = ', eta_voltage
write(*,'(A,F10.2,A)') ' Voltage Efficiency = ', eta_voltage*100.0d0, ' percent'
write(*,'(A,F10.4)') ' Overall Efficiency = ', eta_overall
write(*,'(A,F10.2,A)') ' Overall Efficiency = ', eta_overall*100.0d0, ' percent'
write(*,*)
write(*,'(A)') '--- HEAT & FUEL CONSUMPTION ---------------------------------'
write(*,'(A,F12.4,A)') ' Heat Generation = ', Q_heat, ' W'
write(*,'(A,F12.6,A)') ' Fuel Consumption = ', mdot_fuel, ' g/s'
write(*,'(A,F12.4,A)') ' Fuel Consumption = ', mdot_fuel*3600.0d0, ' g/h'
write(*,*)
! ── Polarization curve sweep ──────────────────────────────
n_sweep = 50
write(*,'(A)') '--- POLARIZATION CURVE --------------------------------------'
write(*,'(A)') ' i[A/cm2] V_cell[V] P_dens[W/cm2]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, n_sweep
i_sw = 0.01d0 + dble(i-1) * (i_lim*0.95d0 - 0.01d0) / dble(n_sweep-1)
if (i_sw > 0.0d0 .and. i0_a > 0.0d0) then
V_act = (R_gas*T_cell/(alpha_a*n_e*F_const)) * &
log(i_sw/i0_a + sqrt((i_sw/i0_a)**2+1.0d0))
else
V_act = 0.0d0
end if
V_ohm = i_sw * r_ohm
if (i_sw < i_lim) then
V_conc = -(R_gas*T_cell/(n_e*F_const))*log(1.0d0-i_sw/i_lim)
else
V_conc = 0.5d0
end if
V_sw = E_Nernst - V_act - V_ohm - V_conc
if (V_sw < 0.0d0) V_sw = 0.0d0
P_sw = V_sw * i_sw
write(*,'(F10.4,4X,F10.6,4X,F10.6)') i_sw, V_sw, P_sw
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Nernst: E = E0 + (dS/(nF))(T-298) + (RT/(nF))ln(PH2*PO2^0.5)'
write(*,'(A)') ' Activation: V_act = (RT/(a*n*F))*arcsinh(i/(2*i0))'
write(*,'(A)') ' Ohmic: V_ohm = i * R_area'
write(*,'(A)') ' Concentration: V_conc = -(RT/(nF))*ln(1 - i/i_lim)'
write(*,'(A)') ' Efficiency: eta = (dG/dH) * (V_cell/E_Nernst)'
end program fuel_cell_thermo
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 fuel_cell_thermo.f90 -o fuel_cell_thermo
2. Execution with input.txt redirection:
fuel_cell_thermo < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 353.15 1.0 0.6 200.0 100 1.5 0.21 0.5 1e-4 0.15 2.0
Parameter Description:
Cell type (1=PEM, 2=SOFC, 3=DMFC, 4=AFC) Cell temperature T [K] System pressure [atm] Current density [A/cm2] Active cell area [cm2] Number of cells in stack Hydrogen partial pressure pH2 [atm] Oxygen partial pressure pO2 [atm] Transfer coefficient alpha Exchange current density i0 [A/cm2] Area-specific resistance r [Ohm-cm2] Limiting current density ilim [A/cm2]