deploy: 8050f33

BattMoTeam · Dec 18, 2024 · b3bc678 · b3bc678
1 parent a8f461f
commit b3bc678
Show file tree

Hide file tree

Showing 76 changed files with 844 additions and 92 deletions.
diff --git a/_images/gas-layer-illustration.png b/_images/gas-layer-illustration.png
diff --git a/_images/runGasSupply_01.png b/_images/runGasSupply_01.png
diff --git a/_images/runGasSupply_02.png b/_images/runGasSupply_02.png
diff --git a/_images/runGasSupply_03.png b/_images/runGasSupply_03.png
diff --git a/_images/runGasSupply_04.png b/_images/runGasSupply_04.png
diff --git a/_images/runGasSupply_05.png b/_images/runGasSupply_05.png
diff --git a/_images/runGasSupply_06.png b/_images/runGasSupply_06.png
diff --git a/_images/runGasSupply_07.png b/_images/runGasSupply_07.png
diff --git a/_images/runGasSupply_08.png b/_images/runGasSupply_08.png
diff --git a/_images/runGasSupply_09.png b/_images/runGasSupply_09.png
diff --git a/_images/runProtonicMembrane_01.png b/_images/runProtonicMembrane_01.png
diff --git a/_images/runProtonicMembrane_02.png b/_images/runProtonicMembrane_02.png
diff --git a/_images/runProtonicMembrane_03.png b/_images/runProtonicMembrane_03.png
diff --git a/_images/runProtonicMembrane_04.png b/_images/runProtonicMembrane_04.png
diff --git a/_images/runProtonicMembrane_05.png b/_images/runProtonicMembrane_05.png
diff --git a/_images/runProtonicMembrane_06.png b/_images/runProtonicMembrane_06.png
diff --git a/_images/runProtonicMembrane_07.png b/_images/runProtonicMembrane_07.png
diff --git a/_sources/publishedExamples/runGasSupply.rst.txt b/_sources/publishedExamples/runGasSupply.rst.txt
@@ -6,8 +6,8 @@ Gas Supply Layer
 ================
 *Generated from runGasSupply.m*
 
+.. include:: runGasSupplyPreamble.rst
 
-Single gas supply layer where a mixture of two gases (here H2O and O2) is injected at the bottom and flows to the top (See illustration below).
 
 json input data
 ===============
@@ -177,7 +177,7 @@ We plot the water mass fractions at three different time step
 
 .. code-block:: matlab
 
-  inds = [1; 5; 10];
+  inds = [1; 2; 4; 8];
   ninds = numel(inds);
   for i = 1 : ninds
       figure
@@ -197,6 +197,9 @@ We plot the water mass fractions at three different time step
 .. figure:: runGasSupply_08.png
   :figwidth: 100%
 
+.. figure:: runGasSupply_09.png
+  :figwidth: 100%
+
 
 
 complete source code can be found :ref:`here<runGasSupply_source>`
