Calculating current for a basic redox model

[nril0001]

Hi All,

so I am trying to model a basic redox and to do that I need to be able to model current at each step including capacitance and resistance. I am able to model the faradic current fine and include the equation effectively in the algebraic section however I run into some issues when attempting to model the capacitive current.

Capacitive current is Cdl*dVdt so the change is voltage with time, however, because Eeff is already a variable I am unable to use the Eeff.diff(pybamm.t) command to access its time derivative. I get the following error:

ModelError: time derivative of the variable found (Eeff [non-dim]') in algebraic equation Current [non-dim]

If I want to be able to model the capacitive current I need to be able to obtain an approximation of the time derivative of Eeff. If this was an array set up then I would simply use:

dEdt = (Eeff[v]-E[v-1])/deltaT eq 1

However from what I have seen this is infeasible with pybamm. I was wondering if anyone knows an alternative approach that will allow for the calculation of an approximate derivative?

I've tested the approach in eq 1 as a post-process approximation and it holds up fairly well against some other simulation software I have on hand, however, if I want to model current effectively I need to be able to calculate the current at each individual step so that I can apply it to my calculation of Eeff.

Eeff = Eapp - Ru*i eq 2

Current variable definitions and implementation have been included below. I am attempting to use a linear voltage sweep. I have also included my post-processing setup for current which is able to approximate the faradaic current fine and also introduce the effect of capacitance as well, however, I am unable to use it to introduce the impact of resistance into the system.

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


    # Create state variables for model
    sc_Ox = pybamm.Variable("O(surf) [non-dim]")
    sc_Red = pybamm.Variable("R(surf) [non-dim]")
    Eeff = pybamm.Variable("Eeff [non-dim]")
    i = pybamm.Variable("Current [non-dim]")

    # Faradaic current (Butler Volmer)     
    i_f = (k0 * (sc_Red) * np.e**((1-alpha) * (Eeff-E0))) - (k0 * (sc_Ox) * np.e**(-alpha * (Eeff-E0)))
    
    #Capacitative current
    dEdt = (((Eapp - Ru*i) - Eeff)/deltaT_nd)
    i_cap = Cdl * dEdt
    

    # PDEs - left hand side is assumed to be time derivative of the PDE
    model.rhs = {
        sc_Ox: i_f, #i_f is the echem contribution, cat_con is chemical contribution
        sc_Red: -i_f,
    }
    
    # algebraic equations (none)
    model.algebraic = {
        i: i_f - i_cap - i,
        Eeff: Eapp - Ru*i  - Eeff,
    }

    # Setting boundary and initial conditions
    model.boundary_conditions = {
    }

    model.initial_conditions = {
        sc_Ox: pybamm.Scalar(1),
        sc_Red: pybamm.Scalar(0),
        i: pybamm.Scalar(0),
        Eeff: E_start,
    }


   # Post-processing calculations of current-------------------------------------------------------------------------------------------------
        c_O = solution["O(surf) [non-dim]"](times_nd)
        E = solution["Applied Voltage [non-dim]"](times_nd)
        Cdl = solution["Cdl"](times_nd[1])
        Ru = solution["Ru"](times_nd[1])
        
        #faradaic current
        I_f = []
        I_f.append(0)
        
        #capacitative current
        I_cdl = []
        I_cdl.append(0)
        
        for v in range(1, self._m):
            
            #generate faradaic current
            dOdt = (c_O[v]-c_O[v-1])/self._deltaT_nd
            I =  dOdt 
            I_f.append(I)
            
            #generate capacitative current
            dEdt = (E[v]-E[v-1])/self._deltaT_nd
            dIdt = (I_f[v]-I_f[v-1])/self._deltaT_nd
            I_Cdl = Cdl * (dEdt - (Ru*dIdt))
            I_cdl.append(I_Cdl)
            
        current = []
        for v in range(0, len(I_f)):
            if v < len(I_f)/2:
                current.append(I_f[v] - I_Cdl)
            else:
                current.append(I_f[v] + I_Cdl)
        
        
        current = np.array(current)

[brosaplanella]

Is Eeff a variable (i.e. you need to solve for it) or a parameters (i.e. known quantity)? If the former, you should refactor the model so you get an ODE for Eeff, if the latter you should be able to get the derivative.

[nril0001]

Eeff is a variable. It is defined as Eeff = Eapp - Ru* i. Where Eapp is the applied potential defined by a voltage sweep equation, Ru is resistance and i is current. Eeff is directly impacted by current and also directly feedbacks back into the faradaic current using Butler Volmer

Eeff = Eapp - Ru*i

i = I_f + I_cp

I_f = (k0 * (sc_Red) * np.e**((1-alpha) * (Eeff-E0))) - (k0 * (sc_Ox) * np.e**(-alpha * (Eeff-E0)))

I_cp = Cdl * dEeff/dt

The only way I can think of realistically defining dEeff/dt at this current point is:

dEeff/dt = (Eeff[v]-Eeff[v-1])/deltaT eq 1

If I was to input

model.rhs = {
sc_Ox: i_f, #i_f is the echem contribution, cat_con is chemical contribution
sc_Red: -i_f,
Eeff: (Eapp - Ru*i - Eeff)/deltaT_nd
}

algebraic equations (none)

model.algebraic = {
    i: i_f - i_cap - i,
}

---

should that operate in a similar way to eq 1 in pybamm? Also if I restructure in this way, I still don't see how I could retrieve the derivative to use in the calculation of i.

Input voltage protocol

Edc_forward = -pybamm.t
Edc_backwards = pybamm.t - 2 * t_reverse
Eapp = E_start + \
    (pybamm.t <= t_reverse) * Edc_forward + \
    (pybamm.t > t_reverse) * Edc_backwards 

[brosaplanella]

But then when does the time derivative come in?

[nril0001]

the time derivative needs to be used in the capacitive current

i_cap = Cdl * dEeff/dt

i = i_f + i_cap

[valentinsulzer]

Does this work?

Eeff = pybamm.Variable("Eeff")
i = (Eapp - Eff) / Ru
i_f = ...
i_cap = i - i_f
dEffdt = i_cap / Cdl
model.rhs = {Eff: dEffdt}

[nril0001]

mathematically it makes sense but implementing it doesn't produce amazing results, unfortunately.