Skip to content

Commit cee062c

Browse files
valentinsulzervalentinsulzer
committed
#492 merge master
2 parents 84867e1 + c75de8c commit cee062c

File tree

19 files changed

+257
-73
lines changed

19 files changed

+257
-73
lines changed

‎.gitignore

+1
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
*.png
55
/local/
66
*.DS_Store
7+
*.mat
78

89
# don't ignore important .txt files
910
!requirements*

‎CHANGELOG.md

+1
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
- Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
66
- Added submodels for operating modes other than current-controlled ([#751](https://github.com/pybamm-team/PyBaMM/pull/751))
7+
- Added optional R(x) distribution in particle models ([#745](https://github.com/pybamm-team/PyBaMM/pull/745))
78
- Generalized importing of external variables ([#728](https://github.com/pybamm-team/PyBaMM/pull/728))
89
- Separated active and inactive material volume fractions ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
910
- Added submodels for tortuosity ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))

‎examples/scripts/SPMe_SOC.py

+60-51
Original file line numberDiff line numberDiff line change
@@ -27,14 +27,14 @@
2727
w = 0.207 / factor
2828
A = h * w
2929
l_s = 2.5e-5
30-
l1d = (l_n + l_p + l_s)
30+
l1d = l_n + l_p + l_s
3131
vol = h * w * l1d
3232
vol_cm3 = vol * 1e6
3333
tot_cap = 0.0
3434
tot_time = 0.0
3535
fig, axes = plt.subplots(1, 2, sharey=True)
3636
I_mag = 1.01 / factor
37-
print('*' * 30)
37+
print("*" * 30)
3838
for enum, I_app in enumerate([-1.0, 1.0]):
3939
I_app *= I_mag
4040
# load model
@@ -44,27 +44,33 @@
4444
# load parameter values and process model and geometry
4545
param = model.default_parameter_values
4646
param.update(
47-
{"Electrode height [m]": h,
48-
"Electrode width [m]": w,
49-
"Negative electrode thickness [m]": l_n,
50-
"Positive electrode thickness [m]": l_p,
51-
"Separator thickness [m]": l_s,
52-
"Lower voltage cut-off [V]": 2.8,
53-
"Upper voltage cut-off [V]": 4.7,
54-
"Maximum concentration in negative electrode [mol.m-3]": 25000,
55-
"Maximum concentration in positive electrode [mol.m-3]": 50000,
56-
"Initial concentration in negative electrode [mol.m-3]": 12500,
57-
"Initial concentration in positive electrode [mol.m-3]": 25000,
58-
"Negative electrode surface area density [m-1]": 180000.0,
59-
"Positive electrode surface area density [m-1]": 150000.0,
60-
"Current function [A]": I_app,
61-
}
47+
{
48+
"Electrode height [m]": h,
49+
"Electrode width [m]": w,
50+
"Negative electrode thickness [m]": l_n,
51+
"Positive electrode thickness [m]": l_p,
52+
"Separator thickness [m]": l_s,
53+
"Lower voltage cut-off [V]": 2.8,
54+
"Upper voltage cut-off [V]": 4.7,
55+
"Maximum concentration in negative electrode [mol.m-3]": 25000,
56+
"Maximum concentration in positive electrode [mol.m-3]": 50000,
57+
"Initial concentration in negative electrode [mol.m-3]": 12500,
58+
"Initial concentration in positive electrode [mol.m-3]": 25000,
59+
"Negative electrode surface area density [m-1]": 180000.0,
60+
"Positive electrode surface area density [m-1]": 150000.0,
61+
"Current function [A]": I_app,
62+
}
6263
)
6364
param.process_model(model)
6465
param.process_geometry(geometry)
6566
s_var = pybamm.standard_spatial_vars
66-
var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
67-
s_var.r_n: 5, s_var.r_p: 10}
67+
var_pts = {
68+
s_var.x_n: 5,
69+
s_var.x_s: 5,
70+
s_var.x_p: 5,
71+
s_var.r_n: 5,
72+
s_var.r_p: 10,
73+
}
6874
# set mesh
6975
mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
7076
# discretise model
@@ -74,58 +80,61 @@
7480
t_eval = np.linspace(0, 0.2, 100)
7581
sol = model.default_solver.solve(model, t_eval)
7682
var = "Positive electrode average extent of lithiation"
77-
xpext = pybamm.ProcessedVariable(model.variables[var],
78-
sol.t, sol.y, mesh=mesh)
83+
xpext = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
7984
var = "Negative electrode average extent of lithiation"
80-
xnext = pybamm.ProcessedVariable(model.variables[var],
81-
sol.t, sol.y, mesh=mesh)
85+
xnext = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
8286
var = "X-averaged positive particle surface concentration"
83-
xpsurf = pybamm.ProcessedVariable(model.variables[var],
84-
sol.t, sol.y, mesh=mesh)
87+
xpsurf = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
8588
var = "X-averaged negative particle surface concentration"
86-
xnsurf = pybamm.ProcessedVariable(model.variables[var],
87-
sol.t, sol.y, mesh=mesh)
88-
time = pybamm.ProcessedVariable(model.variables["Time [h]"],
89-
sol.t, sol.y, mesh=mesh)
89+
xnsurf = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
90+
time = pybamm.ProcessedVariable(
91+
model.variables["Time [h]"], sol.t, sol.y, mesh=mesh
92+
)
9093
# Coulomb counting
9194
time_hours = time(sol.t)
9295
dc_time = np.around(time_hours[-1], 3)
9396
# Capacity mAh
9497
cap = np.absolute(I_app * 1000 * dc_time)
9598
cap_time = np.absolute(I_app * 1000 * time_hours)
96-
axes[enum].plot(cap_time,
97-
xnext(sol.t), 'r-', label='Average Neg')
98-
axes[enum].plot(cap_time,
99-
xpext(sol.t), 'b-', label='Average Pos')
100-
axes[enum].plot(cap_time,
101-
xnsurf(sol.t), 'r--', label='Surface Neg')
102-
axes[enum].plot(cap_time,
103-
xpsurf(sol.t), 'b--', label='Surface Pos')
104-
axes[enum].set_xlabel('Capacity [mAh]')
99+
axes[enum].plot(cap_time, xnext(sol.t), "r-", label="Average Neg")
100+
axes[enum].plot(cap_time, xpext(sol.t), "b-", label="Average Pos")
101+
axes[enum].plot(cap_time, xnsurf(sol.t), "r--", label="Surface Neg")
102+
axes[enum].plot(cap_time, xpsurf(sol.t), "b--", label="Surface Pos")
103+
axes[enum].set_xlabel("Capacity [mAh]")
105104
handles, labels = axes[enum].get_legend_handles_labels()
106105
axes[enum].legend(handles, labels)
107106
if I_app < 0.0:
108-
axes[enum].set_ylabel('Extent of Lithiation, Elecrode Ratio: '
109-
+ str(e_ratio))
110-
axes[enum].title.set_text('Charge')
107+
axes[enum].set_ylabel(
108+
"Extent of Lithiation, Elecrode Ratio: " + str(e_ratio)
109+
)
110+
axes[enum].title.set_text("Charge")
111111
else:
112-
axes[enum].title.set_text('Discharge')
113-
print('Applied Current', I_app, 'A', 'Time',
114-
dc_time, 'hrs', 'Capacity', cap, 'mAh')
112+
axes[enum].title.set_text("Discharge")
113+
print(
114+
"Applied Current",
115+
I_app,
116+
"A",
117+
"Time",
118+
dc_time,
119+
"hrs",
120+
"Capacity",
121+
cap,
122+
"mAh",
123+
)
115124
tot_cap += cap
116125
tot_time += dc_time
117126

118-
print('Anode : Cathode thickness', e_ratio)
119-
print('Total Charge/Discharge Time', tot_time, 'hrs')
120-
print('Total Capacity', np.around(tot_cap, 3), 'mAh')
127+
print("Anode : Cathode thickness", e_ratio)
128+
print("Total Charge/Discharge Time", tot_time, "hrs")
129+
print("Total Capacity", np.around(tot_cap, 3), "mAh")
121130
specific_cap = np.around(tot_cap, 3) / vol_cm3
122-
print('Total Capacity', specific_cap, 'mAh.cm-3')
131+
print("Total Capacity", specific_cap, "mAh.cm-3")
123132
capacities.append(tot_cap)
124133
specific_capacities.append(specific_cap)
125134

126135
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
127136
ax1.plot(thicknesses / l_p, capacities)
128137
ax2.plot(thicknesses / l_p, specific_capacities)
129-
ax1.set_ylabel('Capacity [mAh]')
130-
ax2.set_ylabel('Specific Capacity [mAh.cm-3]')
131-
ax2.set_xlabel('Anode : Cathode thickness')
138+
ax1.set_ylabel("Capacity [mAh]")
139+
ax2.set_ylabel("Specific Capacity [mAh.cm-3]")
140+
ax2.set_xlabel("Anode : Cathode thickness")

‎examples/scripts/SPMe_step.py

-1
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,6 @@
4949
# plot
5050
voltage = pybamm.ProcessedVariable(
5151
model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh
52-
5352
)
5453
step_voltage = pybamm.ProcessedVariable(
5554
model.variables["Terminal voltage [V]"], step_solution.t, step_solution.y, mesh=mesh

‎examples/scripts/compare_extrapolations.py

-1
Original file line numberDiff line numberDiff line change
@@ -29,4 +29,3 @@
2929
solutions = [sim_lin.solution, sim_quad.solution]
3030
plot = pybamm.QuickPlot(models, sim_lin.mesh, solutions)
3131
plot.dynamic_plot()
32-

‎examples/scripts/compare_lithium_ion_particle_distribution.py

+82
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
#
2+
# Compare lithium-ion battery models
3+
#
4+
import argparse
5+
import numpy as np
6+
import pybamm
7+
8+
parser = argparse.ArgumentParser()
9+
parser.add_argument(
10+
"--debug", action="store_true", help="Set logging level to 'DEBUG'."
11+
)
12+
args = parser.parse_args()
13+
if args.debug:
14+
pybamm.set_logging_level("DEBUG")
15+
else:
16+
pybamm.set_logging_level("INFO")
17+
18+
# load models
19+
options = {"thermal": "isothermal"}
20+
models = [
21+
pybamm.lithium_ion.DFN(options, name="standard DFN"),
22+
pybamm.lithium_ion.DFN(options, name="particle DFN"),
23+
]
24+
25+
26+
# load parameter values and process models and geometry
27+
params = [models[0].default_parameter_values, models[1].default_parameter_values]
28+
params[0]["Typical current [A]"] = 1.0
29+
params[0].process_model(models[0])
30+
31+
32+
params[1]["Typical current [A]"] = 1.0
33+
34+
35+
def negative_distribution(x):
36+
return 1 + x
37+
38+
39+
def positive_distribution(x):
40+
return 1 + (x - (1 - models[1].param.l_p))
41+
42+
43+
params[1]["Negative particle distribution in x"] = negative_distribution
44+
params[1]["Positive particle distribution in x"] = positive_distribution
45+
params[1].process_model(models[1])
46+
47+
# set mesh
48+
var = pybamm.standard_spatial_vars
49+
var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5}
50+
51+
# discretise models
52+
for param, model in zip(params, models):
53+
# create geometry
54+
geometry = model.default_geometry
55+
param.process_geometry(geometry)
56+
mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
57+
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
58+
disc.process_model(model)
59+
60+
# solve model
61+
solutions = [None] * len(models)
62+
t_eval = np.linspace(0, 0.3, 100)
63+
for i, model in enumerate(models):
64+
solutions[i] = model.default_solver.solve(model, t_eval)
65+
66+
67+
output_variables = [
68+
"Negative particle surface concentration",
69+
"Electrolyte concentration",
70+
"Positive particle surface concentration",
71+
"Current [A]",
72+
"Negative electrode potential [V]",
73+
"Electrolyte potential [V]",
74+
"Positive electrode potential [V]",
75+
"Terminal voltage [V]",
76+
"Negative particle distribution in x",
77+
"Positive particle distribution in x",
78+
]
79+
80+
# plot
81+
plot = pybamm.QuickPlot(models, mesh, solutions, output_variables=output_variables)
82+
plot.dynamic_plot()

‎examples/scripts/heat_equation.py

+1-5
Original file line numberDiff line numberDiff line change
@@ -122,11 +122,7 @@ def T_exact(x, t):
122122
label="Numerical" if i == 0 else "",
123123
)
124124
plt.plot(
125-
xx,
126-
T_exact(xx, t),
127-
"-",
128-
color=color,
129-
label="Exact (t={})".format(plot_times[i]),
125+
xx, T_exact(xx, t), "-", color=color, label="Exact (t={})".format(plot_times[i])
130126
)
131127
plt.xlabel("x", fontsize=16)
132128
plt.ylabel("T", fontsize=16)

‎input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv

+1
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Negative electrode OCP [V],[function]graphite_mcmb2528_ocp_Dualfoil1998,
1111
Negative electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
1212
Negative electrode active material volume fraction,0.7,,assuming zero binder volume fraction
1313
Negative particle radius [m],1E-05,Scott Moura FastDFN,
14+
Negative particle distribution in x,1,,
1415
Negative electrode surface area density [m-1],180000,Scott Moura FastDFN,
1516
Negative electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
1617
Negative electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,

‎input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv

+1
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Positive electrode OCP [V],[function]lico2_ocp_Dualfoil1998,
1111
Positive electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
1212
Positive electrode active material volume fraction,0.7,,assuming zero binder volume fraction
1313
Positive particle radius [m],1E-05,Scott Moura FastDFN,
14+
Positive particle distribution in x,1,,
1415
Positive electrode surface area density [m-1],150000,Scott Moura FastDFN,
1516
Positive electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
1617
Positive electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,

‎pybamm/__init__.py

+6-1
Original file line numberDiff line numberDiff line change
@@ -100,7 +100,12 @@ def version(formatted=False):
100100
from .expression_tree.interpolant import Interpolant
101101
from .expression_tree.input_parameter import InputParameter
102102
from .expression_tree.parameter import Parameter, FunctionParameter
103-
from .expression_tree.broadcasts import Broadcast, PrimaryBroadcast, FullBroadcast
103+
from .expression_tree.broadcasts import (
104+
Broadcast,
105+
PrimaryBroadcast,
106+
FullBroadcast,
107+
ones_like,
108+
)
104109
from .expression_tree.scalar import Scalar
105110
from .expression_tree.variable import Variable
106111
from .expression_tree.independent_variable import (

‎pybamm/expression_tree/broadcasts.py

+23
Original file line numberDiff line numberDiff line change
@@ -184,3 +184,26 @@ def evaluate_for_shape(self):
184184
)
185185

186186
return child_eval * vec
187+
188+
189+
def ones_like(*symbols):
190+
"""
191+
Create a symbol with the same shape as the input symbol and with constant value '1',
192+
using `FullBroadcast`.
193+
194+
Parameters
195+
----------
196+
symbols : :class:`Symbol`
197+
Symbols whose shape to copy
198+
"""
199+
# Make a symbol that combines all the children, to get the right domain
200+
# that takes all the child symbols into account
201+
sum_symbol = symbols[0]
202+
for sym in symbols:
203+
sum_symbol += sym
204+
205+
# Just return scalar 1 if symbol has no domain (no broadcasting necessary)
206+
if sum_symbol.domain == []:
207+
return pybamm.Scalar(1)
208+
else:
209+
return FullBroadcast(1, sum_symbol.domain, sum_symbol.auxiliary_domains)

‎pybamm/models/submodels/particle/fickian/base_fickian_particle.py

-10
Original file line numberDiff line numberDiff line change
@@ -35,16 +35,6 @@ def _flux_law(self, c, T):
3535
def _unpack(self, variables):
3636
raise NotImplementedError
3737

38-
def set_rhs(self, variables):
39-
40-
c, N, _ = self._unpack(variables)
41-
42-
if self.domain == "Negative":
43-
self.rhs = {c: -(1 / self.param.C_n) * pybamm.div(N)}
44-
45-
elif self.domain == "Positive":
46-
self.rhs = {c: -(1 / self.param.C_p) * pybamm.div(N)}
47-
4838
def set_boundary_conditions(self, variables):
4939

5040
c, _, j = self._unpack(variables)

‎pybamm/models/submodels/particle/fickian/fickian_many_particles.py

+22
Original file line numberDiff line numberDiff line change
@@ -46,8 +46,30 @@ def get_coupled_variables(self, variables):
4646
N_s = self._flux_law(c_s, T_k)
4747

4848
variables.update(self._get_standard_flux_variables(N_s, N_s))
49+
50+
if self.domain == "Negative":
51+
x = pybamm.standard_spatial_vars.x_n
52+
R = pybamm.FunctionParameter("Negative particle distribution in x", x)
53+
variables.update({"Negative particle distribution in x": R})
54+
55+
elif self.domain == "Positive":
56+
x = pybamm.standard_spatial_vars.x_p
57+
R = pybamm.FunctionParameter("Positive particle distribution in x", x)
58+
variables.update({"Positive particle distribution in x": R})
4959
return variables
5060

61+
def set_rhs(self, variables):
62+
63+
c, N, _ = self._unpack(variables)
64+
65+
if self.domain == "Negative":
66+
R = variables["Negative particle distribution in x"]
67+
self.rhs = {c: -(1 / (R ** 2 * self.param.C_n)) * pybamm.div(N)}
68+
69+
elif self.domain == "Positive":
70+
R = variables["Positive particle distribution in x"]
71+
self.rhs = {c: -(1 / (R ** 2 * self.param.C_p)) * pybamm.div(N)}
72+
5173
def _unpack(self, variables):
5274
c_s = variables[self.domain + " particle concentration"]
5375
N_s = variables[self.domain + " particle flux"]

0 commit comments

Comments
 (0)