Merge pull request #3116 from pybamm-team/msmr

MSMR model

CHANGELOG.md

@@ -1,11 +1,11 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
-
 - The parameter "Ambient temperature [K]" can now be given as a function of position `(y,z)` and time `t`. The "edge" and "current collector" heat transfer coefficient parameters can also depend on `(y,z)` ([#3257](https://github.com/pybamm-team/PyBaMM/pull/3257))
 - Spherical and cylindrical shell domains can now be solved with any boundary conditions ([#3237](https://github.com/pybamm-team/PyBaMM/pull/3237))
 - Processed variables now get the spatial variables automatically, allowing plotting of more generic models ([#3234](https://github.com/pybamm-team/PyBaMM/pull/3234))
 - Numpy functions now work with PyBaMM symbols (e.g. `np.exp(pybamm.Symbol("a"))` returns `pybamm.Exp(pybamm.Symbol("a"))`). This means that parameter functions can be specified using numpy functions instead of pybamm functions. Additionally, combining numpy arrays with pybamm objects now works (the numpy array is converted to a pybamm array) ([#3205](https://github.com/pybamm-team/PyBaMM/pull/3205))
+- Implement the MSMR model ([#3116](https://github.com/pybamm-team/PyBaMM/pull/3116))
 
 ## Bug fixes
 

docs/source/api/models/lithium_ion/electrode_soh.rst

@@ -8,4 +8,8 @@ Electrode SOH models
 
 .. autofunction:: pybamm.lithium_ion.get_min_max_stoichiometries
 
+.. autofunction:: pybamm.lithium_ion.get_initial_ocps
+
+.. autofunction:: pybamm.lithium_ion.get_min_max_ocps
+
 .. footbibliography::

docs/source/api/models/lithium_ion/index.rst

@@ -9,5 +9,6 @@ Lithium-ion Models
   mpm
   dfn
   newman_tobias
+  msmr
   yang2017
   electrode_soh

docs/source/api/models/lithium_ion/msmr.rst

@@ -0,0 +1,7 @@
+Multi-Species Multi-Reaction (MSMR) Model
+=========================================
+
+.. autoclass:: pybamm.lithium_ion.MSMR
+    :members:
+
+.. footbibliography::

docs/source/api/models/submodels/interface/kinetics/butler_volmer.rst

@@ -1,5 +1,5 @@
-Butler Volumer
-==============
+Butler Volmer
+=============
 
 .. autoclass:: pybamm.kinetics.SymmetricButlerVolmer
     :members:

docs/source/api/models/submodels/interface/kinetics/index.rst

@@ -10,5 +10,6 @@ Kinetics
   marcus
   no_reaction
   tafel
+  msmr_butler_volmer
   total_main_kinetics
   inverse_kinetics/index

docs/source/api/models/submodels/interface/kinetics/msmr_butler_volmer.rst

@@ -0,0 +1,5 @@
+MSMR Butler Volmer
+==================
+
+.. autoclass:: pybamm.kinetics.MSMRButlerVolmer
+    :members:

docs/source/api/models/submodels/interface/open_circuit_potential/index.rst

@@ -6,3 +6,4 @@ Open-circuit potential models
   base_ocp
   current_sigmoid_ocp
   single_ocp
+  msmr_ocp

docs/source/api/models/submodels/interface/open_circuit_potential/msmr_ocp.rst

@@ -0,0 +1,8 @@
+MSMR Open Circuit Potential
+===========================
+
+
+.. autoclass:: pybamm.open_circuit_potential.MSMROpenCircuitPotential
+    :members:
+
+.. footbibliography::

docs/source/api/models/submodels/particle/index.rst

@@ -8,3 +8,4 @@ Particle
   fickian_diffusion
   polynomial_profile
   x_averaged_polynomial_profile
+  msmr_diffusion

docs/source/api/models/submodels/particle/msmr_diffusion.rst

@@ -0,0 +1,7 @@
+MSMR Diffusion
+==============
+
+.. autoclass:: pybamm.particle.MSMRDiffusion
+    :members:
+
+.. footbibliography::

docs/source/examples/index.rst

@@ -59,6 +59,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
     notebooks/models/lead-acid.ipynb
     notebooks/models/lithium-plating.ipynb
     notebooks/models/MPM.ipynb
+    notebooks/models/MSMR.ipynb
     notebooks/models/pouch-cell-model.ipynb
     notebooks/models/rate-capability.ipynb
     notebooks/models/SEI-on-cracks.ipynb

docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb

@@ -199,6 +199,7 @@
     "V_hat = pybamm.Parameter(\"Partial molar volume [m3.mol-1]\")\n",
     "c_inf = pybamm.Parameter(\"Bulk electrolyte solvent concentration [mol.m-3]\")\n",
     "\n",
+    "\n",
     "def D(cc):\n",
     "    return pybamm.FunctionParameter(\"Diffusivity [m2.s-1]\", {\"Solvent concentration [mol.m-3]\": cc})"
    ]
@@ -485,9 +486,11 @@
     "    {\"SEI layer\": {xi: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}}\n",
     ")\n",
     "\n",
+    "\n",
     "def Diffusivity(cc):\n",
     "    return cc * 10**(-12)\n",
     "\n",
+    "\n",
     "# parameter values (not physically based, for example only!)\n",
     "param = pybamm.ParameterValues(\n",
     "    {\n",
@@ -565,6 +568,7 @@
     "L_0_eval = param.evaluate(L_0)\n",
     "xi = np.linspace(0, 1, 100) # dimensionless space\n",
     "\n",
+    "\n",
     "def plot(t):\n",
     "    _, (ax1, ax2) = plt.subplots(1, 2 ,figsize=(10,5))\n",
     "    ax1.plot(solution.t, L_out(solution.t) * 1e6)\n",
@@ -581,6 +585,7 @@
     "    plt.tight_layout()\n",
     "    plt.show()\n",
     "    \n",
+    "\n",
     "import ipywidgets as widgets\n",
     "widgets.interact(plot, t=widgets.FloatSlider(min=0,max=solution.t[-1],step=0.1,value=0));"
    ]

docs/source/examples/notebooks/models/DFN.ipynb

@@ -62,22 +62,22 @@
     "\n",
     "#### Current:\n",
     "$$\n",
-    "i_{\\text{e,n}}\\big|_{x=0} = 0, \\quad i_{\\text{e,p}}\\big|_{x=1}=0,  \\\\\n",
+    "i_{\\text{e,n}}\\big|_{x=0} = 0, \\quad i_{\\text{e,p}}\\big|_{x=L}=0,  \\\\\n",
     "\\phi_{\\text{e,n}}\\big|_{x=L_{\\text{n}}} = \\phi_{\\text{e,s}}\\big|_{x=L_{\\text{n}}}, \\quad i_{\\text{e,n}}\\big|_{x=L_{\\text{n}}} = i_{\\text{e,s}}\\big\\vert_{x=L_{\\text{n}}} = I, \\\\ \n",
-    "\\phi_{\\text{e,s}}\\big|_{x=1-L_{\\text{p}}} = \\phi_{\\text{e,p}}\\big|_{x=1-L_{\\text{p}}}, \\quad \n",
-    "    i_{\\text{e,s}}\\big|_{x=1-L_{\\text{p}}} = i_{\\text{e,p}}\\big|_{x=1-L_{\\text{p}}} = I.\n",
+    "\\phi_{\\text{e,s}}\\big|_{x=L-L_{\\text{p}}} = \\phi_{\\text{e,p}}\\big|_{x=L-L_{\\text{p}}}, \\quad \n",
+    "    i_{\\text{e,s}}\\big|_{x=L-L_{\\text{p}}} = i_{\\text{e,p}}\\big|_{x=L-L_{\\text{p}}} = I.\n",
     "$$\n",
     "\n",
     "#### Concentration in the electrolyte:\n",
     "$$\n",
-    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=1}=0,\\\\ \n",
+    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=L}=0,\\\\ \n",
     "c_{\\text{e,n}}\\big|_{x=L_{\\text{n}}} = c_{\\text{e,s}}|_{x=L_{\\text{n}}}, \\quad N_{\\text{e,n}}\\big|_{x=L_{\\text{n}}}=N_{\\text{e,s}}\\big|_{x=L_{\\text{n}}}, \\\\\n",
-    "c_{\\text{e,s}}|_{x=1-L_{\\text{p}}}=c_{\\text{e,p}}|_{x=1-L_{\\text{p}}}, \\quad N_{\\text{e,s}}\\big|_{x=1-L_{\\text{p}}}=N_{\\text{e,p}}\\big|_{x=1-L_{\\text{p}}}.\n",
+    "c_{\\text{e,s}}|_{x=L-L_{\\text{p}}}=c_{\\text{e,p}}|_{x=L-L_{\\text{p}}}, \\quad N_{\\text{e,s}}\\big|_{x=L-L_{\\text{p}}}=N_{\\text{e,p}}\\big|_{x=L-L_{\\text{p}}}.\n",
     "$$\n",
     "\n",
     "####  Concentration in the electrode active material:\n",
     "$$\n",
-    "N_{\\text{s,k}}\\big|_{r_{\\text{k}}=0} = 0, \\quad \\text{k} \\in \\text{n, p}, \\quad \\ \\ - N_{\\text{s,k}}\\big|_{r_{\\text{k}}=1} = \\frac{j_{\\text{k}}}{F}, \\quad \\text{k} \\in \\text{n, p}.\n",
+    "N_{\\text{s,k}}\\big|_{r_{\\text{k}}=0} = 0, \\quad \\text{k} \\in \\text{n, p}, \\quad \\ \\ N_{\\text{s,k}}\\big|_{r_{\\text{k}}=R_{\\text{k}}} = \\frac{j_{\\text{k}}}{F}, \\quad \\text{k} \\in \\text{n, p}.\n",
     "$$\n",
     "\n",
     "#### Reference potential:\n",
@@ -87,8 +87,8 @@
     "#### And the initial conditions:\n",
     "\n",
     "$$\n",
-    "c_{\\text{s,k}}(x,r,0) = c_{\\text{s,k,0}}, \\quad \\phi_{\\text{s,n}}(x,0) = 0, \\quad \\phi_{\\text{s,p}}(x,0) = \\phi_{\\text{s,p,0}}, \\\\ \\text{k} \\in \\text{n, p},\\\\\n",
-    "\\phi_{\\text{e,k}}(x,0) = \\phi_{\\text{e,0}}, \\quad c_{\\text{e,k}}(x,0) = 1, \\\\ \\text{k} \\in \\text{n, s, p}. \n",
+    "c_{\\text{s,k}}(x,r,0) = c_{\\text{s,k,0}}, \\quad \\text{k} \\in \\text{n, p},\\\\\n",
+    "c_{\\text{e,k}}(x,0) = c_{\\text{e,0}}, \\quad \\text{k} \\in \\text{n, s, p}. \n",
     "$$\n"
    ]
   },
@@ -269,7 +269,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pybamm",
+   "display_name": "dev",
    "language": "python",
    "name": "python3"
   },
@@ -287,7 +287,7 @@
   },
   "vscode": {
    "interpreter": {
-    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
+    "hash": "bca2b99bfac80e18288b793d52fa0653ab9b5fe5d22e7b211c44eb982a41c00c"
    }
   }
  },

docs/source/examples/notebooks/models/SPMe.ipynb

@@ -25,28 +25,26 @@
     "\n",
     "ii) At the centre of each particle the standard no-flux condition is imposed, and the flux on the surface of the particle is simply the current $I$ divided by the thickness of the electrode $L_{\\text{k}}$, as in the SPM. Since lithium is transferred between the electrolyte and particles, the flux through the particle surface also enters the electrolyte diffusion equation as a source/sink term. There is no transfer of lithium between the electrolyte and current collectors, which leads to no flux boundary conditions on the lithium concentration in the electrolyte $c_{\\text{e,k}}$ at either end of the cell. \n",
     "\n",
-    "iii) We must also impose initial conditions which correspond to setting an initial concentration in each particle $c_{\\text{s,k}}(t=0) = c_{\\text{s,k,0}}$, and to having no deviation from the initial (uniform) lithium concentration in the electrolyte $c_{\\text{e,k}}(t=0) = 0$.  \n",
+    "iii) We must also impose initial conditions which correspond to setting an initial concentration in each particle $c_{\\text{s,k}}(t=0) = c_{\\text{s,k,0}}$, and to having no deviation from the initial (uniform) lithium concentration in the electrolyte $c_{\\text{e,k}}(t=0) =  c_{\\text{e,0}}$.  \n",
     "\n",
     "\n",
     "The model equations for the SPMe read: \n",
     "\n",
     "\n",
     "#### Particles: \n",
     "$$\n",
-    "\\mathcal{C}_{\\text{k}} \\frac{\\partial c_{\\text{s,k}}}{\\partial t} = -\\frac{1}{r_{\\text{k}}^2} \\frac{\\partial}{\\partial r_{\\text{k}}} \\left(r_{\\text{k}}^2 N_{\\text{s,k}}\\right), \\\\\n",
+    "\\frac{\\partial c_{\\text{s,k}}}{\\partial t} = -\\frac{1}{r_{\\text{k}}^2} \\frac{\\partial}{\\partial r_{\\text{k}}} \\left(r_{\\text{k}}^2 N_{\\text{s,k}}\\right), \\\\\n",
     "N_{\\text{s,k}} = -D_{\\text{s,k}}(c_{\\text{s,k}}) \\frac{\\partial c_{\\text{s,k}}}{\\partial r_{\\text{k}}}, \\quad \\text{k} \\in \\text{n, p},\n",
     "$$\n",
     "\n",
     "$$\n",
-    "N_{\\text{s,k}}\\big|_{r_{\\text{k}}=0} = 0, \\quad \\text{k} \\in \\text{n, p}, \\quad \\ \\ - \\frac{a_{R, \\text{k}}\\gamma_{\\text{k}}}{\\mathcal{C}_{\\text{k}}} N_{\\text{s,k}}\\big|_{r_{\\text{k}}=1} = \n",
+    "N_{\\text{s,k}}\\big|_{r_{\\text{k}}=0} = 0, \\quad \\text{k} \\in \\text{n, p}, \\quad \\ \\ N_{\\text{s,k}}\\big|_{r_{\\text{k}}=R_{\\text{k}}} = \n",
     "\\begin{cases}\n",
-    "\t\t  \\frac{I}{L_{\\text{n}}}, \\quad &\\text{k}=\\text{n}, \\\\ \n",
-    "\t\t  -\\frac{I}{L_{\\text{p}}}, \\quad &\\text{k}=\\text{p}, \n",
+    "\t\t  \\frac{I}{Fa_{\\text{n}}L_{\\text{n}}}, \\quad &\\text{k}=\\text{n}, \\\\ \n",
+    "\t\t  -\\frac{I}{Fa_{\\text{p}}L_{\\text{p}}}, \\quad &\\text{k}=\\text{p}, \n",
     "\\end{cases} \\\\\n",
-    "c_{\\text{s,k}}(r_{\\text{k}},0) = c_{\\text{s,k,0}}, \\quad \\text{k} \\in \\text{n, p},\n",
     "$$\n",
-    "\n",
-    "where $D_{\\text{s,k}}$ is the diffusion coefficient in the solid, $N_{\\text{s,k}}$ denotes the flux of lithium ions in the solid particle within the region $\\text{k}$, and $r_{\\text{k}} \\in[0,1]$ is the radial coordinate of the particle in electrode $\\text{k}$. All other relevant parameters are given in the table at the end of this notebook.\n",
+    "where $D_{\\text{s,k}}$ is the diffusion coefficient in the solid, $N_{\\text{s,k}}$ denotes the flux of lithium ions in the solid particle within the region $\\text{k}$, and $r_{\\text{k}} \\in[0,R_{\\text{k}}]$ is the radial coordinate of the particle in electrode $\\text{k}$. All other relevant parameters are given in the table at the end of this notebook.\n",
     "\n",
     "\n",
     "#### Electrolyte: \n",
@@ -66,10 +64,10 @@
     "$$\n",
     "\n",
     "$$\n",
-    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=1}=0, \\\\\n",
+    "N_{\\text{e,n}}\\big|_{x=0} = 0, \\quad N_{\\text{e,p}}\\big|_{x=L}=0, \\\\\n",
     "c_{\\text{e,k}}(x,0) = 0, \\quad \\text{k} \\in \\text{n, s, p},\n",
     "$$\n",
