% 0D_SOEC_model.

% Simple 0D Solid Oxide Electrolysis Cell (SOEC) model (H2O electrolysis)

% - Steady-state evaluation for given current density

% - Simple dynamic temperature ODE

% - Includes: Nernst, activation (Tafel), ohmic, concentration

overpotentials

% - Outputs cell voltage, outlet compositions, and temperature evolution

% Usage:

% - Edit the 'inputs' section below to set operating conditions.

% - Run the script. It will compute steady-state cell voltage for the

% prescribed current density and then simulate a short transient for T.

% Notes:

% - This is a 0D single-cell model (no spatial discretization).

% - It treats electrochemical production/consumption via Faraday's law.

% - Mass transport losses are modelled with a simple limiting-current

approach.

% - The model currently handles only H2O <-> H2 + 1/2 O2 electrolysis.

To

% include CO2 co-electrolysis, extend the stoichiometry and Faradaic

split.

clear; close all; clc

F = 96485.3329; R = 8.314462618; % constants

%% ---------------------- INPUTS ----------------------

% Cell / geometry

A_cell = 0.01; % active area [m^2]

n_cells = 1; % number of repeating cells (stack) - affects area

LUT Group Confidential - Other information (3Y)

 % Operating conditions

T0 = 900 + 273.15; % operating temperature [K]

P = 1e5; % pressure [Pa]

% Gas inlet molar flows (mol/s)

M_dot_H2O_in = 0.01; % steam inlet mol/s

M_dot_H2_in = 0.0; % hydrogen inlet (purge / recirculation)

M_dot_inert = 0.0; % inert (N2) inlet

% Electrochemical operating point

i_app = 0.5; % applied current density [A/cm^2] -> convert later

i_app = i_app * 1e4; % [A/m^2]

% Electrolyte / materials

L_elec = 10e-6; % electrolyte thickness [m]

sigma_ion_ref = 1.0; % ionic conductivity at T_ref [S/m]

T_ref = 1000 + 273.15; % reference temperature for sigma

Ea_sigma = 40e3; % activation energy for ionic conductivity [J/mol]

% Kinetics (exchange current density)

i0_ref = 1e4; % reference exchange current density [A/m^2] at T_ref

Ea_i0 = 80e3; % activation energy for i0 [J/mol]

alpha = 0.5; % transfer coefficient

% Concentration / mass transport

i_lim = 2e5; % limiting current density [A/m^2] (simple input)

% Thermal / energy

LUT Group Confidential - Other information (3Y)

 m_cell = 0.1; % mass of the cell (kg) - for transient T

Cp_cell = 500; % heat capacity [J/kg-K]

U_loss = 10; % overall heat loss coefficient [W/K] to ambient

T_amb = 298.15; % ambient temperature [K]

%% ---------------------- PREPARE & DERIVED ----------------------

A_tot = A_cell * n_cells; % total active area [m^2]

% Compute temperature-dependent properties

sigma_ion = sigma_ion_ref * exp(-Ea_sigma/R*(1./T0 - 1./T_ref));

i0 = i0_ref * exp(-Ea_i0/R*(1./T0 - 1./T_ref));

% Stoichiometry for H2O electrolysis: H2O + 2e- -> H2 + O2- (per 2 e-)

n_e = 2; % electrons per H2 produced

%% ---------------------- MASS BALANCE (0D, steady-state) ----------------------

% Reaction molar rate (mol/s) from applied current:

R_reaction = (i_app * A_tot) / (n_e * F); % mol H2 produced per second

% Outlet flows (steady-state, neglecting accumulation)

M_dot_H2_out = M_dot_H2_in + R_reaction;

M_dot_H2O_out = M_dot_H2O_in - R_reaction;

if M_dot_H2O_out < 0

warning('Consumed more H2O than inlet — increase inlet steam or

reduce current');

M_dot_H2O_out = max(M_dot_H2O_out, 1e-12);

end

M_dot_O2_out = 0.5 * R_reaction;

M_dot_total_out = M_dot_H2_out + M_dot_H2O_out + M_dot_O2_out +

M_dot_inert;

LUT Group Confidential - Other information (3Y)

 X_H2 = M_dot_H2_out / M_dot_total_out;

X_H2O = M_dot_H2O_out / M_dot_total_out;

X_O2 = M_dot_O2_out / M_dot_total_out;

p_H2 = X_H2 * P; p_H2O = X_H2O * P; p_O2 = X_O2 * P;

%% ---------------------- ELECTROCHEMICAL MODEL ----------------------

% 1) Nernst (reversible) potential for H2O/H2