diff --git a/_sources/publishedExamples/runGasSupplyPreamble.rst.txt b/_sources/publishedExamples/runGasSupplyPreamble.rst.txt
@@ -0,0 +1,92 @@
+Model overview
+==============
+
+We present here a simulation model for a single gas layer. The complete electroyser has two such layers, see
+illustration below.
+
+.. figure:: /img/gas-layer-illustration.png  
+   :target: ../_images/gas-layer-illustration.png
+   :align: center
+
+   Illustration of the electrolyser with the two gas layer domains
+
+We consider the gas layer alone, without coupling with the electrodes. We will choose here the gas layer used at the
+anode. We inject a mixture of two gases (H2O and O2). We model the gas supply layer as a porous media. The flow is then driven by a pressure gradient and the effect of diffusion.
+
+Governing equations
+-------------------
+
+The governing equations are given by mass conservation. We have
+
+.. math::
+   \frac{\partial\rho x_i}{\partial t} + \nabla\cdot\left(j_{\text{conv},i} + j_{\text{diff}, i}\right) = 0.
+
+where :math:`\rho` denotes the density, :math:`x_i` the mass fraction of the component :math:`i`, for
+:math:`i=\{1,2\}`. The convection and diffusion fluxes are denoted with :math:`j_{\text{conv},i}` and
+:math:`j_{\text{diff}, i}`, respectively.
+
+We use a simple equation of state for the gases, as we assume that they behave as ideal gas.
+
+For the momentum equation, we use the Darcy approximation which covers the case of a flow in a porous media. The Darcy velocity is given by
+
+.. math::
+   u = -\frac{K}{\mu}\nabla p,
+
+where :math:`K` denotes the permeability and :math:`mu` the viscosity. We are working on a model which uses the Stokes equations, which can be used for a free flow.
+
+The convection term for each species is given by
+
+.. math::
+   u_{\text{conv}, i} =  \rho x_i u.
+
+For the diffusion term, we use a Fick diffusion model and the diffusion flux for each component is given by
+
+.. math::
+   \hat j_{\text{diff}, i} = -D_i \nabla c_i
+
+where :math:`c_i` is the concentration of the specie :math:`i` and :math:`\hat j_{\text{diff}, i}` is the flux in
+:math:`mol/m^2/s`. We convert this flux in a mass flux, using the assumption that the gas are ideal gas, to obtain
+
+.. math::
+   j_{\text{diff}, i} = -D_i \nabla (\rho x_i)
+
+Our governing equations can now be rewritten using the explicit expression for the fluxes,
+
+.. math::
+   \frac{\partial\rho x_i}{\partial t} + \nabla\cdot\left(\rho x_i u - D_i\nabla (\rho x_i)\right) = 0
+
+This system of equation is discretized using finite volume and solved in time using backward Euler, for stability. The
+unknowns are the pressure and one of the mass fractions.
+
+Boundary conditions
+-------------------
+
+We have implemented boundary equations of three types, which can be used for both inlet and outlet:
+
+* **total mass rate and composition**
+* **pressure and composition**
+* **pressure and zero mass fraction derivative condition**
+
+The boundary conditions are applied on a subset of external faces in the grid (see drawing below). The two first types
+of boundary conditions are standard.
+
+When we impose the **rate and composition**, the given rate corresponds to the sum of all the input mass rates of every face
+which belongs to the boundary condition (It means that each boundary face individually can see different mass rates). A
+given composition corresponds to given mass fractions. This condition is use for the inlet in the example below, see
+json input :battmofile:`here <ProtonicMembrane/jsonfiles/gas-supply.json#16>`.
+
+The **pressure and composition** boundary condition corresponds to a standard Dirichlet condition. We would typically use
+such condition at the outlet. Even for a convection dominated flow, the composition at the boundary is then given. This
+could then result in an artifact in the outlet region because, in an experimental setup, we do not control the
+composition directly at the outlet, and we do not want to model and simulate the flow outside the gas layer
+region. Therefore, we introduced the third boundary condition.
+
+When we use a *zero mass fraction derivative* boundary condition, we impose that the *derivative* of the mass fraction
+in the normal direction at the boundary is zero, :math:`\frac{\partial x_i}{\partial n} = 0`. In this way, we model a
+situation where the outlet is connected to a long pipe which is itself connected to the ambient room, where we have a
+given composition. In the pipe, the flow is dominated by convection and we cannot observe the separation effect of
+diffusion (due to the difference in the value of the diffusion coefficient for each specie) at the beginning of the pipe
+. Hence, the vanishing condition of the derivative of the mass fractions. This condition is implemented in the example
+below, see :battmofile:`here <ProtonicMembrane/jsonfiles/gas-supply.json#19>`.
+
+
diff --git a/_sources/publishedExamples/runGasSupply_source.rst.txt b/_sources/publishedExamples/runGasSupply_source.rst.txt
@@ -10,9 +10,6 @@ Source code for runGasSupply
 
   %% Gas Supply Layer
   %
