Fuel cell system efficiency increasing at increasing altitude

function fuel_cell_gui

% Initialize Figure 1 (Fuel Cell Performance)

fig = uifigure('Name', 'Fuel Cell Performance', 'Position', [100, 100, 1000, 650]);

% Axes for Graphs in Figure 1

ax1 = uiaxes(fig, 'Position', [50, 400, 400, 200]); % Polarization & Net Power

yyaxis(ax1, 'right'); ylabel(ax1, 'Power (W)'); % Power axis

yyaxis(ax1, 'left'); ylabel(ax1, 'Voltage (V)');

ax2 = uiaxes(fig, 'Position', [550, 400, 400, 200]); % H2 Flow Rates

ax3 = uiaxes(fig, 'Position', [50, 150, 400, 200]); % Air Flow Rates

ax4 = uiaxes(fig, 'Position', [550, 150, 400, 200]); % Efficiency

% Altitude Input for ISA Calculations

altitude_label = uilabel(fig, 'Position', [50, 90, 150, 20], 'Text', 'Altitude (m):');

altitude_input = uieditfield(fig, 'numeric', 'Position', [180, 90, 100, 25], 'Limits', [0, 20000], 'Value', 0);

altitude_input.ValueChangedFcn = @(s, e) updateFigures(); % Update figures when altitude changes

% Pressure Ratio Slider (Affects Both Temperature and Pressure)

pressure_slider = uislider(fig, 'Position', [100, 80, 300, 3], 'Limits', [1, 3], 'Value', 1.5);

pressure_slider.ValueChangedFcn = @(s, e) updateFigures(); % Update both figures

pressure_label = uilabel(fig, 'Position', [50, 90, 100, 20], 'Text', 'Pressure Ratio');

% Current Input

current_label = uilabel(fig, 'Position', [50, 20, 150, 20], 'Text', 'Input Current (A):');

current_input = uieditfield(fig, 'numeric', 'Position', [180, 20, 100, 25], 'Limits', [0.001, 3], 'Value', 1);

current_input.ValueChangedFcn = @(s, e) updateFigures(); % Update both figures

% Compressor Cycle Type Dropdown

compressor_cycle_label = uilabel(fig, 'Position', [50, 120, 150, 20], 'Text', 'Compressor Cycle:');

compressor_cycle_dropdown = uidropdown(fig, 'Position', [180, 120, 100, 25], ...

'Items', {'Isentropic', 'Polytropic', 'Multistage'}, 'Value', 'Isentropic');

compressor_cycle_dropdown.ValueChangedFcn = @(s, e) updateFigures(); % Update both figures

% Initialize Figure 2 (Fuel Cell Number of Cells)

fig_cells = uifigure('Name', 'Fuel Cell Number of Cells', 'Position', [1200, 100, 1000, 650]);

% Axes for Graphs in Figure 2

ax_power = uiaxes(fig_cells, 'Position', [50, 400, 400, 200]); % Power Curves

ax_air_flow = uiaxes(fig_cells, 'Position', [50, 150, 400, 200]); % Air Flow Rate

ax_efficiency = uiaxes(fig_cells, 'Position', [550, 150, 400, 200]); % System Efficiency

% Number of Cells Input

num_cells_label = uilabel(fig_cells, 'Position', [50, 20, 150, 20], 'Text', 'Number of Cells:');

num_cells_input = uieditfield(fig_cells, 'numeric', 'Position', [180, 20, 100, 25], 'Limits', [1, 10000], 'Value', 1000);

num_cells_input.ValueChangedFcn = @(s, e) updateFigure2(); % Update Figure 2 based on number of cells

% Shared variables for Figure 1 and Figure 2

V_cell = []; % Voltage per cell

i_values = []; % Current values

P_fuel_cell = []; % Fuel cell power (single cell)

P_comp = []; % Compressor power (single cell)

P_net = []; % Net power (single cell)

V_H2 = []; % Hydrogen flow rate

V_air = []; % Air flow rate

efficiency = []; % System efficiency

I_input = []; % Selected current input

