A Hybrid Backward Euler Control Volume Method to Solve the Concentration-Dependent Solid-State Diffusion Problem in Battery Modeling

Kudakwashe Chayambuka^1,2,3, Grietus Mulder^1,2, Dmitri L. Danilov^3,4, Peter H. L. Notten^3,4,5*

¹VITO, Mol, Belgium.
²EnergyVille, Genk, Belgium.
³Eindhoven University of Technology, Eindhoven, Netherlands.
⁴Forschungszentrum Jülich, Fundamental Electrochemistry (IEK-9), Jülich, Germany.
⁵University of Technology Sydney, Broadway, Sydney, NSW, Australia.

DOI: 10.4236/jamp.2020.86083                

Abstract

Several efficient analytical methods have been developed to solve the solid-state diffusion problem, for constant diffusion coefficient problems. However, these methods cannot be applied for concentration-dependent diffusion coefficient problems and numerical methods are used instead. Herein, grid-based numerical methods derived from the control volume discretization are presented to resolve the characteristic nonlinear system of partial differential equations. A novel hybrid backward Euler control volume (HBECV) method is presented which requires only one iteration to reach an implicit solution. The HBECV results are shown to be stable and accurate for a moderate number of grid points. The computational speed and accuracy of the HBECV, justify its use in battery simulations, in which the solid-state diffusion coefficient is a strong function of the concentration.

Keywords

Solid-State Diffusion, Implicit Methods, Backward Euler

1. Introduction

Since the discovery of the electrochemical intercalation reaction of lithium in layered titanium disulfide (TiS₂) by Whittingham in 1976, lithium-ion batteries (LIB) have grown in popularity and application to surpass all other electrical energy storage devices [1] [2]. This groundbreaking discovery demonstrated the reversible shuttling of lithium ions in solid layered materials and sparked the development of intercalation electrode materials for room temperature secondary batteries. However, due to the low voltage of TiS₂ as a cathode material (approximately 2 V vs. Li⁺/Li) and the lack of a safe anode the time, the Li/TiS₂ battery was not commercialized.

The intercalation mechanism would later be revisited in the search of a safe anode, which completely avoids the danger of lithium dendrites. These needle-like structures protrude from the metallic anode surface during platting, causing internal short circuits and battery-related fires. The coupling of two intercalation electrode materials is therefore, the fundamental basis of the ubiquitous LIB which are used in most consumer electronics and have recently captured public attention due to their application in electric vehicles. It is therefore meritorious, for the 2019 Nobel Prize in chemistry, to be awarded to John Goodenough, M. Stanley Whittingham and Akira Yoshino, for developing the LIB [3]. We should also mention the nickel metal hydride batteries (NiMH), the emerging sodium-ion batteries (SIB) and potassium-ion batteries, which likewise make use of intercalation active materials.

Inasmuch as intercalation of energy carrying species inside electrodes is fundamental to the operation of the aforementioned battery chemistries, the mathematical modeling of this mechanism is equally critical to the success of physics-based battery models. The time-dependent radial transport of species inside spherical active particles is governed by Fick’s second law

$\frac{\partial c}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(D{r}^2\frac{\partial c}{\partial r}\right),$ (1)

where $c=c\left(r,t\right)$ is the concentration of the intercalated species [mol·m⁻³], $D=D\left(c\right)$ is the solid-state diffusion coefficient [m²·s⁻¹], r is the radial distance from the center of the particle [m] and t is time [s]. Equation (1) does not assume a constant diffusion coefficient, which indeed, may vary as a function of the concentration. This has consequences on the complexity of the numerical solution method as we shall later explore in detail. For a physics-based battery model, we should highlight that the ultimate goal is to derive the particle surface concentration for a given surface flux boundary condition (Neumann boundary condition) using a numerical solution method that is both robust and computationally inexpensive. The surface concentration is the most important parameter which governs reaction kinetics and the state-of-charge in electrochemical battery models.

Several efficient methods for solving Equation (1) in the case of a constant diffusion coefficient exist and have been applied in many battery models [4] [5] [6]. A constant diffusion coefficient is the customary assumption in battery models, a simplification which allows analytical and numerical solutions from heat transfer theory to be used [7]. Such numerical methods are detailed in the classic reviews by Subramanian et al. [8] and Zeng et al. [9]. In our assessment, Liu’s pseudo steady-state method (PSS) resolves the seemingly antagonistic requirements of computational speed and numerical accuracy [4]. A recent modification to the PSS method, provides faster and stable results for long time simulations [5].