E0 = 0; % for H2O/H2 at standard (we'll use Nernst formulation below)

% Nernst potential for reaction: H2O <-> H2 + 0.5 O2

E_rev = (R*T0)/(n_e*F) * log( (p_H2 * sqrt(p_O2)) / p_H2O );

% 2) Activation overpotential (Tafel-like approximate)

% Use inverse Butler-Volmer approx: eta_act = (RT/alphaF) * asinh(i/(2*i0))

eta_act = (R*T0/(alpha*F)) * asinh(i_app./(2*i0));

% 3) Ohmic loss

eta_ohm = i_app .* (L_elec ./ (sigma_ion));

% 4) Concentration loss (simple limiting-current expression)

if i_app >= 0.999*i_lim

eta_conc = 10; % very large (practically unreachable)

else

eta_conc = - (R*T0/(n_e*F)) * log(1 - i_app./i_lim);

end

% Total cell voltage (electrolyzer mode: V = E_rev + sum(overpotentials))

V_cell = E_rev + eta_act + eta_ohm + eta_conc;

LUT Group Confidential - Other information (3Y)

 %% ---------------------- HEAT GENERATION (steady-state) ----------------------

% Reversible heat (entropy term) approximation: q_rev_per_mol = -T *

dS_rxn

% Use Gibbs relation: dG = dH - T dS => T dS = dH - dG

% For simplicity, we approximate TdS term via standard molar values is

omitted;

% instead include irreversible heat from overpotentials: Q_irr = i*A *

sum(etas)

Q_irrev = i_app * A_tot * (eta_act + eta_ohm + eta_conc); % [W]

% Also include electrochemical reaction enthalpy term if desired (omitted

here)

%% ---------------------- OUTPUTS ----------------------

fprintf('--- 0D SOEC steady-state result ---\n');

fprintf('T = %.1f K (%.1f C)\n', T0, T0-273.15);

fprintf('Applied current density = %.2f A/cm^2\n', i_app/1e4);

fprintf('Cell voltage = %.3f V (E_rev = %.3f V)\n', V_cell, E_rev);

fprintf('Overpotentials: eta_act=%.3f V, eta_ohm=%.3f V, eta_conc=%.3f

V\n', eta_act, eta_ohm, eta_conc);

fprintf('Outlet molar flows (mol/s): H2=%.5f, H2O=%.5f, O2=%.5f\n',

M_dot_H2_out, M_dot_H2O_out, M_dot_O2_out);

fprintf('Mole fractions: X_H2=%.3f, X_H2O=%.3f, X_O2=%.3f\n', X_H2,

X_H2O, X_O2);

fprintf('Irreversible heat generation ~ %.2f W\n', Q_irrev);

%% ---------------------- TRANSIENT TEMPERATURE (simple ODE)

----------------------

% dT/dt = (Q_irrev - U_loss*(T-T_amb)) / (m*Cp)

Tspan = [0 2000];

Tinit = T0;

LUT Group Confidential - Other information (3Y)

 opts = odeset('RelTol',1e-6,'AbsTol',1e-8);

[Tsol, Ysol] = ode15s(@(t,y)

tempODE(t,y,Q_irrev,U_loss,T_amb,m_cell,Cp_cell), Tspan, Tinit, opts);

figure; plot(Tsol-273.15, 'LineWidth',1.5); xlabel('Time step'); ylabel('T

(°C)'); title('Simple transient of cell temperature'); grid on;

%% ---------------------- FUNCTIONS ----------------------

function dTdt = tempODE(~, T, Q_irrev, U_loss, T_amb, m_cell, Cp_cell)

% Simple lumped thermal ODE

Q_loss = U_loss * (T - T_amb);

dTdt = (Q_irrev - Q_loss) / (m_cell * Cp_cell);

end

LUT Group Confidential - Other information (3Y)