idx_closest = []; % Closest index for selected current

% Initialize plots for both figures

updateFigures();

% Function to update both figures

function updateFigures()

updateFigure1(); % Update Figure 1

updateFigure2(); % Update Figure 2

end

% Function to update Figure 1 (independent of number of cells)

function updateFigure1()

% Constants

F = 96485;

R = 8.314;

E0 = 1.225;

P_base = 101325;

alpha = 0.5;

in = 0.01;

i0 = 2.5E-7;

omega = 0.091;

a = 3.5E-4;

b = 2.95;

stoich_ratio = 3;

gamma_air = 1.4;

R_air = 287;

% T_in = 298; % Ambient Temperature in Kelvin

% P_in = 101325; % Ambient Pressure in Pascals

% Get Input Values

P_comp_ratio = pressure_slider.Value;

I_input = current_input.Value;

compressor_cycle = compressor_cycle_dropdown.Value;

altitude = altitude_input.Value;

% Get Ambient Pressure and Temperature from ISA Model

[P_in, T_in] = getPressureAndTemperatureAtAltitude(altitude);

% Current Array

i_values = linspace(0.001, 2.2, 100); % Update current range to 3A

% Initialize Results

V_cell = zeros(size(i_values));

V_H2 = zeros(size(i_values));

V_air = zeros(size(i_values));

P_net = zeros(size(i_values));

P_fuel_cell = zeros(size(i_values));

P_comp = zeros(size(i_values));

efficiency = zeros(size(i_values));

% Conversion Factor for SLPM (Standard Liters per Minute)

conversion_factor = 22.414 * 60;

% Iterative Loop to Update Temperature and Polarization Curve

T_out = T_in; % Initial guess for outlet temperature

for iter = 1:10 % Perform 10 iterations (can be adjusted for convergence)

% Calculate Values

for idx = 1:length(i_values)

i = i_values(idx);

E_cell = E0 + (R * T_out) / (2 * F) * log((P_in * P_comp_ratio / P_base)^(1/2));

% Updated voltage equation with Frano-Barbir concentration loss

% Frano-Barbir loss equation with water vapor consideration

% V = a * (x_O2_amb * P_amb / x_O2_exit * P_exit) * exp(b * i)

x_O2_amb = 0.21; % Oxygen mole fraction (assumed ambient air)

x_O2_exit = 0.21; % Oxygen mole fraction (assumed exit)

V_loss = a * (x_O2_amb * P_in / (x_O2_exit * P_in)) * exp(b * i);

S_Volt = E_cell - (R * T_out) / (2 * alpha * F) * log(max((i + in), 1e-12) / i0) - i * omega - V_loss;

S_Power = i * S_Volt;

P_fuel_cell(idx) = S_Power; % Fuel cell power for a single cell

% Mass flow rates

M_H2 = 2.016e-3; % kg/mol

M_air = 28.97e-3; % kg/mol

mdot_H2 = (i * 1 / (2 * F)) * M_H2;

mdot_air = (i * stoich_ratio / (0.21 * 4 * F)) * M_air;

% Convert mol/s to SLPM

V_H2(idx) = mdot_H2 * conversion_factor; % Flow rate in SLPM

V_air(idx) = mdot_air * conversion_factor; % Flow rate in SLPM

% Compressor power calculation based on selected cycle type

switch compressor_cycle

case 'Isentropic'

compressor_efficiency = 0.5;

Cp_air = 1005; % J/kg-K (instead of R_air)

%P_comp(idx) = (mdot_air * Cp_air * T_in) / compressor_efficiency * ((P_comp_ratio^((gamma_air - 1) / gamma_air)) - 1);

P_comp(idx) = (mdot_air * R_air * T_in) / (compressor_efficiency * (gamma_air - 1)) * ((P_comp_ratio)^((gamma_air - 1) / gamma_air) - 1);

case 'Polytropic'

n = 1.3; % Polytropic exponent

P_comp(idx) = (mdot_air * R_air * T_in) / (n - 1) * ((P_comp_ratio)^((n - 1) / n) - 1);