A challenging problem which has remained relatively unscrutinized for many years, is how to efficiently simulate a variable diffusion coefficient [10] [11] [12]. This problem deserves increased attention because of the need to address stress effects in particles, [13] phase separation in two-phase materials and temperature effects occurring at high discharge rates. Overwhelming evidence corroborates that a constant diffusion coefficient is a rare exception in intercalation particles. For example, the diffusion coefficient in nickel hydroxide particles used in NiMH batteries, varies by 3 orders of magnitude over the span of the state-of-charge [14] [15]. In addition, the diffusion coefficient in Na₃V₂(PO₄)₂F₃ active particles in SIB varies by 2 orders of magnitude, even though ex situ XRD patens reveal negligible variation in lattice parameters, a characteristic of solid-solution processes [16] [17].

While incorporating a concentration-dependent diffusion coefficient in electrochemical models leads to more accurate simulation results, the lack of an analytical exact solution method for Equation (1), for a general, nonlinear form of the diffusion coefficient, means one has to rely on numerical methods. [18] Traditionally, this is solved by the explicit and implicit finite difference methods (FDM) [6] [19]. However, explicit schemes are conditionally stable and therefore, computationally expensive. For a grid spacing $\text{Δ}x$, the von Neumann stability condition for the explicit scheme

$\text{Δ}t\le \frac{\text{Δ}{x}^2}{2D},$ (2)

imposes a strong time step restriction for particle sizes of the order of 0.5 - 5 µm in radius. In general, FDM do not conserve a perfect mass balance and this error is propagated to the overall battery simulation at long times [20]. Other grid-based methods, with a perfect mass conservation include, the finite element method (FEM) and the finite volume method (FVM). While both methods are renowned for their robustness, they have the inherent disadvantage of not calculating concentrations at specific node points, in particular, at the surface boundary [20] [21]. Instead, one obtains volume averaged concentrations within discrete volume elements. Despite the additional computations and approximation errors arising from approximating the surface concentration, these methods have been successfully applied to battery simulations, with variable diffusion coefficients.

To address the shortcomings of the FEM and FVM, Zeng et al. [21] proposed the control volume method (CVM). The CVM is a class of finite volume discretization, which computes concentrations at node points. Therefore, surface concentrations are obtained directly in the CVM. Compared to the FVM, the CVM has a higher accuracy for a given number of mesh points and is more suited for battery modeling. [21] While these authors focused on the Crank-Nicolson time domain discretization, which sometimes produces oscillatory solutions, we herein present a backward Euler control volume method (BECV) to resolve the spherical diffusion problem with a variable diffusion coefficient. This method incorporates all the advantages of the CVM, with the added advantage of being stable and easier to implement. Because obtaining a fully-implicit solution involves a series of iterative steps, a hybrid backward Euler control volume method (HBECV) is herein introduced for the first time. The HBECV method is based on the linearization of the functional form of the diffusion coefficient and obtains the implicit solution in a single iteration. Using the HBECV, computationally efficient, accurate and unconditionally stable results are obtained for the surface concentration. While this work focuses on the single particle diffusion problem, the solution extension to a full, physics-based battery model is trivial.

The challenge with all grid based methods is that, the number of states in a system of coupled differential equations dramatically increases, when the number of spatially distributed grid points increases [6] [21]. In order to reduce the number of grid points, while maintaining a high accuracy, it is recommended to implement a non-regular grid spacing, which allocates more points towards the surface boundary, where the concentration gradients will be steep initially. Therefore, we present the CVM for non-uniform grid spacing and further perform a mesh optimization study in which factors affecting optimum mesh spacing are identified. It is shown that the diffusion length is a primary factor, which determines the optimum grid spacing.

2. Model

The second order, time-dependent, partial differential equation for species diffusion inside a spherical particle, is governed by Fick’s second law which is expressed in Equation (1). Two Neumann-type boundary conditions are relevant for a battery simulation, at the surface and at the center of the particle. At the surface, due to interfacial electrochemical reactions, the rate of species transport across the surface is expressed as a flux. This boundary condition is written as

$-D\frac{\partial c}{\partial r}=J,\text{at}r=R,\forall t,$ (3)

where J is the interfacial flux of species [mol·m⁻²·s⁻¹] and R is the radius of the particle [m].

We define J as positive for species diffusing out of the particle. J can also be expressed in terms of the current density and its magnitude is the same at all points of the surface, i.e. it is uniform. In turn, this implies spherical symmetry.

At the center of the particle, due to the flux symmetry, the flux is zero

$-D\frac{\partial c}{\partial r}=0,\text{at}r=0,\forall t.$ (4)

To obtain the concentration profile inside the particle, we apply the CVM. For Neumann-type boundary conditions, the CVM performs better compared to the FDM because mass conservation breaks down in the FDM.

Since the diffusion problem of Equations (1), (3)-(4) is spherically symmetric, the magnitude of flux depends only on r. Consider a set of N discretization points on r, such that

${\left\{r_i\right\}}_i^N,\text{for}1\le i\le N$ (5)

where the index i defines node positions of CVM, and N is a non-zero natural number of grid points. Accordingly, $r_1=0$ and $r_N=R$ are the particle center and surface, respectively.

