deploy: c5960bf

_sources/publishedExamples/runProtonicMembrane.rst.txt

@@ -0,0 +1,153 @@
+.. _runProtonicMembrane:
+
+
+=======================
+Protonic Membrane model
+=======================
+*Generated from runProtonicMembrane.m*
+
+.. include:: runProtonicMembranePreamble.rst
+
+
+Load and parse input from given json files
+==========================================
+The source of the json files can be seen in :battmofile:`protonicMembrane<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
+
+.. code-block:: matlab
+
+  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
+
+.. code-block:: matlab
+
+  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
+
+.. code-block:: matlab
+
+  model = ProtonicMembraneCell(inputparams);
+
+The model is equipped for simulation using the following command (this step may become unnecessary in future versions)
+
+.. code-block:: matlab
+
+  model = model.setupForSimulation();
+
+
+Initial state setup
+===================
+We setup the initial state using a default setup included in the model
+
+.. code-block:: matlab
+
+  state0 = model.setupInitialState();
+
+
+Schedule 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 :battmofile:`here<ProtonicMembrane/protonicMembrane.json#86>`).
+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).
+
+.. code-block:: matlab
+
+  schedule = model.Control.setupSchedule(inputparams.jsonstruct);
+
+We change the default tolerance
+
+.. code-block:: matlab
+
+  model.nonlinearTolerance = 1e-8;
+
+
+Simulation
+==========
+We run the simulation
+
+.. code-block:: matlab
+
+  [~, states, report] = simulateScheduleAD(state0, model, schedule);
+
+
+Plotting
+========
+We setup som shortcuts for convenience
+
+.. code-block:: matlab
+
+  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.
+
+.. code-block:: matlab
+
+  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.
+
+.. code-block:: matlab
+
+  state = states{end};
+  state = model.addVariables(state, schedule.control);
+  
+  figure(1)
+  plot(xc, state.(elyte).pi)
+  title('Electromotive potential (\pi)')
+  xlabel('x / m')
+  ylabel('\pi / V')
+  
+  figure(2)
+  plot(xc, state.(elyte).pi - state.(elyte).phi)
+  title('Electronic chemical potential (E)')
+  xlabel('x / m')
+  ylabel('E / V')
+  
+  figure(3)
+  plot(xc, state.(elyte).phi)
+  title('Electrostatic potential (\phi)')
+  xlabel('x / m')
+  ylabel('\phi / V')
+  
+  figure(4)
+  plot(xc, log(state.(elyte).sigmaEl))
+  title('logarithm of conductivity (\sigma)')
+  xlabel('x / m')
+  xlabel('log(\sigma/Siemens)')
+
+.. figure:: runProtonicMembrane_01.png
+  :figwidth: 100%
+
+.. figure:: runProtonicMembrane_02.png
+  :figwidth: 100%
+
+.. figure:: runProtonicMembrane_03.png
+  :figwidth: 100%
+
+.. figure:: runProtonicMembrane_04.png
+  :figwidth: 100%
+
+
+
+complete source code can be found :ref:`here<runProtonicMembrane_source>`

_sources/publishedExamples/runProtonicMembranePreamble.rst.txt