-    "where $D_{\\text{e}}$ is the diffusion coefficient in the solid, $N_{\\text{e,k}}$ denotes the flux of lithium ions in the electrolyte within the region $\\text{k}$, and $x\\in[0,1]$ is the macroscopic through-cell distance. This equation is also solved subject to continuity of concentration and flux at the electrode/separator interfaces.\n",
+    "where $D_{\\text{e}}$ is the diffusion coefficient in the solid, $N_{\\text{e,k}}$ denotes the flux of lithium ions in the electrolyte within the region $\\text{k}$, and $x\\in[0,L]$ is the macroscopic through-cell distance. This equation is also solved subject to continuity of concentration and flux at the electrode/separator interfaces.\n",
     "\n",
     "### Voltage Expression\n",
     "The voltage is obtained from the expression: \n",
@@ -90,7 +88,7 @@
     "where\n",
     "$$\n",
     "\\bar{c}_{\\text{e,n}} = \\frac{1}{L_{\\text{n}}}\\int_0^{L_{\\text{n}}} c_{\\text{e,n}} \\, \\text{d}x, \\quad\n",
-    "\\bar{c}_{\\text{e,p}} = \\frac{1}{L_{\\text{p}}}\\int_{1-L_{\\text{p}}}^{1} c_{\\text{e,p}} \\, \\text{d}x.\n",
+    "\\bar{c}_{\\text{e,p}} = \\frac{1}{L_{\\text{p}}}\\int_{L-L_{\\text{p}}}^{L} c_{\\text{e,p}} \\, \\text{d}x.\n",
     "$$\n",
     "\n",
     "More details can be found in [[3]](#References)."