Each i^th node point in Equation (5) can be assigned a control volume element to it. For a spherical geometry, and for $2\le i\le N-1$, such a control volume element is a shell whose external faces (boundaries) are located halfway between adjacent nodes.

Let ${\text{IB}}_i$ and ${\text{OB}}_i$ define the i^th inner control volume boundary and outer control volume boundary, respectively, according to

${\text{IB}}_i=r_{i-1}+\frac{\Delta r_{i-1}}{2},\text{for}2\le i\le N-1$ (6)

${\text{OB}}_i=r_i+\frac{\Delta r_i}{2},\text{for}2\le i\le N-1$ (7)

where $\Delta r_i=r_{i+1}-r_i$ is the spacing between two adjacent node points. In this way, the i^th inner control volume element is a shell imbedding node point i.

At the boundaries, exceptions arise. Both ${\text{IB}}_1$ and ${\text{OB}}_N$ are not located between node points but right at node points $r_1$ and $r_N$, respectively. This implies, ${\text{IB}}_1=0$ and ${\text{OB}}_N=R$. On the other hand, ${\text{OB}}_1$ and ${\text{IB}}_N$ are located between adjacent node points and can be expressed as

${\text{OB}}_1=\frac{\Delta r_1}{2},$ (8)

${\text{IB}}_N=r_{N-1}+\frac{\Delta r_{N-1}}{2}.$ (9)

Figure 1 illustrates the CVM discretization. Black dots and solid black lines represent node points while white dots and dotted black lines represent control volume boundaries. A magnified view of the discretized particle illustrates the 3D nature of spherical shells arising from inner control volume discretization. Several important features of the control volume discretization should be noted:

1) Concentrations are calculated at the node points only. No concentrations are calculated at the boundaries between control volume elements.

2) The concentration profile between nodes is assumed to be linear.

3) Concentration gradients are calculated at spherical shell boundaries using concentration values from adjacent nodes.

4) Interior boundaries of the control volume shells are located halfway between adjacent nodes.

Let $v_i$ denote the control volume element at node point i. Assume $m_i\left(t\right)$, the amount of electrochemically active species [mol] inside $v_i$ at arbitrary time t, can be expressed as a product of the concentration at node point i and volume of corresponding spherical shell. This means

Figure 1. Diffusion in a spherical particle illustrating the control volume discretization along the particle radius. The solid black lines and solid black dots represent grid point i, while the dotted black line and white dots represent control volume boundaries. In the magnified view, i is surrounded by an imaginary control volume between the outer boundary ${\text{OB}}_i$ and inner boundary ${\text{IB}}_i$. The flux is defined positive for species diffusing out of the particle.

$m_i\left(t\right)={\displaystyle \underset{v_i}{\int }c\left(r\right)\text{d}V\left(r\right)}\approx c\left(r_i\right){\displaystyle \underset{v_i}{\int }\text{d}V\left(r\right)}=c\left(r_i\right)U_i,$ (10)

where $U_i$ is the total volume of $v_i$ [m³].

Mass conservation law in the absence of source term(s) dictates that any change in $m_i$ corresponds to the net-flux via control of the volume inner and outer boundaries, i.e.

${\Delta }_t{m}_i\left(t\right)={\Delta }_t{c}_i{U}_i=\left(J_i{A}_i-J_{i-1}{A}_{i-1}\right)\Delta t,\text{for}1\le i\le N$ (11)

where $J_i$ represents mass flux though control volume boundaries [mol·m⁻²·s⁻¹], $\Delta t$ is the time step [s] and the symbol ${\Delta }_t$ denotes a change of a variable in time. $A_i$ is the surface area of the i^th control volume boundary [m²] which is defined as

$A_i=4\pi {\left(r_i+\frac{\Delta r_i}{2}\right)}^2,\text{for}1\le i\le N-1,$ (12)

at the surface boundary ${\text{OB}}_N$ as

$A_N=4\pi R^3,$ (13)

and at the center boundary ${\text{IB}}_1$ as

$A_0=0.$ (14)

Note that, the absence of source term(s) in Equation (11) is due to the lack of internal species production or consumption within the active particles.

Now, according to remarks ii and iii, the fluxes at the control volume boundaries can be derived as follows

$-D\left(c\right){\frac{\partial c}{\partial r}|}_{r=r_i+\frac{\Delta r_i}{2}}\approx -D_{i+1/2}\frac{c_{i+1}-c_i}{r_{i+1}-r_i}=-D_{i+1/2}\frac{c_{i+1}-c_i}{\Delta r_i},$ (15)

where $D_{i+1/2}=D\left(\frac{c_{i+1}+c_i}{2}\right)$ is the concentration-dependent diffusion coefficient at the control volume boundary [m²·s⁻¹]. A half-sum is applied because of remarks ii and iv. Substituting Equation (12) and (15) into Equation (11), gives

