Alkaline Membrane Electrolyser

Generated from runElectrolyser.m

The implementation of a AEM simulation code was motivated by the need for better simulation capabilities for H2 electrolysis and the possibility to validate the code using an implementation we had access to, which is also the one used in [GBK+23]. The mathematical model for AEM is complex andincludes many coupledprocesses . We refer to the papers mentioned above for all the details. Here, we include a figure from [GBK+23] which illustrates the following processes

Two parallel ionic conduction path ways: ionomer and liquid electrolyte,
Water concentration gradient build-up in the membrane inducing water exchange between electrodes,
Water exchange between the liquid and gas phases in the electrolyte,
Hydroxide exchange between the electrolyte and ionomer.

../_images/electolysermodeling.png — Illustration taken from [GBK+23]

We use a model hierarchy which enables us to structure the code in a hopefully clearly readable way. The main model is split into three models

Oxygen Evolution Electrode,
Ionomer Membrane,
Hydrogen Evolution Electrode.

The electrode models share the same structure given by

Porous Transport Layer,
Catalyst Layer,
Exchange Reactions Model,

The governing equations consists of charge and mass conservations equations with exchange source terms, see Figure. They are assembled at each model level for the core parts while the coupling terms are assembled at levels above in the model hierarchy. We use generic discrete spatial differentiation operators that can be applied to any dimension, even if the geometry we have setup for the result in this model is only 1D. Once the equations are discretized in space, we discretize in time and solve the resulting non-linear system of equations using an implicit scheme. This approach is robust independently of the choice oftime step size. We use automatic differentiation to assemble the equations and compute the derivative of the system that are needed in the Newton algorithmused to solve the equations. These computation steps are in fact generic and we rely on the BattMo infrastructure to run them.

We consider a one-dimensional electrolyser cell where the volume fractions of each component are given in Figure below. We increase the current linearly from 0 to 3 A/cm^2. The pressures at the boundary are set to 1 atmosphere and the \(OH^{-}\) concentration to 1 mol/litre.

../_images/electrolyter_volume_fractions.png — Volume fractions used initially in the cell. Illustration picture taken from [GBK+23].

The input datas are given in a json format. A simple way to get acquainted to the format is to look at the examples used in the example below. We provide also json schemas that describes these inputs:

Catalyst Layer CatalystLayer.schema.json
Electrolyser Electrolyser.schema.json
Evolution Electrode EvolutionElectrode.schema.json
Exchange Reaction ExchangeReaction.schema.json
Ionomer Membrane IonomerMembrane.schema.json
Porous Transport Layer PorousTransportLayer.schema.json
Electrolyser Geometry ElectrolyserGeometry.schema.json

Here comes the listing of the code:

Load MRST modules

mrstModule add ad-core matlab_bgl

Setup input

Setup the physical properties for the electrolyser using json input file alkalineElectrolyser.json

jsonstruct= parseBattmoJson('Electrolyser/Parameters/alkalineElectrolyser.json');
inputparams = ElectrolyserInputParams(jsonstruct);

Setup the grids. We consider a 1D model and the specifications can be read from a json input file, here electrolysergeometry1d.json, using setupElectrolyserGridFromJson.

jsonstruct= parseBattmoJson('Electrolyser/Parameters/electrolysergeometry1d.json');
inputparams = setupElectrolyserGridFromJson(inputparams, jsonstruct);

We define shortcuts for the different submodels.

inm = 'IonomerMembrane';
her = 'HydrogenEvolutionElectrode';
oer = 'OxygenEvolutionElectrode';
ptl = 'PorousTransportLayer';
exr = 'ExchangeReaction';
ctl = 'CatalystLayer';

Setup model

model = Electrolyser(inputparams);

Setup the initial condition

We use the default initial setup implemented in the model

[model, initstate] = model.setupBcAndInitialState();

Setup the schedule with the time discretization

We run the simulation over 10 hours, increasing the input current linearly in time.

total = 10*hour;

n   = 100;
dt  = total/n;
dts = rampupTimesteps(total, dt, 5);

We use the function rampupControl to increase the current linearly in time

controlI = -3*ampere/(centi*meter)^2; % if negative, O2 and H2 are produced
tup      = total;
srcfunc  = @(time) rampupControl(time, tup, controlI, 'rampupcase', 'linear');
control  = struct('src', srcfunc);

step = struct('val', dts, 'control', ones(numel(dts), 1));
schedule = struct('control', control, 'step', step);

Setup the non-linear solver

We do only minor modifications here from the standard solver

nls = NonLinearSolver();
nls.verbose = false;
nls.errorOnFailure = false;

model.verbose = false;

Run the simulation

[~, states, report] = simulateScheduleAD(initstate, model, schedule, 'NonLinearSolver', nls, 'OutputMiniSteps', true);

Visualize the results

The results contain only the primary variables of the system (the unknwons that descrive the state of the system). We use the method addVariables to add all the intermediate quantities that are computed to solve the equations but not stored automatically in the result.

for istate = 1 : numel(states)
    states{istate} = model.addVariables(states{istate});
end

We extract the time, voltage and current values for each time step

time = cellfun(@(state) state.time, states);
E    = cellfun(@(state) state.(oer).(ptl).E, states);
I    = cellfun(@(state) state.(oer).(ctl).I, states);

We plot the results for the voltage and current

set(0, 'defaultlinelinewidth', 3)
set(0, 'defaultaxesfontsize', 15)

figure
subplot(2, 1, 1)
plot(time/hour, E)
xlabel('time [hour]');
ylabel('voltage');
title('Polarisation curve');

subplot(2, 1, 2)
plot(time/hour, -I/(1/(centi*meter)^2));
xlabel('time [hour]');
ylabel('Current [A/cm^2]');
title('Input current')

pH distribution plot

We consider the three domains and plot the pH in each of those. We setup the helper structures to iterate over each domain for the plot.

models = {model.(oer).(ptl), ...
          model.(her).(ptl), ...
          model.(inm)};

fields = {{'OxygenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}  , ...
          {'HydrogenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}, ...
          {'IonomerMembrane', 'cOH'}};

h = figure();
set(h, 'position', [10, 10, 800, 450]);
hold on

ntime = numel(time);
times = linspace(1, ntime, 10);
cmap  = cmocean('deep', 10);

for ifield = 1 : numel(fields)

    fd       = fields{ifield};
    submodel = models{ifield};

    x    = submodel.grid.cells.centroids;

    for itimes = 1 : numel(times);

        itime = floor(times(itimes));
        % The method :code:`getProp` is used to recover the value from the state structure
        val   = model.getProp(states{itime}, fd);
        pH    = 14 + log10(val/(mol/litre));

        % plot of pH for the current submodel.
        plot(x/(milli*meter), pH, 'color', cmap(itimes, :));

    end

end

xlabel('x  /  mm');
ylabel('pH');
title('pH distribition in electrolyser')

colormap(cmap)
hColorbar = colorbar;
caxis([0 3]);
hTitle = get(hColorbar, 'Title');
set(hTitle, 'string', 'J (A/cm^2)');

complete source code can be found here