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

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.


Introduction
Since the discovery of the electrochemical intercalation reaction of lithium in layered titanium disulfide (TiS 2 ) by Whittingham in 1976, lithium-ion batteries (LIB) have grown in popularity and application to surpass all other electrical 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 tration. 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 re-quirements 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 3 V 2 (PO 4 ) 2 F 3 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 Δx , the von Neumann stability condition for the explicit scheme 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.

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 where J is the interfacial flux of species [mol·m −2 ·s −1 ] and R is the radius of the 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 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 where the index i defines node positions of CVM, and N is a non-zero natural number of grid points. Accordingly, 1 0 r = and N r 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 1 i N ≤ ≤ − , such a control volume element is a shell whose external faces (boundaries) are located halfway between adjacent nodes.
Let IB i and OB i define the i th inner control volume boundary and outer control volume boundary, respectively, according to where 1 i i i r r r + ∆ = − 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 1 IB and OB N are not located between node points but right at node points 1 r and N r , respectively. This implies, On the other hand, 1 OB and IB N are located between adjacent node points and can be expressed as 2) The concentration profile between nodes is assumed to be linear. . 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 OB i and inner boundary IB i .
The flux is defined positive for species diffusing out of the particle.
Mass conservation law in the absence of source term(s) dictates that any change in i m corresponds to the net-flux via control of the volume inner and outer boundaries, i.e.
at the surface boundary OB N as and at the center boundary 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 ( )  (12) and (15) into Equation (11), gives Journal of Applied Mathematics and Physics In order to eliminate the factor 4π from subsequent derivations, let a normalized volume i V be defined as at the center ( 1 i = ) and at the surface boundary ( i N = ), 1 V and N V are defined as To further economize notations, let variable i K [m 3 ] be introduced Therefore, taking Equations (17) and (20) into account, Equation (16) can be expressed in terms of i K and i V as Rearranging Equation (21) we finally obtain ( ) 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, 1 i = and i N = .
At the center (at 1 i = ), there is zero flux through 1 IB according to Equation (4). Furthermore, the surface area at 1 IB is zero according to Equation (14).
Applying the general mass balance on 1 v gives Finally, Equation (23) expressed in terms of the variable i K , becomes ( ) 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 1 j − represent the previous time step. Therefore, The system of equations, Equations (22), (24) and (26) is thus expressed in the backward Euler scheme as ( ) at the center ( ) and at the surface ( ) Equations (27)

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 and J is a (column) vector, ( ) 0, 0, , 0, J ′  , whose only non-zero entry is at the N th point, corresponding to the surface.
M is a tridiagonal matrix. If 0 i V > , 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

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 where Y is the common factor of the geometric series. In this study, Y varies between 2 and 20. If 10 Y = , 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 5 dt = s is used as reference solution.

Results and Discussion
To investigate the accuracy of HBECV, the set of parameters from Zeng et al. [21] is used. Table 1   is the practical capacity of the electrode material of 160 mAh·g −1 . 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 20 tot k = 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 1 tot k = 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.  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 refer-ence solution, while the cost function is the normalized root of squared deviations between these solutions. 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 For the parameters listed in Table 1, 12 opt Y = and the number of grid points distributed by opt Y is defined as opt N . It is important to further reduce opt N , to a practical number, relevant to full-cell battery modeling. Figure 5 shows the relative error on surface concentration at , as a function of opt N . As expected, the error relative to the fine grid mesh increases as opt N decreases. However, the BECV and the HBECV remarkably converge to the same error when opt N 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.

Conclusions
In this paper, the backward Euler control volume (BECV) and the hybrid back- 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.