$\begin{array}{l}-D_{i-\frac{1}{2}}\frac{c_i-c_{i-1}}{\Delta r_{i-1}}{\left(r_{i-1}+\frac{\Delta r_{i-1}}{2}\right)}^2\Delta t+D_{i+\frac{1}{2}}\frac{c_{i+1}-c_i}{\Delta r_i}{\left(r_i+\frac{\Delta r_i}{2}\right)}^2\Delta t\\ =\frac{U_i}{4\pi }{\Delta }_t{c}_i,\text{for}2\le i\le N-1.\end{array}$ (16)

In order to eliminate the factor 4π from subsequent derivations, let a normalized volume $V_i$ be defined as

$V_i=\frac{1}{3}\left[{\left(r_i+\frac{\Delta r_i}{2}\right)}^3-{\left(r_{i-1}+\frac{\Delta r_{i-1}}{2}\right)}^3\right],\text{for}2\le i\le N-1,$ (17)

at the center ( $i=1$ ) and at the surface boundary ( $i=N$ ), $V_1$ and $V_N$ are defined as

$V_1=\frac{1}{24}\Delta r_1^3,$ (18)

$V_N=\frac{1}{3}\left[R^3-{\left(R-\frac{\Delta r_{N-1}}{2}\right)}^3\right].$ (19)

To further economize notations, let variable $K_i$ [m³] be introduced

$K_i=D_{i+1/2}\frac{\Delta t}{\Delta r_i}{\left(r_i+\frac{\Delta r_i}{2}\right)}^2,\text{for}1\le i\le N-1.$ (20)

Therefore, taking Equations (17) and (20) into account, Equation (16) can be expressed in terms of $K_i$ and $V_i$ as

$-K_{i-1}\left(c_i-c_{i-1}\right)+K_i\left(c_{i+1}-c_i\right)=V_i{\Delta }_t{c}_i,\text{for}2\le i\le N-1.$ (21)

Rearranging Equation (21) we finally obtain

$K_{i-1}{c}_{i-1}-\left(K_{i-1}+K_i\right)c_i+K_i{c}_{i+1}=V_i{\Delta }_t{c}_i,\text{for}2\le i\le N-1.$ (22)

Equation (22) is defined at interior node points. It is possible, starting with the general mass balance expression of Equation (11) and following the steps shown in Equations (15)-(22), to obtain expressions for the two remaining boundary cases, $i=1$ and $i=N$.

At the center (at $i=1$ ), there is zero flux through ${\text{IB}}_1$ according to Equation (4). Furthermore, the surface area at ${\text{IB}}_1$ is zero according to Equation (14). Applying the general mass balance on $v_1$ gives

$D_{1+1/2}\frac{c_2-c_1}{\Delta r_1}{\left(r_1+\frac{\Delta r_1}{2}\right)}^2\Delta t=V_1{\Delta }_t{c}_1,$ (23)

Finally, Equation (23) expressed in terms of the variable $K_i$, becomes

$K_1\left(c_2-c_1\right)=V_1{\Delta }_t{c}_1.$ (24)

At the surface (at $i=N$ ), there is a uniform interfacial flux J. Applying the boundary condition of Equation (3) and the general mass balance on $v_N$, we obtain

$-D_{N-1/2}\frac{c_N-c_{N-1}}{\Delta r_{N-1}}{\left(r_{N-1}+\frac{\Delta r_{N-1}}{2}\right)}^2\Delta t-J{R}^2\Delta t=V_N{\Delta }_t{c}_N,$ (25)

Finally, introducing the variable $K_i$ into Equation (25) gives

$-K_{N-1}\left(c_N-c_{N-1}\right)-J{R}^2\Delta t=V_N{\Delta }_t{c}_N.$ (26)

Up to this point, the temporal discretization is intentionally omitted because the system of equations, Equations (22), (24) and (26) can be solved either by the forward Euler or the backward Euler method.

Let the superscript j represent the current time step and $j-1$ represent the previous time step. Therefore, ${\Delta }_t{c}_i=c_i\left(t\right)-c_i\left(t-\text{Δ}t\right)=c_i^j-c_i^{j-1}$. The system of equations, Equations (22), (24) and (26) is thus expressed in the backward Euler scheme as

$-K_{i-1}^j{c}_{i-1}^j+\left(K_{i-1}^j+K_i^j+V_i\right)c_i^j-K_i^j{c}_{i+1}^j=V_i{c}_i^{j-1},\text{for}2\le i\le N-1$ (27)

at the center

$\left(K_1^j+V_1\right)c_1^j-K_1^j{c}_2^j=V_1{c}_1^{j-1},$ (28)

and at the surface

$-K_{N-1}^j{c}_{N-1}^j+\left(K_{N-1}^j+V_N\right)c_N^j=V_N{c}_N^{j-1}-J{R}^2\Delta t.$ (29)

