<jupyter_start><jupyter_text>A step towards the Single Particle Model In the [previous notebook](./2-a-pde-model.ipynb) we saw how to solve a PDE model in pybamm. Now it is time to solve a real-life battery problem! We consider the problem of spherical diffusion in the negative electrode particle within the single particle model. That is,\begin{equation*}  \frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c),\end{equation*}with the following boundary and initial conditions:\begin{equation*}  \left.\frac{\partial c}{\partial r}\right\vert_{r=0} = 0, \quad \left.\frac{\partial c}{\partial r}\right\vert_{r=R} = -\frac{j}{FD}, \quad \left.c\right\vert_{t=0} = c_0,\end{equation*}where $c$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. In this example we use the following parameters:| Symbol | Units              | Value                                          ||:-------|:-------------------|:-----------------------------------------------|| $R$      | m                | $10 \times 10^{-6}$                            || $D$      | m${^2}$ s$^{-1}$ | $3.9 \times 10^{-14}$                          || $j$      | A m$^{-2}$       | $1.4$                                          || $F$      | C mol$^{-1}$     | $96485$                                        || $c_0$    | mol m$^{-3}$     | $2.5 \times 10^{4}$                            | Note that all battery models in PyBaMM are written in dimensionless form for better numerical conditioning This is discussed further in [the simple SEI model notebook](./5-a-simple-SEI-model.ipynb).  Setting up the modelAs before, we begin by importing the pybamm library into this notebook, along with any other packages we require, and start with an empty `pybamm.BaseModel`<jupyter_code>%pip install pybamm -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import matplotlib.pyplot as plt

model = pybamm.BaseModel()<jupyter_output><empty_output><jupyter_text>We then define all of the model variables and parameters. Parameters are created using the `pybamm.Parameter` class and are given informative names (with units). Later, we will provide parameter values and the `Parameter` objects will be turned into numerical values. For more information please see the [parameter values notebook](../parameter-values.ipynb).<jupyter_code>R = pybamm.Parameter("Particle radius [m]")
D = pybamm.Parameter("Diffusion coefficient [m2.s-1]")
j = pybamm.Parameter("Interfacial current density [A.m-2]")
F = pybamm.Parameter("Faraday constant [C.mol-1]")
c0 = pybamm.Parameter("Initial concentration [mol.m-3]")

c = pybamm.Variable("Concentration [mol.m-3]", domain="negative particle")<jupyter_output><empty_output><jupyter_text>Now we define our model equations, boundary and initial conditions, as in the previous example.<jupyter_code># governing equations
N = -D * pybamm.grad(c)  # flux
dcdt = -pybamm.div(N)
model.rhs = {c: dcdt}  

# boundary conditions 
lbc = pybamm.Scalar(0)
rbc = -j / F / D
model.boundary_conditions = {c: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")}}

# initial conditions 
model.initial_conditions = {c: c0}<jupyter_output><empty_output><jupyter_text>Finally, we add any variables of interest to the dictionary `model.variables`<jupyter_code>model.variables = {
    "Concentration [mol.m-3]": c,
    "Surface concentration [mol.m-3]": pybamm.surf(c),
    "Flux [mol.m-2.s-1]": N,
}<jupyter_output><empty_output><jupyter_text>Using the model In order to discretise and solve the model we need to provide values for all of the parameters. This is done via the `pybamm.ParameterValues` class, which accepts a dictionary of parameter names and values<jupyter_code>param = pybamm.ParameterValues(
    {
        "Particle radius [m]": 10e-6,
        "Diffusion coefficient [m2.s-1]": 3.9e-14,
        "Interfacial current density [A.m-2]": 1.4,
        "Faraday constant [C.mol-1]": 96485,
        "Initial concentration [mol.m-3]": 2.5e4,
    }
)<jupyter_output><empty_output><jupyter_text>Here all of the parameters are simply scalars, but they can also be functions or read in from data (see  [parameter values notebook](../parameter-values.ipynb)). As in the previous example, we define the particle geometry. Note that in this example the definition of the geometry contains a parameter, the particle radius $R$<jupyter_code>r = pybamm.SpatialVariable("r", domain=["negative particle"], coord_sys="spherical polar")
geometry = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": R}}}<jupyter_output><empty_output><jupyter_text>Both the model and geometry can now be processed by the parameter class. This replaces the parameters with the values<jupyter_code>param.process_model(model)
param.process_geometry(geometry)<jupyter_output><empty_output><jupyter_text>We can now set up our mesh, choose a spatial method, and discretise our model<jupyter_code>submesh_types = {"negative particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)}
var_pts = {r: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {"negative particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model);<jupyter_output><empty_output><jupyter_text>The model is now discretised and ready to be solved.  Solving the model As is the previous example, we choose a solver and times at which we want the solution returned.<jupyter_code># solve
solver = pybamm.ScipySolver()
t = np.linspace(0, 3600, 600)
solution = solver.solve(model, t)

# post-process, so that the solution can be called at any time t or space r
# (using interpolation)
c = solution["Concentration [mol.m-3]"]
c_surf = solution["Surface concentration [mol.m-3]"]

# plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))

ax1.plot(solution.t, c_surf(solution.t))
ax1.set_xlabel("Time [s]")
ax1.set_ylabel("Surface concentration [mol.m-3]")

r = mesh["negative particle"].nodes # radial position
time = 1000  # time in seconds
ax2.plot(r * 1e6, c(t=time, r=r), label="t={}[s]".format(time))
ax2.set_xlabel("Particle radius [microns]")
ax2.set_ylabel("Concentration [mol.m-3]")
ax2.legend()

plt.tight_layout()
plt.show()<jupyter_output>2020-05-30 11:11:30,262 - [WARNING] processed_variable.get_spatial_scale(497): No scale set for spatial variable r_n. Using default of 1 [m].
{'fraction_non_alphanumeric': 0.13573193329215566, 'fraction_numerical': 0.019765287214329835, 'mean_word_length': 3.8919939577039275, 'pattern_counts': {'":': 15, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '189525', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2120, 'n_tokens_neox': 2008, 'n_words': 649, 'path': 'examples/notebooks/Creating Models/3-negative-particle-problem.ipynb', 'repo_name': 'DrSOKane/PyBaMM', 'size': '6476'}