Modeling porous electrodes: Part 3 (Newman’s BAND method)

I previously made a post about the basic formulation of porous electrode theory (part 1), and a post about how the choice of kinetic expression determines what the result will look like (part 2). The purpose of this post is to explain how to solve these basic problems using the BAND formulation found in John Newman’s textbook Electrochemical Systems (Newman and Thomas-Alyea in the Third Edition). I also provide Python code so you can get the solution yourself and modify it for your own teaching and learning needs.

As a model problem we’ll use the boundary value problem of Tafel kinetics in a concentration-independent porous electrode in Cartesian coordinates. The reason to choose this one is because Newman and Tobias provide the analytical solution in their 1962 paper “Theoretical Analysis of Current Distribution in Porous Electrodes.” You can also find this solution in Fuller and Harb’s Electrochemical Engineering, in the problems at the end of Chapter 5. The reason to choose a problem with an analytical solution is so you can check your work to know if it’s right! We lay out the problem below, which is the same as that in Newman and Tobias (1962) except one sign convention. This is four equations with four unknowns (i1, i2, Φ1, Φ2):

If you read the electrochemical engineering literature, solving problems such as this with Newman’s BAND method often comes up, and you might wonder what that is. Technically “BAND” is some Fortran code set up to solve 2-point boundary value problems in a way that is convenient for electrochemists. (My own favorite document written about BAND is “Modeling and reactor simulation” by Douglas Bennion in the AIChE symposium series 1983, Vol 79, Num 229, pp 25-36. It’s a resource meant for teachers and is somewhat hard to find so I’ve uploaded it here.) The formulation Newman uses to set up the problem casts each of the equations in this way:

Here c is a variable. The coefficients a, b, and d are those for the second derivative, first derivative, and the variable itself. And g is the constant term. When you linearize equation 2 above and then cast all of the four equations this way, you get something that looks like the tables below. I’ve collected some constants together as ß and P to keep it from looking too messy:

Once you have this, you can (a) build a block tridiagonal matrix, (b) choose initial guesses for the answers, and (c) iterate to a solution by repeatedly solving the problem until your initial guesses and final answers match. Explaining how this works is too long for a blog post, so follow these links for:

If you execute the code you should get a figure like that below, which solves the problem. (I still use Python 2.7 so you might have to do some mild editing and debugging.)

What I use this for

The reason I started writing a BAND-type method in Python is (1) I wanted to learn more about how BAND worked so I could teach it better, and (2) I wanted to model electrochemical systems without having to buy expensive software or use Fortran. It’s not exactly like Newman’s code, which also included a function to solve the matrix problem. Rather, this code lets you input the problem in Newman’s way and build the corresponding block tridiagonal matrix, then it uses SciPy to solve the matrix problem.

When I was working on the paper above in 2016, I wanted to compare experimental data to the Chen + Podlaha + Cheh alkaline battery model developed in the 1990s. Python seems like the best choice for something you can also develop as a teaching tool, so I used that to generate the results in Figure 3 of the paper. Since then I’ve used it to train some MS students getting insight into shallow-cycled MnO2 cathodes. When Zhicheng Lu’s MS thesis appeared on the internet, I received a few questions about the model by email. So this post is meant to help anyone doing similar things.