Protonic Membrane model

Generated from runProtonicMembrane.m

Model overview

We consider the model of a mixed proton and electron conducting membrane, as described in [VST+19] from 2019.

Governing equations

We have three components given by the proton ($H^+$), the $p$ and $n$ type. We want to find expressions for the fluxes for each of the components

We denote by $\phi$ the electrostatic potential. For each of the components $\alpha=\{H^+, p, n\}$, we introduce the electrochemical potential denoted $\bar\mu_\alpha$ and the chemical potential denoted $\mu_\alpha$. The potentials are related with each other through the relation

\[\bar\mu_\alpha = \mu_\alpha + z_\alpha F \phi\]

The fluxes are governed by the gradient of the electrochemical potential. We have

\[j_{\alpha} = -k_\alpha\nabla\bar\mu_\alpha\]

for some coefficient $k_\alpha$ which is not necessarily a constant.

We assume that the chemical potential of $H^+$, that is $\mu_{H^+}$, is constant in the electrolyte and that $k_{H^+}$ is constant. Hence, the current density of $H^+$ is given

\[i_{H^+} = -\sigma_{H^+} \nabla \phi,\]

for a constant conductivity $\sigma_{H^+}$. We denote by $E = -\frac{\mu_n}{F}$ the electronic chemical potential and $\pi = E + \phi$ the electromotive potential. Since we have equilibrium for the n-p reaction, we have the identity

\[d\mu_p + d\mu_n = 0.\]

For the $p$ and $n$ type conductivities, we use the empirical relations

\[\sigma_p(E) = \sigma_p^0\exp\left(\frac{F(E - E_{\text{ref},p})}{RT}\right)\quad\text{ and }\quad\sigma_n(E) = \sigma_n^0\exp\left(\frac{-F(E - E_{\text{ref},n})}{RT}\right).\]

Here $E_{\text{ref},p}$ and $E_{\text{ref},n}$ are two reference potentials (see below).

We consider the steady state. We could introduce later charge and mass capacitors. The unknowns are the functions $\phi(x)$ and $E(x)$ in the electrolyte. The governing equations are given by the mass conservation for the proton and the charge conservation.

The mass conservation for $H^+$ is given by

\[\nabla\cdot j_{H^+} = 0.\]

The current density is given by $i = F (j_{H^+} + j_p - j_n)$ and the textbf{charge conservation equation} is

\[\nabla\cdot i = 0.\]

We introduce the electronic current density $i_{\text{el}} = F(j_p - j_n)$ and we get

\[\nabla\cdot i_{\text{el}} = 0.\]

We can rewrite $i_{\text{el}}$ as

\[i_{\text{el}} = - (\sigma_p(E) + \sigma_n(E))\nabla ( E + \phi ).\]

The governing equations for $\phi(x)$ and $\pi(x)$ are therefore the differential equations

\[\begin{split}\begin{align} -\nabla\cdot(\sigma_{H^+}\nabla\phi) &= 0,\\ -\nabla\cdot((\sigma_p(E) + \sigma_n(E))\nabla \pi) &= 0. \end{align}\end{split}\]

Boundary conditions

We define the over-potential $\eta$ as

\[\eta_\text{elde} = \pi_\text{elde} - \phi_\text{elde} - \text{OCP}_\text{elde}\]

where $\text{OCP}_\text{elde}$ is the open-circuit potential for the given electrode. The value of the $text{OCP}$ at each electrode depend on the composition at the electrode (see [VST+19] for the expressions).

At the anode, we imposte that the proton current is given through the following Buttler-Volmer type expression

\[i_{H^+, \text{an}} = -i_0\frac{e^{-\beta\frac{ z F \eta_{\text{an}}}{RT}} - e^{( 1- \beta)\frac{ z F \eta_{\text{an}}}{RT}}}{ 1+ \frac{i_0}{i_{l,c}}e^{-\beta\frac{ z F \eta_{\text{an}}}{RT}} - \frac{i_0}{i_{l,a}}e^{( 1- \beta)\frac{ z F \eta_{\text{an}}}{RT}}}\]

Here, $i_{l,c}$ and $i_{l,a}$ are given constants. The value of the reference current density $i_0$ is also constant

The total current is given by $i_{\text{an}} = I$ for some constant current $I$.

At the cathode, we impose that the electrostatic potential is equal to zero and a relation between the $H^+$ current and the over-potential that takes a linear form,

\[R_\text{CT} \ i_{H^+, \text{ct}} = \eta_\text{ct},\]

