Super high c(Li+) when using exponential form electrolyte property and low porosity

[RuiheLi]

Dear PyBaMM,

I tried to run 1C discharge with Chen2020 and an exponential form of diffusivity and conductivity (based on Nyman2008) and set the negative porosity to be 0.1.
Super high c(Li+) occurs:

I tried:

1. increase mesh point
2. change tolerance "D_e__c_e" from 10 to 1e-15

But both didn't work. Would you think that this is a numerical problem? I wonder whether any of you have any ideas on dealing with this problem. Any suggestions and ideas will be really appreciated! Big thanks in advance!

Yours Sincerely,
Ruihe Li

Script is pasted here FYI.

# import modules
import pybamm
import numpy as np;import matplotlib.pyplot as plt;
from pybamm import constants,exp,tanh,sqrt;

from pybamm import exp,sqrt
def electrolyte_diffusivity_Nyman2008Exp(c_e, T):
    D_c_e = (
        6 * exp( -1 *(c_e/1000)) 
        * exp(-5/298.15) 
        * exp(-95/298.15*(c_e/1000)) * 1e-10 
    )
    return D_c_e

def electrolyte_conductivity_Nyman2008Exp(c_e, T):
    sigma_e = (
        0.1 * 0.06248 * (1+298.15-0.05559) * 
        (c_e/1e3) * (1 - 3.084 *sqrt(c_e/1e3) 
        + 1.33 *(1+ 0.03633 *(exp(1000/298.15))*c_e/1e3)   ) 
        / (1+(c_e/1e3)**4*( 0.00795 *exp(1000/298.15))) 
    )
    return sigma_e
print(electrolyte_diffusivity_Nyman2008Exp(10,298.15))
print(electrolyte_conductivity_Nyman2008Exp(10,298.15))

# set-up and run model
V_max = 4.2;        V_min = 2.5;
total_cycle = 1; save_at_cycles= 1
Experiment_Long   = pybamm.Experiment( [ (
    f"Discharge at 1 C until {V_min} V (1 minutes period)" , )] * total_cycle ) 
Experiment_s   = pybamm.Experiment( [ (
    f"Discharge at 1 C for 20 minutes (10 seconds period)" , )] * total_cycle ) 

Sim = [];Sol=[]; 

Para_0 = pybamm.ParameterValues("Chen2020")

Para_0['Negative electrode porosity'] =  0.1
var = pybamm.standard_spatial_vars  
mesh_thick = 30; mesh_par=30;
var_pts = {
    var.x_n: int(mesh_thick),  
    var.x_s: int(20),  
    var.x_p: int(mesh_thick),  
    var.r_n: int(mesh_par),  
    var.r_p: int(mesh_par),  }  

pybamm.settings.tolerances['D_e__c_e'] = 10

# Run the first model, with Nyman2008 polynimial form
model = pybamm.lithium_ion.DFN()
c_e = model.variables["Electrolyte concentration [mol.m-3]"]
T = model.variables["Cell temperature [K]"]
model.variables["Electrolyte conductivity [S.m-1]"] =(
    Para_0['Electrolyte conductivity [S.m-1]'](c_e, T))
model.variables["Electrolyte diffusivity [m2.s-1]"] =(
    Para_0['Electrolyte diffusivity [m2.s-1]'](c_e, T))

sim = pybamm.Simulation(
    model, experiment = Experiment_Long,
    parameter_values = Para_0,
    solver = pybamm.CasadiSolver(),
    var_pts=var_pts,)  
try:
    sol = sim.solve(save_at_cycles=save_at_cycles,);
    print(sol.cycles[-1].steps[-1]);  # a way to check whether the solution is finalized 
except:
    print('Fail for electrolyte: Poly')
else:
    Sol.append(sol)   
    print('Succeed for electrolyte: Poly')

# Run the first model, with Nyman2008 exponeitial form
model_2 = pybamm.lithium_ion.DFN()
Para_0.update({
    'Electrolyte conductivity [S.m-1]': 
    electrolyte_conductivity_Nyman2008Exp})
Para_0.update({
    'Electrolyte diffusivity [m2.s-1]': 
    electrolyte_diffusivity_Nyman2008Exp})
c_e = model_2.variables["Electrolyte concentration [mol.m-3]"]
T = model_2.variables["Cell temperature [K]"]
model_2.variables["Electrolyte conductivity [S.m-1]"] =(
    Para_0['Electrolyte conductivity [S.m-1]'](c_e, T))
model_2.variables["Electrolyte diffusivity [m2.s-1]"] =(
    Para_0['Electrolyte diffusivity [m2.s-1]'](c_e, T))

sim_2 = pybamm.Simulation(
    model_2, experiment = Experiment_Long,
    parameter_values = Para_0,
    solver = pybamm.CasadiSolver(),
    var_pts=var_pts,)  
try:
    sol_2 = sim_2.solve(save_at_cycles=save_at_cycles,);
    print(sol_2.cycles[-1].steps[-1]);  # a way to check whether the solution is finalized 
except:
    print('Fail for electrolyte: Exp')
else:
    Sol.append(sol_2)   
    print('Succeed for electrolyte: Exp')

label = ['Nyman - Poly', 'Nyman - Exp', ]
output_variables3 = [
    "Terminal voltage [V]",   
    "Electrolyte potential [V]",
    "Electrolyte concentration [mol.m-3]",
    "Electrolyte flux [mol.m-2.s-1]",
    "Electrolyte diffusivity [m2.s-1]",
    "Electrolyte conductivity [S.m-1]"
]
quick_plot = pybamm.QuickPlot(
    [sol for sol in Sol], 
    output_variables3,label,
    variable_limits='tight',time_unit='minutes',n_rows=2,
    figsize = (12,8)) #     spatial_unit='mm',
quick_plot.dynamic_plot();