Why there are difference in the results between Comsol and PyBaMM?

[ljt999]

Dear PyBaMM team,

First and foremost, I'd like to commend you for the exceptional work on PyBaMM. It runs both fast and conveniently, and I genuinely appreciate the efforts you've put into it.

I've come across a discrepancy that I hope you could give me some suggestions. Despite ensuring identical parameters (geometry, material properties, tolerance settings, etc.) in both Comsol and PyBaMM, I'm not getting consistent results. For your reference, I've configured PyBaMM to use the CasADi solver, the same one used by Comsol. I provide some comparisons of the results below:

Based on these results, could you offer some insights or recommendations regarding the discrepancies observed? I truly value your time and expertise. Thank you in advance!

[brosaplanella]

Can you share a minimum working example of the PyBaMM code?

[brosaplanella]

What is sensitive to parameters is the model (e.g. SPM, DFN...), not the parameter set (e.g. Chen 2020), so it is a matter to find if any parameters are defined slightly different in COMSOL than in PyBaMM, but in theory they should be the exact same set. I suspect the OCVs are defined the same way in COMSOL and in PyBaMM (the only difference might be whether the variable is concentration or stoichiometry), so changing the fits will not give a better agreement.

Is there any way to find which equations is COMSOL exactly solving for (as in, the exact mathematical formulation of the equations)? That would help check discrepancies with the PyBaMM model.

[ljt999]

Your assistance has been immensely helpful to me; I've resolved the issue. I greatly appreciate your help! Thank you so much!

[brosaplanella]

No worries! May I ask what the issue was, so we know it in case someone encounters it again in the future?

[ljt999]

I realized that the definitions of exchange current factors differ between COMSOL and PyBaMM. In COMSOL, I defined it using flux, so I need to multiply by the Faraday constant in PyBaMM. I appreciate your assistance again!

[Brushnie33]

Hello, I was wondering if there is a link for the comsol model that is used for comparison against the pybamm DFN model (these in this folder: https://github.com/pybamm-team/PyBaMM/tree/977dcf9cdfc2fc978da3b1b56b7423454a2a2929/examples/scripts/compare_comsol) shows perfect match. But when I used Okane2022 parameters for comparison, and I also converted reaction rate constant (multiply by 1000^0.5 and divided by Faraday's constant, 1000 is the electrolyte reference concentration). The results show good match for C-rates up to 0.50C, but there is discrepancy when C-rate goes above 1C. What could be the cause for that? @brosaplanella

[brosaplanella]

The simulation was done a few years ago by @rtimms, I believe. Not sure if he still has the file he used for validation. The fact that agreement is good at low C-rates and gets work at higher C-rates suggests some of the parameters related to the dynamics not being correct, so I would double check that. I am also a bit confused on why the curves stop at different voltages. Is it a resolution issue?

[Brushnie33]

@brosaplanella @rtimms I was suspecting that this difference might be coming from different stoichiometric window, in COMSOL, users specify minimum and maximum SOCs for each electrode, which is not explicitly specified in PyBaMM. Please neglect the stopping voltage (I didn't specify the stop voltage in COMSOL so the simulation was more time-based, but I can add that to match). I also did a comparison using half-cell set-up, the results are much better than the full-cell case above, but still see some differences.

[rtimms]

I don't have access to that file anymore, but it was for the Marquis 2019 parameters. If you're using a different parameters set and the half-cells fit well but not the full cell then it is probably the electrode balance.

You can get the stoichiometry windows by doing

parameter_values = pybamm.ParameterValues("OKane2022")
pybamm.lithium_ion.get_min_max_stoichiometries(parameter_values)

We don't set min/max sto directly, but instead calculate it to be consistent with the OCPs, electrode capacities, and voltage limits.

[Brushnie33]

@rtimms I used esoh function to calculate the stoichiometries, and the results still didn't match well. I also find out that when I increase the thickness (from 4.2e-5(Xu2019) to 7.56e-5(Okane2022), about 80% increase) of the positive electrode in the half-cell (also scaled current/capacity, but all other parameters remain the same), then the difference between the half-cell cell voltage in pybamm and comsol increases substantially. So the problem is not likely due to the stoichiometries, but maybe mesh or solver issues. Any thoughts?

[rtimms]

It’s quite hard to say without knowing more details of the implementation and parameters used. One thing to try is to simulate different limiting cases to try and debug the problem (eg low C-rate looks good so capacities and OCPs are probably fine, but things go wrong at higher rate). Try to probe whether this is kinetics, transport etc.

[Brushnie33]

@rtimms After further investigation, I have found that the issue comes from the electrolyte conductivity and diffusivity, it might be the units of c_e in comsol and pybamm are different. When I used constant values for electrolyte conductivity and diffusivity, the results match well for both half-cell and full-cell configurations. If I used an analytical equation or interpolation, there is a big mismatch for higher C rates.