Why Positive electrode/Negative electrode/Electrolyte current density calculated in P2D has one more dimension than others?

[XuboGU]

Hello!

I found some variables in the solution of P2D have one more dimension than the mesh points: like Positive electrode current density [A.m-2], "Negative electrode current density [A.m-2]", and ["Electrolyte current density [A.m-2]"]. I'm wondering why these variables have one more dimension than the meshes point. Do they use different meshes when calculation?

A simple example can be as follows

import pybamm
import os
print(pybamm.__version__)

model = pybamm.lithium_ion.DFN(
    name='DFN_wo_thermal' ,
    options = {
        "particle": "Fickian diffusion", 
        # "cell geometry": "arbitrary",  # default
        # "thermal": "lumped", 
        # "particle mechanics": "swelling only",
    }
)

param = pybamm.ParameterValues('Ai2020')

capacity = 2.28 # Ah
param.update({
    "Current function [A]": capacity*pybamm.InputParameter('C-rate')
})

# set mesh (geometry)
var3 = pybamm.standard_spatial_vars
var_pts3 = {
    var3.x_n: 10,  # ne
    var3.x_s: 10,  # sep
    var3.x_p: 10,   # pe
    var3.r_n: 5,     # ne partical
    var3.r_p: 5,
}  # points

# simulation setting 
sim = pybamm.Simulation(
    model,
    var_pts = var_pts3,
    parameter_values=param,
    # solver=pybamm.CasadiSolver(mode='fast')
)

# solve for C-rates
Crates = [1]  # unit [1/h]
solutions =[]
for Crate in Crates:
    print(f'{Crate} C')
    # If t_eval provided as a list the solution is returned at 100 points within
    # the interval `[t0, tf]`.
    sol = sim.solve(t_eval=[0, 3600/Crate*1], inputs={'C-rate': Crate})
    solutions.append(sol)

solution = solutions[0]
print(solution["Positive electrode current density [A.m-2]"].entries.shape) # Is_p
print(solution["Negative electrode current density [A.m-2]"].entries.shape) # Is_n
print(solution["Electrolyte current density [A.m-2]"].entries.shape) # Ie

print(solution["Positive electrode exchange current density [A.m-2]"].entries.shape) # j_p
print(solution["Negative electrode exchange current density [A.m-2]"].entries.shape) # j_n

print(solution['Electrolyte concentration [mol.m-3]'].entries.shape) # c_e 
print(solution['Electrolyte potential [V]'].entries.shape) # phi_e

The results are like

(11, 100)
(11, 100)
(31, 100)
(10, 100)
(10, 100)
(30, 100)
(30, 100)

It can be found that the Positive electrode current density [A.m-2], "Negative electrode current density [A.m-2]", and ["Electrolyte current density [A.m-2]"] have one more dimension than the mesh points.

[valentinsulzer]

The other dimension is time

[XuboGU]

Hi Valentin. Thanks for your reply :).

I check the outputs carefully. I think the current problem is (x, t) = (30, 100), so the results like 'Electrolyte concentration [mol.m-3]'] etc provides a shape (30, 100), but the variable "Electrolyte current density [A.m-2]", somehow gives one more dimension in x, i.e., (31, 100). The additional dimension does not look like the time ..

[valentinsulzer]

Ah sorry I misunderstood. We use the finite volume method so some variables are evaluated at the center of the nodes, some at the edges. Generally those that are "flux"-like (i.e. gradient of another variable) are evaluated on edges. For example, if you have a mesh from 0 to 1 with just 2 mesh points, variables like "concentration" will be evaluated at (0.25, 0.75) and variables like "current density" will be evaluated at (0, 0.5, 1).

[valentinsulzer]

You can use solution["Electrolyte current density [A.m-2]" as an interpolator to return the values at any points in x you like and it will account for the different variable types (evaluates on centers or evaluates on edges) e.g.

solution["Electrolyte current density [A.m-2](x=np.linspace(0,L_n,100))

[XuboGU]

That makes sense. Thanks for your time, Valentin :)