fix bug in leading surface form conductivity

[rtimms]

Description

Fixes a bug where a factor of electrode surface area to volume ratio is missing in the rhs of the LeadingOrderDifferential conductivity model.

In particular, this fixes the missing semi-circle when simulating EIS

import pybamm
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft

# Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
V = model.variables["Terminal voltage [V]"]
I = model.variables["Current [A]"]
parameter_values = pybamm.ParameterValues("Marquis2019")
frequencies = np.logspace(-4, 2, 30)

# Time domain
I = 50 * 1e-3
number_of_periods = 20
samples_per_period = 16


def current_function(t):
    return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)


parameter_values["Current function [A]"] = current_function

start_time = timer.time()

sim = pybamm.Simulation(
    model, parameter_values=parameter_values, solver=pybamm.ScipySolver()
)

impedances_time = []
for frequency in frequencies:
    # Solve
    period = 1 / frequency
    dt = period / samples_per_period
    t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
    sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency})
    # Extract final two periods of the solution
    time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
    current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
    voltage = sol["Voltage [V]"].entries[-3 * samples_per_period - 1 :]
    # FFT
    current_fft = fft(current)
    voltage_fft = fft(voltage)
    # Get index of first harmonic
    idx = np.argmax(np.abs(current_fft))
    impedance = -voltage_fft[idx] / current_fft[idx]
    impedances_time.append(impedance)

end_time = timer.time()
time_elapsed = end_time - start_time
print("Time domain method: ", time_elapsed, "s")


fig, ax = plt.subplots()
data = np.array(impedances_time)
ax.plot(data.real, -data.imag, "o")
ax_max = max(data.real.max(), -data.imag.max()) * 1.1
plt.axis([0, ax_max, 0, ax_max])
plt.gca().set_aspect("equal", adjustable="box")
ax.set_xlabel(r"$Z_\mathrm{Re}$ [Ohm]")
ax.set_ylabel(r"$-Z_\mathrm{Im}$ [Ohm]")
plt.show()

Before fix:

After fix:

Fixes # (issue)

Type of change

Please add a line in the relevant section of CHANGELOG.md to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

• New feature (non-breaking change which adds functionality)
• Optimization (back-end change that speeds up the code)
• Bug fix (non-breaking change which fixes an issue)

Key checklist:

• No style issues: $ pre-commit run (or $ nox -s pre-commit) (see CONTRIBUTING.md for how to set this up to run automatically when committing locally, in just two lines of code)
• All tests pass: $ python run-tests.py --all (or $ nox -s tests)
• The documentation builds: $ python run-tests.py --doctest (or $ nox -s doctests)

You can run integration tests, unit tests, and doctests together at once, using $ python run-tests.py --quick (or $ nox -s quick).

Further checks:

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works

[brosaplanella]

Looks good, feel free to merge after adding a line to the README.

[Akila1993]

Hello Everyone,

Even though the bug is fixed, the EIS values seem erroneous for experimental data of the LGM50 Cell if you use the Chen2020 parameter set. Where the resistance at 1 kHX is 0.025 ohm in the datasheet and the EIS we are obtaining from PyBaMM is far different. I have attached the EIS plot that was obtained for 0.0001 Hz to 1 kHz and the code snipped.

This result was the same when I used the PyBaMM-EIS tool instead of time domain calculation.

Is there any suggestion or comment on this result?

import pybamm
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft

# Set up
model = pybamm.lithium_ion.SPMe(options={"surface form": "differential"}, name="SPM")
V = model.variables["Terminal voltage [V]"]
I = model.variables["Current [A]"]
parameter_values = pybamm.ParameterValues("Chen2020")
frequencies = np.logspace(-4, 2, 30)

# Time domain
I = 50 * 1e-3
number_of_periods = 20
samples_per_period = 16


def current_function(t):
    return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)


parameter_values["Current function [A]"] = current_function

start_time = timer.time()

sim = pybamm.Simulation(
    model, parameter_values=parameter_values, solver=pybamm.ScipySolver()
)

impedances_time = []
for frequency in frequencies:
    # Solve
    period = 1 / frequency
    dt = period / samples_per_period
    t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
    sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency})
    # Extract final two periods of the solution
    time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
    current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
    voltage = sol["Voltage [V]"].entries[-3 * samples_per_period - 1 :]
    # FFT
    current_fft = fft(current)
    voltage_fft = fft(voltage)
    # Get index of first harmonic
    idx = np.argmax(np.abs(current_fft))
    impedance = -voltage_fft[idx] / current_fft[idx]
    impedances_time.append(impedance)

end_time = timer.time()
time_elapsed = end_time - start_time
print("Time domain method: ", time_elapsed, "s")


fig, ax = plt.subplots()
data = np.array(impedances_time)
ax.plot(data.real, -data.imag, "o")
ax_max = max(data.real.max(), -data.imag.max()) * 1.1
plt.axis([0, ax_max, 0, ax_max])
plt.gca().set_aspect("equal", adjustable="box")
ax.set_xlabel(r"$Z_\mathrm{Re}$ [Ohm]")
ax.set_ylabel(r"$-Z_\mathrm{Im}$ [Ohm]")
plt.show()