From 7c86bebceb618da390146dd567c20fce98e7bfc7 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Wed, 19 Jul 2023 14:49:01 +0100 Subject: [PATCH] fix SEI notebook --- .../6-a-simple-SEI-model.ipynb | 122 +++++++++--------- 1 file changed, 59 insertions(+), 63 deletions(-) diff --git a/docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb b/docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb index e0125ee390..8d46012aec 100644 --- a/docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb +++ b/docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb @@ -46,17 +46,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Dimensional Model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ + "### Dimensional Model\n", + "\n", "In our simple SEI model, we consider a one-dimensional SEI which extends from the surface of a planar negative electrode at $x=0$ until $x=L$, where $L$ is the thickness of the SEI. Since the SEI is porous, there is some electrolyte within the region $x\\in[0, L]$ and therefore some concentration of solvent, $c$. Within the porous SEI, the solvent is transported via a diffusion process according to:\n", "$$\n", - "\\frac{\\partial c}{\\partial t} = - \\nabla \\cdot N, \\quad N = - D(c) \\nabla c \\label{1}\\\\\n", + "\\frac{\\partial c}{\\partial t} = - \\frac{\\partial N}{\\partial x}, \\quad N = - D(c) \\frac{\\partial c}{\\partial x} \\label{1}\\\\\n", "$$\n", "where $t$ is the time, $N$ is the solvent flux, and $D(c)$ is the effective solvent diffusivity (a function of the solvent concentration).\n", "\n", @@ -78,7 +72,22 @@ "$$\n", " R = k c|_{x=0}\n", "$$\n", - "where $k$ is the reaction rate constant (which is in general dependent upon the potential difference across the SEI)." + "where $k$ is the reaction rate constant (which is in general dependent upon the potential difference across the SEI).\n", + "\n", + "### Fixing the moving boundary\n", + "The model above is a moving boundary problem as it is defined in $x\\in[0, L]$. In order to solve it, we need to fix the boundary by introducing the scaling\n", + "$$\n", + " x = L \\xi.\n", + "$$\n", + "Then, applying the chain rule we have\n", + "$$\n", + " \\frac{\\partial}{\\partial x} \\rightarrow \\frac{1}{L} \\frac{\\partial}{\\partial \\xi}, \\quad \\text{ and } \\quad \\frac{\\partial}{\\partial t} \\rightarrow \\frac{\\partial}{\\partial t} - \\frac{L'(t)}{L(t)} \\xi \\frac{\\partial}{\\partial \\xi}.\n", + "$$\n", + "\n", + "Then, (1) can be rewritten as\n", + "$$\n", + " \\frac{\\partial c}{\\partial t} = \\frac{\\hat{V} R}{L} \\xi \\frac{\\partial c}{\\partial \\xi} + \\frac{1}{L^2} \\frac{\\partial}{\\partial \\xi} \\left( D(c) \\frac{\\partial c}{\\partial \\xi}\\right)\n", + "$$" ] }, { @@ -106,8 +115,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "\u001b[33mWARNING: You are using pip version 21.0.1; however, version 21.1.2 is available.\n", - "You should consider upgrading via the '/home/user/Documents/PyBaMM/env/bin/python3.8 -m pip install --upgrade pip' command.\u001b[0m\n", + "\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.0.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.2\u001b[0m\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } @@ -207,7 +217,7 @@ "metadata": {}, "outputs": [], "source": [ - "x = pybamm.SpatialVariable(\"x\", domain=\"SEI layer\", coord_sys=\"cartesian\")\n", + "xi = pybamm.SpatialVariable(\"xi\", domain=\"SEI layer\", coord_sys=\"cartesian\")\n", "c = pybamm.Variable(\"Solvent concentration [mol.m-3]\", domain=\"SEI layer\")\n", "L = pybamm.Variable(\"SEI thickness [m]\")" ] @@ -230,7 +240,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -239,7 +249,7 @@ "\n", "# solvent concentration equation\n", "N = - 1/L * D(c) * pybamm.grad(c)\n", - "dcdt = (V_hat * R) * pybamm.inner(x / L, pybamm.grad(c)) - 1/L * pybamm.div(N)\n", + "dcdt = (V_hat * R) / L * pybamm.inner(xi, pybamm.grad(c)) - 1/L * pybamm.div(N)\n", "\n", "# SEI thickness equation\n", "dLdt = V_hat * R" @@ -255,7 +265,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -279,18 +289,18 @@ "\n", "The boundary condition on the electrode-SEI (x=0) boundary is: \n", "$$\n", - " N|_{x=0} = - R, \\quad N|_{x=0} = - D(c|_{x=0} )\\nabla c|_{x=0}\n", + " N|_{\\xi=0} = - R, \\quad N|_{\\xi=0} = - D(c|_{\\xi=0} )\\left.\\frac{\\partial c}{\\partial \\xi}\\right|_{\\xi=0}\n", "$$\n", "which is a Neumann condition. To implement this boundary condition in pybamm, we must first rearrange the equation so that the gradient of the concentration, $\\nabla c|_{x=0}$, is the subject. Therefore we have\n", "$$\n", - " \\nabla c|_{x=0} = \\frac{R}{D(c|_{x=0} )}\n", + " \\left.\\frac{\\partial c}{\\partial \\xi}\\right|_{\\xi=0} = \\frac{R L}{D(c|_{\\xi=0} )}\n", "$$\n", "which we enter into pybamm as " ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -303,9 +313,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "On the SEI-electrolyte boundary (x=1), we have the boundary condition\n", + "On the SEI-electrolyte boundary ($\\xi=1$), we have the boundary condition\n", "$$\n", - " c|_{x=1} = c_∞\n", + " c|_{\\xi=1} = c_∞\n", "$$" ] }, @@ -319,21 +329,9 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 8, "metadata": {}, - "outputs": [ - { - "ename": "NameError", - "evalue": "name 'c_inf' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m/Users/vsulzer/Documents/Energy_storage/PyBaMM/examples/notebooks/Creating Models/6-a-simple-SEI-model.ipynb Cell 30\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m c_right \u001b[39m=\u001b[39m c_inf\n", - "\u001b[0;31mNameError\u001b[0m: name 'c_inf' is not defined" - ] - } - ], + "outputs": [], "source": [ "c_right = c_inf" ] @@ -348,7 +346,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -384,7 +382,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -402,7 +400,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -434,7 +432,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -467,16 +465,16 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 15, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -484,16 +482,16 @@ "source": [ "# define geometry\n", "geometry = pybamm.Geometry(\n", - " {\"SEI layer\": {x: {\"min\": pybamm.Scalar(0), \"max\": L_0}}}\n", + " {\"SEI layer\": {xi: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}}\n", ")\n", "\n", "def Diffusivity(cc):\n", - " return cc * 10**(-5)\n", + " return cc * 10**(-12)\n", "\n", "# parameter values (not physically based, for example only!)\n", "param = pybamm.ParameterValues(\n", " {\n", - " \"Reaction rate constant [m.s-1]\": 20,\n", + " \"Reaction rate constant [m.s-1]\": 1e-6,\n", " \"Initial thickness [m]\": 1e-6,\n", " \"Partial molar volume [m3.mol-1]\": 10,\n", " \"Bulk electrolyte solvent concentration [mol.m-3]\": 1,\n", @@ -507,7 +505,7 @@ "\n", "# mesh and discretise\n", "submesh_types = {\"SEI layer\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {x: 100}\n", + "var_pts = {xi: 100}\n", "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", " \n", "spatial_methods = {\"SEI layer\": pybamm.FiniteVolume()}\n", @@ -517,13 +515,13 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# solve\n", "solver = pybamm.ScipySolver()\n", - "t = np.linspace(0, 100, 100) # solve for 100s\n", + "t = [0, 100] # solve for 100s\n", "solution = solver.solve(model, t)\n", "\n", "# post-process output variables\n", @@ -541,18 +539,18 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "9a36511202a146fdbc57f80547747cbe", + "model_id": "0f4941d4712049e494267074dca70b4b", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t', max=10.0), Output()), _dom_classes=('widget-inte…" + "interactive(children=(FloatSlider(value=0.0, description='t'), Output()), _dom_classes=('widget-interact',))" ] }, "metadata": {}, @@ -565,20 +563,18 @@ "# plot SEI thickness in microns as a function of t in microseconds\n", "# and concentration in mol/m3 as a function of x in microns\n", "L_0_eval = param.evaluate(L_0)\n", - "x = np.linspace(0, L_0_eval, 100) # dimensionless space\n", + "xi = np.linspace(0, 1, 100) # dimensionless space\n", "\n", "def plot(t):\n", - " t_in_seconds = t / 1e6\n", - " f, (ax1, ax2) = plt.subplots(1, 2 ,figsize=(10,5))\n", - " ax1.plot(solution.t * 1e6, L_out(solution.t) * 1e6)\n", - " ax1.plot(t, L_out(t_in_seconds) * 1e6, 'r.')\n", - " ax1.set_xlim(0, 10)\n", + " _, (ax1, ax2) = plt.subplots(1, 2 ,figsize=(10,5))\n", + " ax1.plot(solution.t, L_out(solution.t) * 1e6)\n", + " ax1.plot(t, L_out(t) * 1e6, 'r.')\n", " ax1.set_ylabel(r'SEI thickness [$\\mu$m]')\n", - " ax1.set_xlabel(r't [$\\mu$s]') \n", + " ax1.set_xlabel(r't [s]') \n", " \n", - " plot_c, = ax2.plot(x * 1e6 * L_out(t_in_seconds), c_out(t_in_seconds, x_in_metres))\n", + " ax2.plot(xi * L_out(t) * 1e6, c_out(t, xi))\n", " ax2.set_ylim(0, 1.1)\n", - " ax2.set_xlim(0, x[-1] * 1e6 * L_out(solution.t[-1])) \n", + " ax2.set_xlim(0, L_out(solution.t[-1]) * 1e6) \n", " ax2.set_ylabel('Solvent concentration [mol.m-3]')\n", " ax2.set_xlabel(r'x [$\\mu$m]')\n", "\n", @@ -586,7 +582,7 @@ " plt.show()\n", " \n", "import ipywidgets as widgets\n", - "widgets.interact(plot, t=widgets.FloatSlider(min=0,max=solution.t[-1]*1e6,step=0.1,value=0));" + "widgets.interact(plot, t=widgets.FloatSlider(min=0,max=solution.t[-1],step=0.1,value=0));" ] }, { @@ -617,7 +613,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -626,7 +622,7 @@ "text": [ "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\n", + "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", "\n" ]