Equations (27)-(29) represent a coupled system of equations since all values at time step j are unknown while values at time step $j-1$ are unknown.

2.1. Solving the Coupled System of Equations

As a first step to finding the solution, the coupled system of Equations (27)-(29) is expressed in matrix form. This is expressed as

$M{c}^j=V{c}^{j-1}-J{R}^2\Delta t,$ (30)

where $c^j$ is a column vector containing concentrations at all node points at time index j, i.e. $c^j=\left(c_1^j,c_2^j,\cdots ,c_N^j\right)$, $M$ is an N-by-N matrix

$\left(\begin{array}{cccccccc}{K}_1^j+V_1& -K_1^j& 0& 0& \cdots & 0& 0& 0\\ -K_1^j& K_1^j+K_2^j+V_2& -K_2& & \cdots & 0& 0& 0\\ 0& -K_2^j& K_2^j+K_3^j+V_3& -K_3^j& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -K_{N-2}^j& K_{N-2}^j+K_{N-1}^j+V_{N-1}& -K_{N-1}^j\\ 0& 0& 0& 0& \cdots & 0& -K_{N-1}^j& K_{N-1}^j+V_N\end{array}\right)$

$V$ is expressed as

$V=\left(\begin{array}{cccc}{V}_1& 0& \cdots & 0\\ 0& V_2& & \\ ⋮& & \ddots & \\ 0& 0& \cdots & V_N\end{array}\right)$

and $J$ is a (column) vector, ${\left(0,0,\cdots ,0,J\right)}^{\prime }$, whose only non-zero entry is at the N^th point, corresponding to the surface.

$M$ is a tridiagonal matrix. If $V_i>0$, a condition which is trivially satisfied by the construction of sequence r as shown in Equation (5), $M$ is strictly (row) diagonally dominant. By the Levy-Desplanques theorem, $M$ is non-singular and therefore invertible [22]. The tridiagonal matrix algorithm (TDMA) can thus be applied to solve Equation (30) as a stable and fast solution method [23]. The TDMA, also known as the Thomas algorithm, is a variant of Gaussian elimination, which is applicable to a diagonally dominant tridiagonal system of N unknowns. Compared to the standard Gaussian elimination or matrix inversion, which require $O\left(N^3\right)$ operations to solve, the TDMA only requires $O\left(N\right)$ operations [24].

For a constant diffusion coefficient, Equation (30) is linear. The solution $c^j$ is rapidly obtained in a single run of the TDMA. However, for the problem posed above, the diffusion coefficient at time step j is not known a priori since D is an implicit function of $c^j$. Equation (23) therefore takes form

$M\left(c^j\right)c^j=V{c}^{j-1}-J{R}^2\Delta t,$ (31)

and represents a non-linear system of equations. Equation (31) can be iteratively solved by various fixed point methods such as Newton’s method, Jacobi, line-by-line and Gauss-Seidel. These iterative methods can be used in combination with the TDMA to obtain convergent solutions [23] [25]. The Newton’s method in particular, attains quadratic convergence if the initial guess is sufficiently close to the solution [21] [25]. Nevertheless, great care must be taken to correctly construct and solve the Newton equation and a poor initial guess may even result in a lack of convergence [25]. In this work, the following iterative scheme (Jacobi) is used

$M\left(c_{k-1}^j\right)c_k^j=V{c}^{j-1}-J{R}^2\Delta t,\text{for}1\le k\le k_{tot}$ (32)

where subscript k denotes the iteration number in a total of $k_{tot}$ iterations. $c_k^j$ is therefore the implicit solution obtained after k iterations. The initial value needed for the first iteration, is defined as

$c_0^j=c^{j-1}.$ (33)

For $k_{tot}\approx 20$, the fully-implicit, BECV solution is obtained. It shall be demonstrated that, due to the initial condition Equation (33), the first iteration of Equation (32) achieves a stable and approximately accurate solution, which is acceptable in many cases. This solution is herein referred to as the HBECV. The HBECV is therefore a linearization of the concentration-dependent $M$, in order to obtain implicit solutions in a single iteration.

2.2. Grid Spacing

In order to accurately determine concentration profiles, many grid points are required. However, more grid points come at considerable computational costs. The accuracy of the CVM, for an economical number of grid points, depends on the spatial distribution of the same. The choice of grid spacing, however, depends on the nature of the problem and boundary conditions. More points are required at the regions where the concentration profile has steep gradients and this holds true at the particle surface boundary.

While the scheme of equations presented in this work, allows for variable spacing of grid points, the literature around grid/mesh optimization is very sparse, this in turn makes it difficult to comprehend the principle factors affecting the optimum grid spacing. In order to evaluate the different grid-point locations, the following geometric spacing equation is applied

$r_i=R\left(1-\frac{Y^{\frac{N-i}{N-1}}-1}{Y-1}\right),1\le i\le N,$ (34)

