Modelling unequal diffusion

[nril0001]

I am trying to model a solution-based redox reaction involving unequal diffusion. The reaction is based on the equations outlined in Compton's: Understanding Voltammetry: Simulation of Electrode Processes. I have been defining the left boundary condition as Bulter Volmer and normalizing the equations in accordance with the tables below. I have tested the model and it works fine when the diffusion coefficients are equal. However, when unequal diffusion coefficients are used the model starts to act weirdly.

k0 = 1e-4
D_s = 1e-9
D_p = 1e-5
Conc S = 1e-6
Conc P = 0

Forward sweep

Under the same conditions, Digielch generates the following plot:

A massive amount of p is being generated where it should not be and I was wondering if you could provide me with some advice on my implementation.

Normalization for the Model

**

**

   # Input voltage protocol
    Edc_forward = -pybamm.t*sigma
    Edc_backwards = pybamm.t*sigma - 2 * t_reverse*sigma
    Eapp = E_start + \
        (pybamm.t <= t_reverse) * Edc_forward + \
        (pybamm.t > t_reverse) * Edc_backwards                 

    # create PyBaMM model object
    model = pybamm.BaseBatteryModel(options=seioptions)

    # Create state variables for model
    #sc is surface concentration
    c_s = pybamm.Variable("S(soln) [non-dim]", domain="solution")
    c_p = pybamm.Variable("P(soln) [non-dim]", domain="solution")
    Eeff = pybamm.Variable("Effective Voltage [non-dim]")
    
    # defining boundary values for S and P
    c_at_electrode_s = pybamm.BoundaryValue(c_s, "left")
    c_at_electrode_p = pybamm.BoundaryValue(c_p, "left")

    # Faradaic current (Butler Volmer)
    i_f = pybamm.BoundaryGradient(c_s, "left")
    butler_volmer = k0 * ((c_at_electrode_p) * pybamm.exp((1 - alpha) * (Eeff - E0))
                        - ((c_at_electrode_s) * pybamm.exp(-alpha * (Eeff - E0))))  
    
    i_cap = Cdl * Eapp.diff(pybamm.t)
    
    i = i_f + i_cap

    # PDEs - left hand side is assumed to be time derivative of the PDE
    model.rhs = {
        c_s: d_S * pybamm.div(pybamm.grad(c_s)),
        c_p: -d_P * pybamm.div(-pybamm.grad(c_p)),
    }
    
    # algebraic equations (none)
    model.algebraic = {
        Eeff: Eapp - Ru*i - Eeff,
    }
    
    # Setting boundary and initial conditions
    model.boundary_conditions = {
        c_s: {
            "right": (cs_nd, "Dirichlet"),
            "left": (-butler_volmer, "Neumann"),                    
        },

        c_p: {
            "right": (cp_nd, "Dirichlet"),
            "left": (butler_volmer, "Neumann"),                 
        } 
    }

    model.initial_conditions = {
        c_s: cs_nd,
        c_p: cp_nd,
        Eeff: E_start,
    }

    # set spatial variables and domain geometry
    x = pybamm.SpatialVariable('x', domain="solution")
    x_max = 6*pybamm.sqrt(d_S*Tmax_nd)

[rtimms]

Looks like you’re missing d_b in one of the boundary conditions, which is why it works when the diffusion coefficients are equal.

[rtimms]

Check how you’ve scaled compared to the textbook as it looks like you’ve retained both diffusion coefficients in the PDE whereas they express things in terms of the ratio of diffusion coefficients. Check how you scaled time is consistent.

[nril0001]

d_S and d_P are a ratio

Dmax = DS_d

#Diffusion coefficients
d_S = DS_d/Dmax
d_P = DP_d/Dmax

c_p: {
"right": (cp_nd, "Dirichlet"),
"left": (butler_volmer/d_P, "Neumann"),
}

Is this the change you are suggesting? because if so then when I made that change

This was the output