line not appear in plot

Below are my code, why the plot not appear. someone help

% PEM Electrolyzer Activation Overvoltage Plot
% Constants
ioa = 2e-9;           % Anode exchange current density (A/cm2)
ioc = 2e-3;           % Cathode exchange current density (A/cm2)
T = 297;              % Temperature (K)
R = 8.314;             % General gas constatnt (J/Kmol)
ctca = 0.5 ;             % charge transfer coefficient anode
ctcc = 0.5 ;           % charge transfer coefficient cathode
z = 2 ;                % stoichiometric coefficient
F = 96485;             % Faradays constant
% Operating current density range
i_density = linspace(0, 1.5, 100);  % Current density range (A/m^2)
% Calculate activation overvoltage
a = R.*T/ctca.*z.*F;
b = R.*T/ctcc.*z.*F;
eta_activation = log((i_density/ioa).^(a))+(-log((i_density/ioc).^(b)))

eta_activation = 1×100

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

% Plot the results
plot(i_density, eta_activation,"b");

[Dyuman Joshi]

The plot is empty because all the values are NaN, see the edit above.

[amirul]

Thank you

Answer 1

[Star Strider]

The ‘i_density’ vector begins with 0 and the log of 0 is -Inf.

Starting it instead with a very small value, and re-writing ‘eta_activation’ to use simple log identities produces finite results.

There is some variation in ‘eta_activation’, as demonstrated by the derivative plot (added) —

% PEM Electrolyzer Activation Overvoltage Plot
% Constants
ioa = 2e-9;           % Anode exchange current density (A/cm2)
ioc = 2e-3;           % Cathode exchange current density (A/cm2)
T = 297;              % Temperature (K)
R = 8.314;             % General gas constatnt (J/Kmol)
ctca = 0.5 ;             % charge transfer coefficient anode
ctcc = 0.5 ;           % charge transfer coefficient cathode
z = 2 ;                % stoichiometric coefficient
F = 96485;             % Faradays constant
% Operating current density range
i_density = linspace(1E-12, 1.5, 100);  % Current density range (A/m^2)
% Calculate activation overvoltage
a = R.*T/ctca.*z.*F;
b = R.*T/ctcc.*z.*F;
eta_activation = log((i_density/ioa).^(a))+(-log((i_density/ioc).^(b)))     % Original

eta_activation = 1×100

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

eta_activation = a*log(i_density/ioa) - b*log(i_density/ioc)                % Rewritten

eta_activation = 1×100

1.0e+10 * 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166 1.3166

deta_activation_di_density = gradient(eta_activation,i_density)             % Derivative

deta_activation_di_density = 1×100

1.0e-03 * 0 0 0 -0.0629 -0.1259 0.0629 0.0629 -0.0629 0.0629 0 -0.0629 0.0629 0 -0.1259 0 0 0 0 0.0629 0 -0.0629 0.0629 0 0.0629 -0.0629 0 0 -0.0629 0.0629 0

% Plot the results
figure
plot(i_density, eta_activation,"b");
xlabel('i\_density')
ylabel('eta\_activation')

figure
plot(i_density, deta_activation_di_density,"g");
xlabel('i\_density')
ylabel('$\frac{d(eta\_activation)}{d(i\_density)}$', 'Interpreter','latex')

.

[VBBV]

@amirul For PEM electrolyzer, the forumula for activation energy is written in a different way.

The constants a and b are calculated incorrectly, see below. As mentioned by @Star Strider, use a small value at starting point,to avoid the Inf values while calculation.

% PEM Electrolyzer Activation Overvoltage Plot
% Constants
ioa = 2e-9;           % Anode exchange current density (A/cm2)
ioc = 2e-3;           % Cathode exchange current density (A/cm2)
T = 297;              % Temperature (K)
R = 8.314;             % General gas constatnt (J/Kmol)
ctca = 0.5 ;             % charge transfer coefficient anode
ctcc = 0.5 ;           % charge transfer coefficient cathode
z = 2 ;                % stoichiometric coefficient
F = 96485;             % Faradays constant
% Operating current density range
i_density = linspace(0, 1.5, 100);  % Current density range (A/m^2)
% Calculate activation overvoltage
a = R.*T/(ctca.*z.*F);  % use the parenthesis for denominator acc to formula 
b = R.*T/(ctcc.*z.*F);
eta_activation = log((i_density/ioa).^(a))-log((i_density/ioc).^(b))

eta_activation = 1×100

NaN 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536 0.3536

% Plot the results
plot(i_density, eta_activation,"b");

[Star Strider]

@VBBV — Thank you!

I would never have caught the errors in ‘a’ and ‘b’ (denominator factors need to be in parentheses).

[amirul]

Thank you very much my friends. Appriciate it

[Star Strider]

Our pleasure!