for a given charge transfer constant $R_\text{CT}$.

Load and parse input from given json files

The source of the json files can be seen in protonicMembrane and 1d-PM-geometry.json

filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', 'protonicMembrane.json');
jsonstruct_material = parseBattmoJson(filename);

filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', '1d-PM-geometry.json');
jsonstruct_geometry = parseBattmoJson(filename);

jsonstruct = mergeJsonStructs({jsonstruct_material, jsonstruct_geometry});

Input structure setup

We setup the input parameter structure which will we be used to instantiate the model

inputparams = ProtonicMembraneCellInputParams(jsonstruct);

% We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
[inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);

Model setup

We instantiate the model for the protonic membrane cell

model = ProtonicMembraneCell(inputparams);

The model is equipped for simulation using the following command (this step may become unnecessary in future versions)

model = model.setupForSimulation();

Initial state setup

We setup the initial state using a default setup included in the model

state0 = model.setupInitialState();

Schedule

We setup the schedule, which means the timesteps and also the control we want to use. In this case we use current control and the current equal to zero (see here here). We compute the steady-state solution so that the time stepping here is more an artifact to reach the steady-state solution. In particular, it governs the pace at which we increase the non-linearity (not detailed here).

schedule = model.Control.setupSchedule(inputparams.jsonstruct);

We change the default tolerance

model.nonlinearTolerance = 1e-8;

Simulation

We run the simulation

[~, states, report] = simulateScheduleAD(state0, model, schedule);

Plotting

We setup som shortcuts for convenience

an    = 'Anode';
ct    = 'Cathode';
elyte = 'Electrolyte';
ctrl  = 'Control';

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

We recover the position of the mesh cell of the discretization grid. This is used when plotting the spatial distribution of some of the variables.

xc = model.(elyte).grid.cells.centroids(:, 1);

We consider the solution obtained at the last time step, which corresponds to the solution at steady-state.

state = states{end};
state = model.addVariables(state, schedule.control);

Plot of electromotive potential

figure
plot(xc, state.(elyte).pi)
title('Electromotive potential (\pi)')
xlabel('x / m')
ylabel('\pi / V')

Plot of electronic chemical potential

figure
plot(xc, state.(elyte).pi - state.(elyte).phi)
title('Electronic chemical potential (E)')
xlabel('x / m')
ylabel('E / V')

Plot of electrostatic potential

figure
plot(xc, state.(elyte).phi)
title('Electrostatic potential (\phi)')
xlabel('x / m')
ylabel('\phi / V')

Plot of the conductivity

figure
plot(xc, log(state.(elyte).sigmaEl))
title('logarithm of conductivity (\sigma)')
xlabel('x / m')
xlabel('log(\sigma/Siemens)')

Evolution of the Faradic efficiency

We increase the current density from 0 to 1 A/cm^2 and plot the faraday efficiency. We sample the current value from 0 to 1 A/cm^2

Is = linspace(0, 1*ampere/((centi*meter)^2), 20);

We run the simulation for each current value and collect the results in the endstates

endstates = {};
for iI = 1 : numel(Is)

    model.Control.I = Is(iI);
    [~, states, report] = simulateScheduleAD(state0, model, schedule, 'NonLinearSolver', nls);

    state = states{end};
    state = model.addVariables(state, schedule.control);

    endstates{iI} = state;

end

We plot the profile of the electromotive potential for the mininum and maximum current values.

figure
hold on

unit = ampere/((centi*meter)^2); % shortcut

state = endstates{1};
plot(xc, state.(elyte).pi, 'displayname', sprintf('I=%g A/cm^2', Is(1)/unit));
state = endstates{end};
plot(xc, state.(elyte).pi, 'displayname', sprintf('I=%g A/cm^2', Is(end)/unit));

title('Electromotive potential (\pi)')
xlabel('x / m')
ylabel('\pi / V')
legend

We retrieve and plot the cell potential

E = cellfun(@(state) state.(an).pi - state.(ct).pi, endstates);

figure
plot(Is/unit, E, '*-');
xlabel('Current density / A/cm^2')
ylabel('E / V')
title('Cell potential');

We retrieve and plot the Faradic efficiency

feff = cellfun(@(state) state.(an).iHp/state.(an).i, endstates);

figure
plot(Is/unit, feff, '*-');
xlabel('Current density / A/cm^2')
ylabel('Faradic efficiency / -')
title('Faradic efficiency');

complete source code can be found here