Capturing a the dynamics of a battery,  with mathematics.

In one of my two undergrad senior year research projects, I present a mean-field model of the nickel-hydroxide (Ni(OH)2) electrode in the Edison rechargeable battery. The model aims to describe the electrode’s dynamics, charge transport, and phase separation. The Poisson-Nernst-Planck equation is used to describe the charge transport, while the Cahn-Hilliard approach is used to model the phase separation. Then, the model equations are solved numerically. Before presenting the project, let’s first discuss an important permilmeniry work I conducted.

Numerically Solving the PNP Equations

While my advisor and I were working on the mathematical model, I also tried to tackle my first numerical problem - solving the Poission Nernst Planck equations. I provided quite a detailed solution, which I will present here.

We would like to solve the following equations (in 1D):

∂p∂t=∇2p+∇(p∇ϕ)∂n∂t=∇2n−∇(n∇ϕ)∇2ϕ=n−p\begin{align} \frac{\partial p}{\partial t} &= \nabla^2 p + \nabla(p\nabla \phi) \\ \frac{\partial n}{\partial t} &= \nabla^2 n - \nabla(n\nabla \phi) \\ \nabla^2 \phi &= n - p \end{align}

Subject to the boundary conditions:

∇p+p∇ϕ=∇n−n∇ϕ=0,x=L,0ϕ(Lx)=1,ϕ(0)=−1\begin{align} \nabla p + p \nabla \phi &= \nabla n - n \nabla \phi = 0, x=L,0 \\ \phi(L_x) &= 1, \phi(0) = -1 \end{align}

To find ϕ\phi for all xx at the current step, I will use the implicit finite difference method. After applying finite difference, the potential ϕ(x,t)\phi(x,t) at the point xm,0<m<M−1x_m, 0<m<M-1 would be:

Δx22(pm−nm)=ϕm−12ϕm+1−12ϕm−1\begin{equation} \frac{\Delta x^2}{2}(p_m - n_m) = \phi_m -\frac{1}{2} \phi_{m+1} -\frac{1}{2} \phi_{m-1} \end{equation}

Therefore, solving for ϕ\phi at a certain time would mean solving the following matrix:

[1−1200...−121−120............0−121−12...00−121][ϕ1ϕ2...ϕM−3ϕM−2]=Δx22[p1−n1+1Δx2ϕ0p2−n2...pM−3−nM−3pM−2−nM−2+1Δx2ϕM]\begin{equation} \begin{aligned} \begin{bmatrix} 1 & -\frac{1}{2} & 0 & 0 & ... \\ -\frac{1}{2} & 1 & -\frac{1}{2} & 0 & ... \\ ... & ... \\ ... & 0 & -\frac{1}{2} & 1 & -\frac{1}{2} \\ ... & 0 & 0 & -\frac{1}{2} & 1 \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ ... \\ \phi_{M-3} \\ \phi_{M-2} \end{bmatrix} = \\ \frac{\Delta x^2}{2} \begin{bmatrix} p_1 - n_1 + \frac{1}{\Delta x^2} \phi_0 \\ p_2 - n_2 \\ ... \\ p_{M-3} - n_{M-3} \\ p_{M-2} - n_{M-2} + \frac{1}{\Delta x^2} \phi_M \end{bmatrix} \end{aligned} \end{equation}

Notice that the LHS matrix is composed of (M−2)(M-2) ones, and (M−3)(M-3) elements −12-\frac{1}{2} on each diagonal.

In python, it looks like this:

def get_phi(phi_0, phi_M, n, p, dx, n_l=0, n_r=0):
    # TODO: add phi_neumann conditions

    # Creating the coeff matrix
    coeff_matrix = np.diag(np.ones(M-2)) + np.diag(-0.5*np.ones(M-3), 1) + np.diag(-0.5*np.ones(M-3), -1)

    # Creating the RHS result vector
    result_vector = np.array([p[m]-n[m] for m in range(2, M-2)])
    right_res = p[M-2] - n[M-2] + phi_M/(dx**2)
    left_res = p[1] - n[1] + phi_0/(dx**2)
    result_vector = ((dx**2)/2)*np.append(np.insert(result_vector, 0, left_res), right_res)

    # Solving for phi and adding the boundry
    phi_trimmed = np.linalg.solve(coeff_matrix, result_vector)
    phi = np.append(np.insert(phi_trimmed, 0, phi_0), phi_M)

    return phi

Together with finding ϕ\phi, we need to iterate p(x,t),n(x,t)p(x,t), n(x,t). This involves finding pmi+1,nmi+1p_m^{i+1}, n_m^{i+1} for each m,j<im, j<i. We assume that we have pmj,nmjp_m^j, n_m^j for every m,j<im, j<i. I will use explicit central difference method to get:

