Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix SEI example notebook #3166

Merged
merged 1 commit into from
Jul 20, 2023
Merged

Fix SEI example notebook #3166

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
122 changes: 59 additions & 63 deletions docs/source/examples/notebooks/creating_models/6-a-simple-SEI-model.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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",
"$$"
]
},
{
Expand Down Expand Up @@ -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"
]
}
Expand Down Expand Up @@ -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]\")"
]
Expand All @@ -230,7 +240,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -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"
Expand All @@ -255,7 +265,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -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": [
Expand All @@ -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",
"$$"
]
},
Expand All @@ -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<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/vsulzer/Documents/Energy_storage/PyBaMM/examples/notebooks/Creating%20Models/6-a-simple-SEI-model.ipynb#X53sZmlsZQ%3D%3D?line=0'>1</a>\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"
]
Expand All @@ -348,7 +346,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -384,7 +382,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -402,7 +400,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -434,7 +432,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -467,33 +465,33 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<pybamm.models.base_model.BaseModel at 0x7f2afc6dc100>"
"<pybamm.models.base_model.BaseModel at 0x7f3a8005b490>"
]
},
"execution_count": 15,
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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": {},
Expand All @@ -565,28 +563,26 @@
"# 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",
" plt.tight_layout()\n",
" 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));"
]
},
{
Expand Down Expand Up @@ -617,7 +613,7 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 16,
"metadata": {},
"outputs": [
{
Expand All @@ -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"
]
Expand Down
Loading