case 'Multistage'

% Assuming 2-stage compression with intercooling

P_comp_ratio_stage = sqrt(P_comp_ratio);

compressor_efficiency = 0.5;

P_comp(idx) = 2 * (mdot_air * R_air * T_in) / (compressor_efficiency * (gamma_air - 1)) * ((P_comp_ratio_stage)^((gamma_air - 1) / gamma_air) - 1);

end

P_net(idx) = S_Power - P_comp(idx); % Net power = fuel cell power - compressor power

efficiency(idx) = (S_Power / (S_Power + P_comp(idx)))*100; % System efficiency

V_cell(idx) = S_Volt;

end

% Update Outlet Temperature

T_out = T_in * (P_comp_ratio^((gamma_air - 1) / gamma_air));

end

% Find closest index for selected current input

[~, idx_closest] = min(abs(i_values - I_input));

% Plot and add red circles for selected values in Figure 1

% Polarization & Power Curves

plot(ax1, i_values, V_cell, 'b', i_values, P_net, 'k', i_values, P_fuel_cell, 'r', 'LineWidth', 2);

hold(ax1, 'on');

plot(ax1, I_input, V_cell(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Voltage

plot(ax1, I_input, P_net(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Net Power

hold(ax1, 'off');

grid(ax1, 'on');

title(ax1, 'Polarization & Power Curves');

xlabel(ax1, 'Current (A)');

ylabel(ax1, 'Voltage (V) & Power (W)');

legend(ax1, 'Voltage', 'Net Power', 'Fuel Cell Power');

% H2 Flow Rates

plot(ax2, i_values, V_H2, 'b', 'LineWidth', 2);

hold(ax2, 'on');

plot(ax2, I_input, V_H2(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for H2 Flow

hold(ax2, 'off');

grid(ax2, 'on');

title(ax2, 'Hydrogen Flow Rate');

xlabel(ax2, 'Current (A)');

ylabel(ax2, 'Flow Rate (SLPM)');

legend(ax2, 'H2 Flow');

% Air Flow Rates

plot(ax3, i_values, V_air, 'b', 'LineWidth', 2);

hold(ax3, 'on');

plot(ax3, I_input, V_air(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Air Flow

hold(ax3, 'off');

grid(ax3, 'on');

title(ax3, 'Air Flow Rate');

xlabel(ax3, 'Current (A)');

ylabel(ax3, 'Flow Rate (SLPM)');

legend(ax3, 'Air Flow');

% Efficiency

plot(ax4, i_values, efficiency, 'b', 'LineWidth', 2);

hold(ax4, 'on');

plot(ax4, I_input, efficiency(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Efficiency

hold(ax4, 'off');

grid(ax4, 'on');

title(ax4, 'Efficiency');

xlabel(ax4, 'Current (A)');

ylabel(ax4, 'Efficiency (%)');

legend(ax4, 'Efficiency');

% Display values in the command window for red markers

fprintf('Selected Current: %.4f A\n', I_input);

fprintf('Voltage at %.4f A: %.4f V\n', I_input, V_cell(idx_closest));

fprintf('Net Power at %.4f A: %.4f W\n', I_input, P_net(idx_closest));

fprintf('Hydrogen Flow Rate at %.4f A: %.4f SLPM\n', I_input, V_H2(idx_closest));

fprintf('Air Flow Rate at %.4f A: %.4f SLPM\n', I_input, V_air(idx_closest));

fprintf('Efficiency at %.4f A: %.4f %%\n\n', I_input, efficiency(idx_closest));

end

% Function to update Figure 2 (based on number of cells)

function updateFigure2()

% Constants

num_cells = num_cells_input.Value;

% Power Calculations

P_fuel_cell_cells = num_cells * P_fuel_cell;

P_net_cells = num_cells * P_net;

V_cell_cells = V_cell;

V_H2_cells = num_cells * V_H2;

V_air_cells = num_cells * V_air;

efficiency_cells = efficiency; % Efficiency based on number of cells

% Plotting the Results for Power Curves (Figure 2)

plot(ax_power, i_values, P_fuel_cell_cells, 'r', i_values, P_net_cells, 'k', 'LineWidth', 2);

hold(ax_power, 'on');

plot(ax_power, I_input, P_fuel_cell_cells(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Power

hold(ax_power, 'off');

grid(ax_power, 'on');

title(ax_power, 'Power Curves (Number of Cells)');

xlabel(ax_power, 'Current (A)');

ylabel(ax_power, 'Power (W)');

legend(ax_power, 'Fuel Cell Power', 'Net Power');

% Air Flow Rates (Figure 2)

plot(ax_air_flow, i_values, V_air_cells, 'b', 'LineWidth', 2);

hold(ax_air_flow, 'on');

plot(ax_air_flow, I_input, V_air_cells(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Air Flow

hold(ax_air_flow, 'off');

grid(ax_air_flow, 'on');

title(ax_air_flow, 'Air Flow Rate (Number of Cells)');

xlabel(ax_air_flow, 'Current (A)');

ylabel(ax_air_flow, 'Flow Rate (SLPM)');

legend(ax_air_flow, 'Air Flow');

% Efficiency (Figure 2)

plot(ax_efficiency, i_values, efficiency_cells, 'b', 'LineWidth', 2);

hold(ax_efficiency, 'on');

plot(ax_efficiency, I_input, efficiency_cells(idx_closest), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Red marker for Efficiency

hold(ax_efficiency, 'off');

grid(ax_efficiency, 'on');

title(ax_efficiency, 'Efficiency (Number of Cells)');

xlabel(ax_efficiency, 'Current (A)');

ylabel(ax_efficiency, 'Efficiency (%)');

legend(ax_efficiency, 'Efficiency');

end

% Function to calculate ISA-based pressure and temperature at altitude

function [P, T] = getPressureAndTemperatureAtAltitude(h)

% Constants for ISA (International Standard Atmosphere)

T0 = 288.15; % Standard temperature in Kelvin

P0 = 101325; % Standard pressure in Pa

L = 0.00649; % Temperature lapse rate

R = 287; % Specific gas constant for air (J/(kg*K))

g0 = 9.81; % Gravitational acceleration (m/s^2)

T = T0 - L * h; % Temperature at altitude

P = P0 * (1 - L * h / T0)^(g0 / (R * L)); % Pressure at altitude

end

end

Answer 1

[Sam Chak]

Hi @Vyom Desai

We are uncertain where to begin in reviewing the code, as all the essential calculative steps for determining system efficiency (that depends on S_Power and P_comp) have not been provided, which prevents us from making a comparison.

Setting that aside, if it seems illogical for efficiency to increase proportionally with altitude, please verify whether a 'minus' sign is missing in some of the terms associated with altitude, such as air pressure and temperature in the code. Both air pressure and temperature should decrease as altitude increases.

When these two variables decrease, do E_cell, V_loss, S_Volt, and P_comp behave as expected?

efficiency(idx) = (S_Power / (S_Power + P_comp(idx)))*100;  % System efficiency

Since you imply that the efficiency should decrease as altitude increases, if your given efficiency formula is correct, we can Sherlock Holmesly deduce that the S_Power should decrease and P_comp should increase as altitude increases.

How can this be verified? I recommend plotting the graph and evaluating the results from the perspective of a fuel cell expert.

[Vyom Desai]

From the polarisation curve, a voltage is provided. That is multiplied by the current to calculate S_Power. This is the power of the system. To calculate the P_Comp is based on which thermodynamic cycle is being used for the compressor. For both of these cases, they are just equations.

When the variables change, E_Cell, V_loss and S_Volt all change as expected as they are affected by the temperature. P_Comp is different, that is directly affected by the required mass flow rate of the fuel cell. in theory as altitude increases, the mass flow rate increases hence the P_comp. but this is not the case.

[Sam Chak]

Hi @Vyom Desai

Can you edit the code (make it GUI-free) for directly plotting out air pressure, temperature, E_Cell, V_loss and S_Volt on MATLAB forum? It helps to see the changes as altitude increases.