%% PEM electrolyser with Gas Supply
%
%% json input data
%
% We load the input data that is given by json structures. The physical properties of the supply gas layer is given in
% the json file :battmofile:`gas-supply-whole-cell.json <ProtonicMembrane/jsonfiles/gas-supply-whole-cell.json>`
clear jsonstruct_material
filename = 'ProtonicMembrane/jsonfiles/gas-supply-whole-cell.json';
jsonstruct_material.GasSupply = parseBattmoJson(filename);
filename = 'ProtonicMembrane/jsonfiles/protonicMembrane.json';
jsonstruct_material.Electrolyser = parseBattmoJson(filename);
jsonstruct_material = removeJsonStructFields(jsonstruct_material, ...
{'Electrolyser', 'Electrolyte', 'Nx'}, ...
{'Electrolyser', 'Electrolyte', 'xlength'});
%%
% The geometrical property of the cell is given in :battmofile:`2d-cell-geometry.json <ProtonicMembrane/jsonfiles/2d-cell-geometry.json>`
%
filename = 'ProtonicMembrane/jsonfiles/2d-cell-geometry.json';
jsonstruct_geometry = parseBattmoJson(filename);
%%
% The initial state is given in :battmofile:`gas-supply-initialization.json <ProtonicMembrane/jsonfiles/gas-supply-initialization.json>`
%
filename = 'ProtonicMembrane/jsonfiles/gas-supply-initialization.json';
jsonstruct_initialization.GasSupply = parseBattmoJson(filename);
%%
% The time stepping parameters are given in :battmofile:`cell-timestepping.json <ProtonicMembrane/jsonfiles/cell-timestepping.json>`
%
filename = 'ProtonicMembrane/jsonfiles/cell-timestepping.json';
jsonstruct_timestepping = parseBattmoJson(filename);
%%
% We merge all the json struct we have created to obtain a complete json structure that can be use to setup the model
%
jsonstruct = mergeJsonStructs({jsonstruct_material , ...
jsonstruct_geometry , ...
jsonstruct_initialization, ...
jsonstruct_timestepping});
%%
% We change the value of the current density to the value we want to use in this simulation
jsonstruct.Electrolyser.Control.currentDensity = 0.1*ampere/((centi*meter)^2);
%% Input parameter setup
%
% We setup the input parameter structure which will we be used to instantiate the model
%
inputparams = ProtonicMembraneCellInputParams(jsonstruct);
%%
%
% We setup the grid, by calling the function :battmo:`setupProtonicMembraneGasLayerGrid`
%
[inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);
%% Model setup
%
% We setup the model for the full cell (electrolyser and a gas supply layer).
model = ProtonicMembraneCell(inputparams);
%%
% The model is equipped for simulation using the following command (this step may become unnecessary in future versions)
%
model = model.setupForSimulation();
%% Grid plots
%
% We plot the different regions with their grids.
%%
% To resolve the non-linearity in the electrolyte, we need a very fine mesh. Plotting the full mesh including the edges
% for this domain will hide its structure. Therefore we set the option :code:`edgecolor` to :code:`none`.
%
figure('position', [337, 757, 3068, 557])
hold on
plotGrid(model.grid, 'edgecolor', 'none')
plotGrid(model.Electrolyser.grid, 'facecolor', 'red', 'edgecolor', 'none')
plotGrid(model.GasSupply.grid, 'facecolor', 0.9*[0 1 1], 'edgealpha', 0.2)
%%
% We plot the boundary and interface faces
s = spring(10);
opts = {'linewidth', 10, ...
'edgecolor', s(8, :)};
%%
% we plot the inlet faces
figure('position', [337, 757, 3068, 557])
plotGrid(model.Electrolyser.grid, 'facecolor', 'red', 'edgecolor', 'none')
plotGrid(model.GasSupply.grid, 'facecolor', 'none', 'edgealpha', 0.2)
plotFaces(model.GasSupply.grid, model.GasSupply.couplingTerms{1}.couplingfaces, opts{:});
%%
% we add the outlet faces
plotFaces(model.GasSupply.grid, model.GasSupply.couplingTerms{2}.couplingfaces, opts{:});
%%
% we add the interface faces
plotFaces(model.GasSupply.grid, model.couplingTerm.couplingfaces(:, 1) , opts{:});
%% Setup initial state
%
% The initial state for the cell is used using the initial state of the gas layer that have been loaded (see
% :battmofile:`here <ProtonicMembrane/jsonfiles/gas-supply-initialization.json>`)
initstate = model.setupInitialState(jsonstruct);
%% Setup schedule
%
% We setup the time stepping. We are computing the steady state solution, using the the time evolution equation. Note
% that we also use a numerical method which gradually introduces the high nonlinearities inherent to the problem (in
% particular the exponential dependence of the electronic conductivity in the membrane). Therefore, the intermediate
% solutions (i.e. those computed before the final step) should be considered with care.
schedule = model.setupSchedule(jsonstruct);
%% Setup nonlinear solver
%
% We use an adaptive time stepping. In this case, the solver will try to keep the number of timesteps equal to 5, based
% on the timestepping history.
ts = IterationCountTimeStepSelector('targetIterationCount', 5);
nls = NonLinearSolver();
nls.timeStepSelector = ts;
nls.maxIterations = 15;
%% Start simulation
%
% We start the simulation
[~, states, report] = simulateScheduleAD(initstate, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls);
%% plotting setup
%
%
close all
set(0, 'defaultlinelinewidth', 3);
set(0, 'defaultaxesfontsize', 15);
set(0, 'defaultfigureposition', [1290, 755, 1275, 559])
elyser = 'Electrolyser';
elyte = 'Electrolyte';
gs = 'GasSupply';
state = states{end};
state = model.addVariables(state);
%% Electrolyte results
% We plot the electrostatic potential :math:`\phi` (we need to remove the plotting of the grid edges, otherwise due to
% the grid refinement, they will hide completely the color of the field)
%
figure
plotCellData(model.(elyser).(elyte).grid, state.(elyser).(elyte).phi, 'edgecolor', 'none');
title('Electrostatic Potential \phi / V')
colorbar
%%
%
% We want also to plot the result as a surface plot. Then, we need to manipulate the output and make it suitable to the
% :code:`surf` function of matlab (we do not explain those details here)
N = gen.nxElectrolyser;
xc = model.(elyser).(elyte).grid.cells.centroids(1 : N, 1);
X = reshape(model.(elyser).(elyte).grid.cells.centroids(:, 1), N, [])/(milli*meter);
Y = reshape(model.(elyser).(elyte).grid.cells.centroids(:, 2), N, [])/(milli*meter);
figure
val = state.(elyser).(elyte).phi;
Z = reshape(val, N, []);
surf(X, Y, Z, 'edgecolor', 'none');
title('\phi / V')
xlabel('x / mm')
ylabel('y / mm')
view(45, 31)
colorbar
%%
%
% We produce a surface plot of the electromotive potential :math:`\pi`.
figure
val = state.(elyser).(elyte).pi;
Z = reshape(val, N, []);
surf(X, Y, Z, 'edgecolor', 'none');
title('pi / V')
xlabel('x / mm')
ylabel('y / mm')
view(45, 31)
colorbar
%% Gas Layer results
%
% We plot the H2O mass fraction in the gas diffusion layer. We observe how the water is being first consumed close to
% the inlet (inlet is at bottom, outlet at top).
%
figure
plotCellData(model.(gs).grid, state.(gs).massfractions{1}, 'edgecolor', 'none');
title('H2O mass fraction / -')
colorbar
%%
% We plot the same result as a surface plot
%
N = gen.nxGasSupply;
X = reshape(model.(gs).grid.cells.centroids(:, 1), N, [])/(milli*meter);
Y = reshape(model.(gs).grid.cells.centroids(:, 2), N, [])/(milli*meter);
val = state.(gs).massfractions{1};
Z = reshape(val, N, []);
surf(X, Y, Z, 'edgecolor', 'none');
colorbar
title('Mass Fraction H2O');
xlabel('x / mm')
ylabel('y / mm')
view([50, 51]);
%%
%
% We plot the O2 mass fraction
%
figure('position', [1290, 755, 1275, 559])
val = state.(gs).massfractions{2};
Z = reshape(val, N, []);
surf(X, Y, Z, 'edgecolor', 'none');
colorbar
title('Mass Fraction O2');
xlabel('x / mm')
ylabel('y / mm')
view([50, 51]);
%% Interface results
%
% Current density in the anode. We recover the current on each face along the anode.
%
i = state.Electrolyser.Anode.i;
%%
%
% From the current values at each face, we compute the current density by dividing with the face areas.
ind = model.Electrolyser.couplingTerms{1}.couplingfaces(:, 2);
yc = model.Electrolyser.Electrolyte.grid.faces.centroids(ind, 2);
areas = model.Electrolyser.Electrolyte.grid.faces.areas(ind);
u = ampere/((centi*meter)^2);
i = (i./areas)/u; % We convert to A/cm^2
figure
plot(yc/(milli*meter), i);
title('Current in Anode / A/cm^2')
xlabel('height / mm')
%%
%
% We do the same but now for the proton current.
%
iHp = state.Electrolyser.Anode.iHp;
ind = model.Electrolyser.couplingTerms{1}.couplingfaces(:, 2);
yc = model.Electrolyser.Electrolyte.grid.faces.centroids(ind, 2);
areas = model.Electrolyser.Electrolyte.grid.faces.areas(ind);
u = ampere/((centi*meter)^2);
iHp = (iHp./areas)/u;
figure
plot(yc/(milli*meter), iHp);
title('iHp in Anode / A/cm^2')
xlabel('height / mm')
%%
%
% By dividing the proton current density with the total current density, we obtain the Faradic efficiency.
% drivingForces.src = schedule.control.src;
% state = model.evalVarName(state, 'Electrolyser.Anode.iHp', {{'drivingForces', drivingForces}});
i = state.Electrolyser.Anode.i;
iHp = state.Electrolyser.Anode.iHp;
ind = model.Electrolyser.couplingTerms{1}.couplingfaces(:, 2);
yc = model.Electrolyser.Electrolyte.grid.faces.centroids(ind, 2);
figure
plot(yc/(milli*meter), iHp./i);
title('Faradic efficiency')
xlabel('height / mm')