@@ -248,7 +246,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pybamm",
+   "display_name": "dev",
    "language": "python",
    "name": "python3"
   },
@@ -266,7 +264,7 @@
   },
   "vscode": {
    "interpreter": {
-    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
+    "hash": "bca2b99bfac80e18288b793d52fa0653ab9b5fe5d22e7b211c44eb982a41c00c"
    }
   }
  },

docs/source/examples/notebooks/models/jelly-roll-model.ipynb

@@ -280,23 +280,36 @@
    "outputs": [],
    "source": [
     "# define spiral \n",
+    "\n",
+    "\n",
     "def spiral_pos_inner(t):\n",
     "    return r0 - eps * delta + eps * t / (2 * pi)\n",
+    "\n",
+    "\n",
     "def spiral_pos_outer(t):\n",
     "    return r0 + eps * delta + eps * t / (2 * pi)\n",
     "\n",
+    "\n",
     "def spiral_neg_inner(t):\n",
     "    return r0 - eps * delta + eps / 2 + eps * t / (2 * pi)\n",
+    "\n",
+    "\n",
     "def spiral_neg_outer(t):\n",
     "    return r0 + eps * delta + eps / 2 + eps * t / (2 * pi)\n",
     "\n",
+    "\n",
     "def spiral_am1_inner(t):\n",
     "    return r0 + eps * delta + eps * t / (2 * pi)\n",
+    "\n",
+    "\n",
     "def spiral_am1_outer(t):\n",
     "    return r0 - eps * delta + eps / 2 + eps * t / (2 * pi)\n",
     "\n",
+    "\n",
     "def spiral_am2_inner(t):\n",
     "    return r0 + eps * delta + eps / 2 + eps * t / (2 * pi)\n",
+    "\n",
+    "\n",
     "def spiral_am2_outer(t):\n",
     "    return r0 - eps * delta + eps + eps * t / (2 * pi)"
    ]

docs/source/examples/notebooks/models/loss_of_active_materials.ipynb

@@ -312,6 +312,7 @@
     "def current_LAM(i, T):\n",
     "    return -1e-10 * (abs(i) + 1e3 * abs(i) ** 0.5)\n",
     "\n",
+    "\n",
     "model = pybamm.lithium_ion.DFN(\n",
     "    options=\n",
     "    {\n",