where Y is the common factor of the geometric series. In this study, Y varies between 2 and 20. If $Y=10$, the logarithmic spacing is obtained, a notable choice in the preceding publication [21]. To evaluate the error in each value of Y, as a function of D and R, a solution obtained from a linear spaced grid of 501 points and $dt=5$ s is used as reference solution.

3. Results and Discussion

To investigate the accuracy of HBECV, the set of parameters from Zeng et al. [21] is used. Table 1 shows the parameters for Li(Ni_1/3Mn_1/3Co_1/3)O₂ particles used in the referenced work and also in this work.

The concentration-dependent diffusion coefficient for Li(Ni_1/3Mn_1/3Co_1/3)O₂ used by Wu et al. (see ref. [26]) is

$D\left(c\right)=D^{ref}\left[1+100{\left(\frac{{\hat{C}}_{theo}\left(c_{\max }-c\right)}{{\hat{C}}_{prac}{c}_{\max }}\right)}^{\frac{3}{2}}\right],$ (35)

where $D^{ref}$ is defined as the reference diffusion coefficient of 2 × 10⁻¹⁶ m²·s⁻¹, ${\hat{C}}_{theo}$ is the theoretical capacity of the electrode material (277.84 mAh·g⁻¹) and ${\hat{C}}_{prac}$ is the practical capacity of the electrode material of 160 mAh·g⁻¹.

Figure 2 illustrates the agreement between our results and literature. After a discharge time t = 400 s, the surface concentration nearly reaches c_max and thus the end of discharge. For the parameters in Table 1, a dense mesh of 501 grid points is needed to eliminate errors due to spatial discretization. These results validate the BECV method whose solutions are obtained after $k_{tot}=20$ iterations of the Jacobi method.

It is interesting to evaluate the effect of reducing k_tot because the computation speed increases with less iteration. Figure 3 shows that the relative error in surface concentration progressively increases when k_tot decreases. Furthermore, approximately 10 iterations of the Jacobi method are necessary to obtain a fully-implicit BECV solution. With regards to fast simulations, the HBECV solution obtained when $k_{tot}=1$ is interesting. From the results, the relative error on the surface concentration of the HBECV method is approximately 0.1%, which is good enough for practical purposes and moreover the results are stable. This validates the HBECV and justifies its use for fast, stable and practically accurate results.

Table 1. Parameters for modeling diffusion in spherical particle with variable diffusion coefficient.

Figure 2. Comparison of the results of the simulated concentration gradients obtained in this work, using the BECV method and results from Zeng et al. [21] Results are obtained using 501 uniformly spaced grid points, using the parameters of Table 1 and the concentration-dependent diffusion coefficient expression, Equation (35).

Figure 3. Relative error in surface concentration over 400 s simulation as a function of the number of iterations. The relative error of the HBECV method is compared to the iterative implicit BECV method. 501 uniformly spaced grid points and parameters in Table 1 are used to calculate the surface concentration. As k_tot, the total number of iterations per time-step increases, the solution converges to the reference solution of 20 iterations.

A uniform grid of 501 points is nevertheless impractical for use in battery simulations where the number of grid points is limited. The first step in reducing the number of grid points is to determine optimal grid spacing. This is herein determined by changing the value of Y in a model with 301 grid points. The optimization of Y is performed using the uniform grid of 501 points as the reference solution, while the cost function is the normalized root of squared deviations between these solutions.

Figure 4 illustrates the dependence of optimum Y ( $Y_{opt}$ ) on the dimensionless parameter $\frac{\sqrt{4D^{ref}T}}{R}$. This parameter represents a scaled diffusion length since the mass diffusion length is known to be equal to $\sqrt{4D^{ref}T}$. Two regions can be easily identified in Figure 4.

1) high diffusion length region where $\frac{\sqrt{4D^{ref}T}}{R}>0.1$ and $3\le Y_{opt}\le 6$,

2) low diffusion length region where $\frac{\sqrt{4D^{ref}T}}{R}\le 0.1$ and $Y_{opt}>6$.

For a high diffusion length, which is more than 10% of R, an evenly spaced grid defined by low Y values should be used. On the other hand, for a low diffusion length, which is less than 10% of R, a logarithmic grid spacing, defined by high Y values, is more appropriate. The immediate conclusion is that more points are needed close to the surface only when the diffusion length is low. This highlights the importance of carefully selecting the grid spacing for a given $D^{ref}$ and R values and not rely on an intuitive feeling of the grid spacing.

For the parameters listed in Table 1, $Y_{opt}=12$ and the number of grid points distributed by $Y_{opt}$ is defined as $N_{opt}$. It is important to further reduce $N_{opt}$, to a practical number, relevant to full-cell battery modeling. Figure 5 shows the relative error on surface concentration at $t=400\text{s}$ and $dt=0.1\text{s}$, as a function of $N_{opt}$. As expected, the error relative to the fine grid mesh increases as $N_{opt}$ decreases. However, the BECV and the HBECV remarkably converge to the same error when $N_{opt}$ decreases. This implies is that, the spatial discretization error exceeds the linearization error of the HBECV when the number of optimally spaced grid points is below 21. Therefore, based on the results shown in Figure 5, the HBECV should be used in battery simulations, instead of the more computationally expensive BECV.