-  % Single gas supply layer where a mixture of two gases (here H2O and O2) is injected at the bottom and flows to the top
-  % (See illustration below).
-  %
   
   %% json input data
   %
@@ -144,7 +141,7 @@ Source code for runGasSupply
   %%
   % We plot the water mass fractions at three different time step
   
-  inds = [1; 5; 10];
+  inds = [1; 2; 4; 8];
   ninds = numel(inds);
   for i = 1 : ninds
       figure

diff --git a/_sources/publishedExamples/runProtonicMembrane.rst.txt b/_sources/publishedExamples/runProtonicMembrane.rst.txt
@@ -11,7 +11,7 @@ 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>`
+The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
 
 .. code-block:: matlab
 
@@ -31,8 +31,11 @@ We setup the input parameter structure which will we be used to instantiate the
 .. code-block:: matlab
 
   inputparams = ProtonicMembraneCellInputParams(jsonstruct);
-  
-  % We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
+
+We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
+
+.. code-block:: matlab
+
   [inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);
 
 
@@ -62,19 +65,13 @@ We setup the initial state using a default setup included in the model
 
 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).
+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 :battmofile:`here<ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
+We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
 
 .. code-block:: matlab
 
   schedule = model.Control.setupSchedule(inputparams.jsonstruct);
 
-We change the default tolerance
-
-.. code-block:: matlab
-
-  model.nonlinearTolerance = 1e-8;
-
 
 Simulation
 ==========
@@ -87,7 +84,7 @@ We run the simulation
 
 Plotting
 ========
-We setup som shortcuts for convenience
+We setup som shortcuts for convenience and introduce plotting options
 
 .. code-block:: matlab
 
@@ -105,7 +102,7 @@ We recover the position of the mesh cell of the discretization grid. This is use
 
   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.
+We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second line adds to the state variable all the variables that are derived from our primary unknowns.
 
 .. code-block:: matlab
 

diff --git a/_sources/publishedExamples/runProtonicMembranePreamble.rst.txt b/_sources/publishedExamples/runProtonicMembranePreamble.rst.txt
@@ -6,21 +6,23 @@ We consider the model of a mixed proton and electron conducting membrane, as des
 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
+In the membrane, we have three components or *species* given by the proton (:math:`H^+`), and the :math:`p` and :math:`n`
+type charge carriers. We need expressions for the fluxes for each of those. The governing equations will then be given
+by charge and mass conservation equations.
 
 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
+   \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
+   j_{\alpha} = -k_\alpha\nabla\bar\mu_\alpha,
 
 for some coefficient :math:`k_\alpha` which is not necessarily a constant.
 
@@ -43,18 +45,17 @@ For the :math:`p` and :math:`n` type conductivities,  we use the empirical relat
 
    \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).   
+Here :math:`E_{\text{ref},p}` and :math:`E_{\text{ref},n}` are two reference potentials. 
 
 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
+The **mass conservation equation** for :math:`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
+The total current density is given by :math:`i = F (j_{H^+} + j_p - j_n)` and the **charge conservation equation** is
 
 .. math::
 
@@ -72,7 +73,8 @@ We can rewrite :math:`i_{\text{el}}` as
    
    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
+We finally obtain the governing equations for :math:`\phi(x)` and :math:`\pi(x)` as the following system of differential
+equations
 
 .. math::
 
@@ -91,17 +93,18 @@ We define the over-potential :math:`\eta` as
    \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).
+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 
+Here, :math:`i_{l,c}` and :math:`i_{l,a}` are given constants. The value of the reference current density :math:`i_0` is
+also constant.
 
-The total current is given by :math:`i_{\text{an}} = I` for some constant current $I$.
+The total current is given by :math:`i_{\text{an}} = I` for some constant current :math:`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,

diff --git a/_sources/publishedExamples/runProtonicMembrane_source.rst.txt b/_sources/publishedExamples/runProtonicMembrane_source.rst.txt
@@ -11,7 +11,7 @@ Source code for runProtonicMembrane
   %% 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