@@ -0,0 +1,114 @@
+Model overview
+==============
+
+We consider the model of a mixed proton and electron conducting membrane, as described in :cite:`V_llestad_2019` from 2019.
+
+Governing equations
+-------------------
+
+We have three components given by the proton (:math:`H^+`), the :math:`p` and :math:`n` type. We want to find expressions for the fluxes for each of the components
+
+We denote by :math:`\phi` the electrostatic potential. For each of the components :math:`\alpha=\{H^+, p, n\}`, we
+introduce the electrochemical potential denoted :math:`\bar\mu_\alpha` and the chemical potential denoted
+:math:`\mu_\alpha`. The potentials are related with each other through the relation
+
+.. math::
+
+   \bar\mu_\alpha = \mu_\alpha + z_\alpha F \phi
+
+The fluxes are governed by the gradient of the electrochemical potential. We have
+
+.. math::
+   
+   j_{\alpha} = -k_\alpha\nabla\bar\mu_\alpha
+
+for some coefficient :math:`k_\alpha` which is not necessarily a constant.
+
+We assume that the chemical potential of :math:`H^+`, that is :math:`\mu_{H^+}`, is constant in the electrolyte and that
+:math:`k_{H^+}` is constant. Hence, the current density of :math:`H^+` is given
+
+.. math::          
+
+    i_{H^+} = -\sigma_{H^+} \nabla \phi,
+
+for a constant conductivity :math:`\sigma_{H^+}`. We denote by :math:`E = -\frac{\mu_n}{F}` the electronic chemical potential and :math:`\pi = E + \phi` the electromotive potential. Since we have equilibrium for the n-p reaction, we have the identity
+
+.. math::
+
+   d\mu_p + d\mu_n = 0.
+
+For the :math:`p` and :math:`n` type conductivities,  we use the empirical relations
+
+.. math::
+
+   \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 :math:`E_{\text{ref},p}` and :math:`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 :math:`\phi(x)` and :math:`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
+
+.. math::
+
+   \nabla\cdot j_{H^+} = 0.
+
+The current density is given by :math:`i = F (j_{H^+} + j_p - j_n)` and the \textbf{charge
+conservation equation} is
+
+.. math::
+
+   \nabla\cdot i = 0.
+   
+We introduce the electronic current density :math:`i_{\text{el}} = F(j_p - j_n)` and we get
+
+.. math::
+   
+   \nabla\cdot i_{\text{el}} = 0.
+
+We can rewrite :math:`i_{\text{el}}` as
+
+.. math::
+   
+   i_{\text{el}} = - (\sigma_p(E) + \sigma_n(E))\nabla ( E + \phi ).
+
+The governing equations for :math:`\phi(x)` and :math:`\pi(x)` are therefore the differential equations
+
+.. math::
+
+   \begin{align}
+   -\nabla\cdot(\sigma_{H^+}\nabla\phi) &= 0,\\
+   -\nabla\cdot((\sigma_p(E) + \sigma_n(E))\nabla \pi) &= 0.
+   \end{align}   
+
+Boundary conditions
+-------------------
+
+We define the over-potential :math:`\eta` as
+
+.. math::
+    
+   \eta_\text{elde} = \pi_\text{elde} - \phi_\text{elde} - \text{OCP}_\text{elde}
+
+where :math:`\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 :cite:`V_llestad_2019` for the expressions).
+
+At the anode, we imposte that the proton current  is given through the following Buttler-Volmer type expression
+
+.. math::
+   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, :math:`i_{l,c}` and :math:`i_{l,a}` are given constants. The value of the reference current density $i_0$ is also constant 
+
+The total current is given by :math:`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 :math:`H^+` current and the over-potential that takes a linear form,
+
+.. math::
+   
+   R_\text{CT} \ i_{H^+, \text{ct}} = \eta_\text{ct},
+
+for a given charge transfer constant :math:`R_\text{CT}`.
+

_sources/publishedExamples/runProtonicMembrane_source.rst.txt

@@ -0,0 +1,115 @@
+:orphan:
+
+.. _runProtonicMembrane_source:
+
+Source code for runProtonicMembrane
+-----------------------------------
+
+.. code:: matlab
+
+
+  %% Protonic Membrane model
+  
+  
+  %% Load and parse input from given json files
+  % The source of the json files can be seen in :battmofile:`protonicMembrane<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
+  % :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/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 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 :battmofile:`here<ProtonicMembrane/protonicMembrane.json#86>`).
+  %
+  % 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);
+  
+  figure(1)
+  plot(xc, state.(elyte).pi)
+  title('Electromotive potential (\pi)')
+  xlabel('x / m')
+  ylabel('\pi / V')
+  
+  figure(2)
+  plot(xc, state.(elyte).pi - state.(elyte).phi)
+  title('Electronic chemical potential (E)')
+  xlabel('x / m')
+  ylabel('E / V')
+  
+  figure(3)
+  plot(xc, state.(elyte).phi)
+  title('Electrostatic potential (\phi)')
+  xlabel('x / m')
+  ylabel('\phi / V')
+  
+  figure(4)
+  plot(xc, log(state.(elyte).sigmaEl))
+  title('logarithm of conductivity (\sigma)')
+  xlabel('x / m')
+  xlabel('log(\sigma/Siemens)')
+  
+  
+  
+  
+  
+