Figure 4. Grid optimization, effect of the geometric factor Y and the dimensionless factor $\frac{\sqrt{4D^{ref}T}}{R}$. The red squares show the optimum geometric factor, $Y_{opt}$ for selected $D^{ref}$ and R values. The bars represent variance within 1% of $Y_{opt}$. Inset showing the distribution of 6 grid points for $Y=2$, and $Y=10$.

Figure 5. Relative error on surface concentration at 400 s of the BECV and the HBECV as a function of $N_{opt}$, the number of optimized grid points. The optimum geometric factor, $Y_{opt}=12$ and time step, $dt=0.1\text{s}$ are used. The reference solution is obtained from 501 uniform grid points and $dt=0.1\text{s}$.

4. Conclusions

In this paper, the backward Euler control volume (BECV) and the hybrid backward Euler control volume (HBECV) methods are presented as efficient numerical methods to resolve the solid-state spherical diffusion problem for a variable diffusion coefficient. The implicit scheme of nonlinear equations is herein shown to be strictly diagonally dominant and efficiently solved by the tridiagonal matrix algorithm (TDMA).

The fully-implicit BECV is shown to require more computations, approximately 10 iterations of the TDMA, in order to reach a convergent solution; however, the HBECV only requires one iteration. The HBECV achieves this by linearization of the implicit scheme of equations. Although, in comparison to the BECV method, the error in surface concentration in the HBECV method is around 0.1%, for a fine grid of 501 points, the error difference decreases when the number of grid points decreases. For a course grid of 21 optimally spaced grid points, the error in surface concentrations converges to the same order of magnitude, which demonstrates that the error due to spatial discretization outweighs the error due to the linearization, as introduced by the HBECV method.

This work further explores the parameters governing optimal grid spacing. The diffusion length emerged as a guiding parameter for selecting optimum grid spacing. It is found that low diffusion length problems require more points distributed near the surface, compared to high diffusion length problems wherein an evenly spaced grid more appropriate.

