The Finite Volume Method

How PyBaMM uses the finite volume method to discretize battery PDEs into ODEs and DAEs for numerical solution.

PyBaMM uses the Finite Volume Method (FVM) to convert systems of partial differential equations (PDEs) into systems of ordinary differential equations (ODEs) and differential algebraic equations (DAEs).

Overview

The FVM discretizes the domain into a series of “volumes” and tracks the average value of each state variable $u$ across each volume.

                          +------------+------------+------------+
                          |            |            |            |
                          |   u(i-1)   |    u(i)    |   u(i+1)   |
                          |            |            |            |
                          +------------+------------+------------+
                                ^            ^            ^
                              x(i-1)       x(i)        x(i+1)
                                     ^            ^
                                  x(i-1/2)     x(i+1/2)
                                     |            |
                                  J(i-1/2)     J(i+1/2)

Each volume $v_i$ is centered at position $x_i$, with interfaces at $x_{i-\frac{1}{2}}$ and $x_{i+\frac{1}{2}}$. The state variable $u_i$ represents the average of $u$ over that volume:

$$u_i = \frac{1}{v_i} \int_\Omega u$$

Flux Calculation

Fluxes at each interface are calculated via the central difference method. For the interface at $x_{i-\frac{1}{2}}$, the flux $J_{i-\frac{1}{2}}$ is given by:

$$J_{i-\frac{1}{2}} = D_{i-\frac{1}{2}} \frac{u_i - u_{i-1}}{x_i - x_{i-1}}$$

where $D_{i-\frac{1}{2}}$ is the harmonic mean of $D_i$ and $D_{i-1}$.

We use the harmonic mean to find the values of transport properties at interfaces, and the arithmetic mean for other properties.

Rate of Change

The rate of change of $u_i$ is calculated by taking the boundary integral of the fluxes in and out of each volume:

$$\frac{du_i}{dt} = \frac{1}{v_i} \oint N(u) \cdot \mathbf{n} \, dS$$

For the 1D case, this simplifies to:

$$\frac{du_i}{dt} = \frac{J_{i-\frac{1}{2}} - J_{i+\frac{1}{2}}}{\Delta x_i}$$

Extrapolation and Boundary Values

One important consequence of using the finite volume method is that we do not directly calculate the value of state variables at the boundaries of the domain. For example, when simulating a battery, we do not directly calculate the potential of the positive electrode at the current collector (i.e., the terminal voltage of the cell). Instead, we must extrapolate from the nearest volume.

Extrapolation Methods

Constant

Uses the value from the nearest volume directly

Linear

Linearly extrapolates from the two volumes nearest the boundary (default in PyBaMM)

Quadratic

Fits a quadratic polynomial to the three volumes nearest the boundary

All of these extrapolations incur some error, which depends on:

• The steepness of the state variable near the boundary
• The size of the mesh

For maximum accuracy, we recommend using the quadratic extrapolation. Linear extrapolation is the default in PyBaMM for legacy reasons.