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