As demonstrated in this work, the HBECV is accurate, stable and easy to implement. In practical battery simulations, where the number of grid points is minimized, the HBECV performs as well as the BECV method. Therefore, given the aforementioned advantages of the HBECV, this method is a justified choice for modeling the second order partial differential equation of solid-state diffusion with a concentration-dependent diffusion coefficient.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Eisler, M.N. (2016) Cold War Computers, California Supercars, and the Pursuit of Lithium-Ion Power. Physics Today, 69, 30-36. https://doi.org/10.1063/PT.3.3296
[2]	Whittingham, M.S. (1976) Electrical Energy Storage and Intercalation Chemistry. Science, 192, 1126-1127. https://doi.org/10.1126/science.192.4244.1126
[3]	NobelPrize.org. Nobel Media AB (2019) The Nobel Prize in Chemistry 2019. https://www.nobelprize.org/prizes/chemistry/2019/press-release/
[4]	Liu, S. (2006) An Analytical Solution to Li/Li+ Insertion into a Porous Electrode. Solid State Ionics, 177, 53-58. https://doi.org/10.1016/j.ssi.2005.09.053
[5]	Chayambuka, K., Mulder, G., Danilov, D.L. and Notten, P.H.L. (2019) A Modified Pseudo-Steady-State Analytical Expression for Battery Modeling. Solid State Communications, 296, 49-53. https://doi.org/10.1016/j.ssc.2019.04.011
[6]	Ramadesigan, V., Boovaragavan, V., Pirkle, J.C. and Subramanian, V.R. (2010) Efficient Reformulation of Solid-Phase Diffusion in Physics-Based Lithium-Ion Battery Models. Journal of The Electrochemical Society, 157, A854-A860. https://doi.org/10.1149/1.3425622
[7]	Subramanian, V.R., Diwakar, V.D. and Tapriyal, D. (2005) Efficient Macro-Micro Scale Coupled Modeling of Batteries. Journal of The Electrochemical Society, 152, A2002-A2008. https://doi.org/10.1149/1.2032427
[8]	Subramanian, V.R., Ritter, J.A. and White, R.E. (2001) Approximate Solutions for Galvanostatic Discharge of Spherical Particles I. Constant Diffusion Coefficient. Journal of The Electrochemical Society, 148, E444-E449. https://doi.org/10.1149/1.1409397
[9]	Zhang, Q. and White, R.E. (2007) Comparison of Approximate Solution Methods for the Solid Phase Diffusion Equation in a Porous Electrode Model. Journal of Power Sources, 165, 880-886. https://doi.org/10.1016/j.jpowsour.2006.12.056
[10]	Mao, J., Tiedemann, W. and Newman, J. (2014) Simulation of Temperature Rise in Li-Ion Cells at Very High Currents. Journal of Power Sources, 271, 444-454. https://doi.org/10.1016/j.jpowsour.2014.08.033
[11]	Kazemi, N., Danilov, D.L., Haverkate, L., Dudney, N.J., Unnikrishnan, S. and Notten, P.H.L. (2019) Modeling of All-Solid-State Thin-Film Li-Ion Batteries: Accuracy Improvement. Solid State Ionics, 334, 111-116. https://doi.org/10.1016/j.ssi.2019.02.003
[12]	Raijmakers, L.H.J., Danilov, D.L., Eichel, R.-A. and Notten, P.H.L. (2020) An Advanced All-Solid-State Li-Ion Battery Model. Electrochimica Acta, 330, 135147. https://doi.org/10.1016/j.electacta.2019.135147
[13]	Christensen, J. and Newman, J. (2006) Stress Generation and Fracture in Lithium Insertion Materials. Journal of Solid State Electrochemistry, 10, 293-319. https://doi.org/10.1007/s10008-006-0095-1
[14]	Motupally, S., Streinz, C.C. and Weidner, J.W. (1995) Proton Diffusion in Nickel Hydroxide Films Measurement of the Diffusion Coefficient as a Function of State of Charge. Journal of the Electrochemical Society, 142, 1401-1408. https://doi.org/10.1149/1.2048589
[15]	Motupally, S., Streinz, C.C. and Weidner, J.W. (1998) Proton Diffusion in Nickel Hydroxide Prediction of Active Material Utilization. Journal of The Electrochemical Society, 145, 29-34. https://doi.org/10.1149/1.1838205
[16]	Shakoor, R.A., Seo, D.-H., Kim, H., Park, Y.-U., Kim, J., Kim, S.-W., Gwon, H., Lee, S. and Kang, K. (2012) A Combined First Principles and Experimental Study on Na3V2(PO4)2F3 for Rechargeable Na Batteries. Journal of Materials Chemistry, 22, 20535-20541. https://doi.org/10.1039/c2jm33862a
[17]	Liu, Z., Hu, Y.-Y., Dunstan, M.T., Huo, H., Hao, X., Zou, H., Zhong, G., Yang, Y. and Grey, C.P. (2014) Local Structure and Dynamics in the Na Ion Battery Positive Electrode Material Na3V2(PO4)2F3. Chemistry of Materials, 26, 2513-2521. https://doi.org/10.1021/cm403728w
[18]	Ford Versypt, A.N. and Braatz, R.D. (2014) Analysis of Finite Difference Discretization schemes for DIFFUSION in Spheres with Variable Diffusivity. Computers & Chemical Engineering, 71, 241-252. https://doi.org/10.1016/j.compchemeng.2014.05.022
[19]	Langtangen, H.P. (2013) Finite Difference Methods for Diffusion Processes. University of Oslo, Oslo.
[20]	Urisanga, P.C., Rife, D., De, S. and Subramanian, V.R. (2015) Efficient Conservative Reformulation Schemes for Lithium Intercalation. Journal of The Electrochemical Society, 162, A852-A857. https://doi.org/10.1149/2.0061506jes
[21]	Zeng, Y., Albertus, P., Klein, R., Chaturvedi, N., Kojic, A., Bazant, M.Z. and Christensen, J. (2013) Efficient Conservative Numerical Schemes for 1D Nonlinear Spherical Diffusion Equations with Applications in Battery Modeling. Journal of The Electrochemical Society, 160, A1565-A1571. https://doi.org/10.1149/2.102309jes
[22]	Taussky, O. (1949) A Recurring Theorem on Determinants. The American Mathematical Monthly, 56, 672-676. https://doi.org/10.1080/00029890.1949.11990209
[23]	Patankar, S. (2018) Numerical Heat Transfer and Fluid Flow. CRC Press, Boca Raton. https://doi.org/10.1201/9781482234213
[24]	Higham, N.J. (1986) Efficient Algorithms for Computing the Condition Number of a Tridiagonal Matrix. SIAM Journal on Scientific and Statistical Computing, 7, 150-165. https://doi.org/10.1137/0907011
[25]	Conte, S.D. and De Boor, C. (2017) Elementary Numerical Analysis: An Algorithmic Approach. SIAM, Philadelphia. https://doi.org/10.1137/1.9781611975208
[26]	Wu, S.-L., Zhang, W., Song, X., Shukla, A.K., Liu, G., Battaglia, V. and Srinivasan, V. (2012) High Rate Capability of Li(Ni1/3Mn1/3Co1/3)O2 Electrode for Li-Ion Batteries. Journal of the Electrochemical Society, 159, A438-A444. https://doi.org/10.1149/2.062204jes