+  % The source of the json files can be seen in :battmofile:`protonicMembrane.json<ProtonicMembrane/jsonfiles/protonicMembrane.json>` and
   % :battmofile:`1d-PM-geometry.json<ProtonicMembrane/jsonfiles/1d-PM-geometry.json>`
   
   filename = fullfile(battmoDir(), 'ProtonicMembrane', 'jsonfiles', 'protonicMembrane.json');
@@ -27,6 +27,7 @@ Source code for runProtonicMembrane
   
   inputparams = ProtonicMembraneCellInputParams(jsonstruct);
   
+  %%
   % We setup the grid, which is done by calling the function :battmo:`setupProtonicMembraneCellGrid`
   [inputparams, gen] = setupProtonicMembraneCellGrid(inputparams, jsonstruct);
   
@@ -44,17 +45,14 @@ Source code for runProtonicMembrane
   
   %% 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>`).
+  % control and the current equal to zero (see :battmofile:`here<ProtonicMembrane/jsonfiles/protonicMembrane.json#118>`).
   %
-  % 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).
+  % We compute the steady-state solution and the time stepping here does not correspond to time values but should be seen
+  % as step-wise increase of the effect of the non-linearity (in particular in the expression of the conductivity which
+  % includes highly nonlineaer effect with the exponential terms. We do not detail here the method).
   
   schedule = model.Control.setupSchedule(inputparams.jsonstruct);
   
-  %%
-  % We change the default tolerance
-  model.nonlinearTolerance = 1e-8;
-  
   %% Simulation
   % We run the simulation
   
@@ -64,7 +62,7 @@ Source code for runProtonicMembrane
   %
   
   %%
-  % We setup som shortcuts for convenience
+  % We setup som shortcuts for convenience and introduce plotting options
   an    = 'Anode';
   ct    = 'Cathode';
   elyte = 'Electrolyte';
@@ -79,7 +77,8 @@ Source code for runProtonicMembrane
   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.
+  % We consider the solution obtained at the last time step, which corresponds to the solution at steady-state. The second
+  % line adds to the state variable all the variables that are derived from our primary unknowns.
   state = states{end};
   state = model.addVariables(state, schedule.control);
   

diff --git a/advancedtopics.html b/advancedtopics.html
@@ -179,6 +179,11 @@
 </ul>
 </li>
 <li class="toctree-l1"><a class="reference internal" href="publishedExamples/runGasSupply.html">Gas Supply</a><ul>
+<li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-overview">Model overview</a><ul>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#governing-equations">Governing equations</a></li>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#boundary-conditions">Boundary conditions</a></li>
+</ul>
+</li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#json-input-data">json input data</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#input-parameter-setup">Input parameter setup</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-setup">Model setup</a></li>

diff --git a/app.html b/app.html
@@ -179,6 +179,11 @@
 </ul>
 </li>
 <li class="toctree-l1"><a class="reference internal" href="publishedExamples/runGasSupply.html">Gas Supply</a><ul>
+<li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-overview">Model overview</a><ul>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#governing-equations">Governing equations</a></li>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#boundary-conditions">Boundary conditions</a></li>
+</ul>
+</li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#json-input-data">json input data</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#input-parameter-setup">Input parameter setup</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-setup">Model setup</a></li>

diff --git a/app_calculations.html b/app_calculations.html
@@ -178,6 +178,11 @@
 </ul>
 </li>
 <li class="toctree-l1"><a class="reference internal" href="publishedExamples/runGasSupply.html">Gas Supply</a><ul>
+<li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-overview">Model overview</a><ul>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#governing-equations">Governing equations</a></li>
+<li class="toctree-l3"><a class="reference internal" href="publishedExamples/runGasSupply.html#boundary-conditions">Boundary conditions</a></li>
+</ul>
+</li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#json-input-data">json input data</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#input-parameter-setup">Input parameter setup</a></li>
 <li class="toctree-l2"><a class="reference internal" href="publishedExamples/runGasSupply.html#model-setup">Model setup</a></li>