pmi+1=pmi+ΔtΔx2×[pm+1i+pm−1i−2pmi+14(pm+1i−pm−1i)(ϕm+1i−ϕm−1i)+pmi(ϕm+1i+ϕm−1i−2ϕmi)]\begin{equation*} \begin{aligned} p_m^{i+1} =& p_m^{i} + \frac{\Delta t}{\Delta x^2} \times \\ &[ p_{m+1}^i + p_{m-1}^i - 2p_m^i \\ &+ \frac{1}{4}(p_{m+1}^i-p_{m-1}^i)(\phi_{m+1}^i-\phi_{m-1}^i) \\ &+ p_m^i(\phi_{m+1}^i+\phi_{m-1}^i-2\phi_m^i)] \end{aligned} \end{equation*}
nmi+1=nmi+ΔtΔx2×[nm+1i+nm−1i−2nmi−14(nm+1i−nm−1i)(ϕm+1i−ϕm−1i)−nmi(ϕm+1i+ϕm−1i−2ϕmi)]\begin{equation*} \begin{aligned} n_m^{i+1} =& n_m^{i} + \frac{ \Delta t}{\Delta x^2} \times \\ & [ n_{m+1}^i + n_{m-1}^i - 2n_m^i \\ &- \frac{1}{4}(n_{m+1}^i-n_{m-1}^i)(\phi_{m+1}^i-\phi_{m-1}^i) \\ &- n_m^i(\phi_{m+1}^i+\phi_{m-1}^i-2\phi_m^i) ] \end{aligned} \end{equation*}

This can be implemented nicely using NumPy:

def iterate_concentration(c, phi, dx, dt, positive, flux_l=0, flux_r=0, neumann=False):
    """
    Assume c and phi are 1XM vectors 
    c_pp - c-prime-prime, cp - c-prime, phi_pp - phi-prime-prime, phi_p - phi-prime
    
    """
    sign = 1 if positive else -1

    c_1, c_M = c[0], c[M-1]
    phi_1, phi_M = phi[0], phi[M-1]
    # c_prev_1, c_prev_M = c_prev[0], c_prev[M-1]


    c_rolled = np.roll(c, 1)
    c_rolled_back = np.roll(c, -1)
    np.put(c_rolled, [0, -1], [c_1, c_M]) # Creating c_{m+1}
    np.put(c_rolled_back, [0, -1], [c_1, c_M]) # Creating c_{m-1}

    phi_rolled = np.roll(phi, 1)
    phi_rolled_back = np.roll(phi, -1)
    np.put(phi_rolled, [0, -1], [phi_1, phi_M]) # Creating phi_{m+1}
    np.put(phi_rolled_back, [0, -1], [phi_1, phi_M]) # Creating phi_{m+1}

    c_pp = c_rolled + c_rolled_back - 2*c
    c_p = c_rolled - c_rolled_back
    
    phi_pp = phi_rolled + phi_rolled_back - 2*phi
    phi_p = phi_rolled - phi_rolled_back
    # print(phi_p, phi_pp)

    c_iter = c + dt/(dx**2)*(c_pp + (1/4)*sign*c_p*phi_p + sign*c*phi_pp)
    
    if neumann:
        c_iter[0] = (dx*flux_l - c_iter[1])/(sign*(phi[1]-phi[0])-1)
        c_iter[-1] = c_iter[-2] - sign*c_iter[-2]*(phi[-1]-phi[-2]) + dx*flux_r


    return c_iter

This eventually, results in this beautiful animation

Nickel-Hydroxide Project

We started from a different set of mathematical equations, but quickly realized that we have to modify the mathematical model to our specific setup. Eventually, the model equations took the following form:

0=−n+∂x(ϵ∂xϕ)\begin{equation*} 0 = -n + \partial_x (\epsilon \partial_x \phi) \end{equation*}
∂n∂t=Dn∂x⋅(∂xn−n∂xϕ)−α2×((1−u)n−(1+u)(1−n))\begin{equation*} \begin{aligned} \frac{\partial n}{\partial t} =& D_n \partial_x \cdot (\partial_x n - n \partial_x \phi) - \frac{\alpha}{2} \times \\ &\big( (1-u)n - (1+u)(1-n) \big) \end{aligned} \end{equation*}
∂u∂t=Du∂x⋅(1−u2)×[11−u2∂xu+βW′′(u)∂xu+λ∂xψ]+α2((1−u)n−(1+u)(1−n))\begin{equation*} \begin{aligned} \frac{\partial u}{\partial t} =& D_u \partial_x \cdot (1-u^2) \times \\ & \left[ \frac{1}{1-u^2} \partial_x u + \beta W''(u) \partial_x u \\+ \lambda \partial_x \psi \right] \\ & + \frac{\alpha}{2} \big( (1-u)n - (1+u)(1-n) \big) \end{aligned} \end{equation*}
ψ=∂x2u\begin{equation*} \psi = \partial_x ^2 u \end{equation*}

Essentialy, the equations describe the motion of the ions inside the solution, affected by diffusion, electric forces, and concentration gradients. A full description of the model, its parameters, and the code, can be found in my GitHub.

An essential part of the dynamics is the phase-transition the material undergoes (modeled with Cahn-Hilliard equations). We can observe this dynamic process in action in the following figure:

The figure was created in MatLab and it shows a Space-time plot with an input electron flux (to the electrode) from the right. The different phases are Purple - Ni(OH)2{Ni(OH)_2}, Blue - NiOOH{NiOOH}. Over time, the Ni(OH)2{Ni(OH)_2} phase takes over the material and the interface moves to the left edge.

After submitting the project and continuing with my studies, my advisor informed me that he and his colleagues further developed my work into a 2D model, achieving the exact behavior they sought. I hope to continue supporting their efforts and contribute to their research, aiming to publish preliminary results in the upcoming couple of years.