Abstract
Lithium-ion batteries are typically modeled using porous electrode theory coupled with various transport and reaction mechanisms with an appropriate discretization or approximation for the solid phase. One of the major difficulties in simulating Li-ion battery models is the need for simulating solid-phase diffusion in a second dimension 
. It increases the complexity of the model as well as the computation time/cost to a great extent. Traditional approach toward solid-phase diffusion leads to more difficulties, with the use of emerging cathode materials, which involve phase changes and thus moving boundaries. A computationally efficient representation for solid-phase diffusion is discussed in this paper. The operating condition has a significant effect on the validity, accuracy, and efficiency of various approximations for the solid-phase diffusion. This paper compares approaches available today for solid-phase reformulation and provides two efficient forms for constant and varying diffusivities in the solid phase. This paper introduces an efficient method of an eigenfunction based Galerkin collocation and a mixed order finite difference method for approximating/representing solid-phase concentration variations within the active materials of porous electrodes for a pseudo-two-dimensional model for lithium-ion batteries.
Export citation and abstract BibTeX RIS
Electrochemical power sources are expected to play a vital role in the future in automobiles, power storage, military, mobile, and space applications. Lithium-ion chemistry has been identified as a good candidate for high power/high energy secondary batteries. Significant progress has been made toward modeling and understanding of lithium-ion batteries using physics-based first-principles models. First-principles-based battery models typically solve electrolyte concentration, electrolyte potential, solid-state potential, and solid-state concentration in the porous electrodes1, 2 and electrolyte concentration and electrolyte potential in the separator. These models are based on transport phenomena, electrochemistry, and thermodynamics. These models are represented by coupled nonlinear partial differential equations (PDEs) in one to two dimensions and are typically solved numerically and require a few minutes to hours to simulate.
Even when one-dimensional transport in the macroscale 
is considered, the continuum models that are used to describe the electrochemical behavior of lithium-ion batteries, involve three coupled nonlinear PDEs (in 
) in each porous electrode and two coupled nonlinear PDEs (in 
) in the separator. For predicting the thermal behavior, one has to add an additional equation for temperature. In addition, solid-state diffusion should be solved in the pseudo second dimension 
in the electrode. Li-ions diffuse (intercalate) into and out of the solid particles of porous electrodes in the pseudo second dimension. Hence, in addition to the equations in the 
direction, solid-state diffusion should be solved in the pseudo dimension 
in the porous electrodes. This diffusion in the microscale is typically modeled using Fick's law of diffusion. One of the major difficulties in the electrochemical engineering models is the inclusion of solid-phase diffusion in a second dimension 
, which increases the complexity of the model as well as the computation time/cost to a great extent. Traditional simulation approaches toward solid-phase diffusion lead to more difficulties with the use of emerging cathode materials, which involve phase changes and thus moving boundaries. Concentration variations in the solid phase is governed by Fick's law of diffusion, and the same in spherical coordinates is given as
where 
. Equation 1 can be converted to a dimensionless form using the following dimensionless variables and parameters
with the boundary conditions
This paper discusses two computationally efficient representations for the solid-phase diffusion. An efficient eigenfunction based Galerkin collocation method is introduced and discussed in the paper. Further, a mixed order finite difference (FD) method with optimal node spacing is introduced, which can be used to reduce the computational cost/time significantly even with varying diffusivities in the solid phase. The operating condition has a significant effect on the validity, accuracy, and efficiency of various approximations for the solid-phase diffusion. The discretization and solver scheme used in time is also a significant factor in deciding the best possible approximation for the solid phase. This paper also compares various methods1–6 for approximating/representing solid-phase concentration variations within the active materials of porous electrodes for a full-order pseudo-two-dimensional (2D) model for lithium-ion batteries. A comparison among these available methods along with a brief mention about their merits and usability is made to identify the best possible method and incorporate in a full-order pseudo-2D model.1
Existing Approximations and the Need for Efficient Reformulation
Porous electrode models of Li-ion batteries often use approximations to eliminate the time-consuming calculation in the second dimension 
for the solid-phase diffusion. These methods include the Duhamel's superposition method,1 diffusion length method,3 the polynomial approximation method,4 the pseudosteady-state (PSS) approach by Liu5 and the penetration depth analysis and mixed order finite element approach.6
Each of the above listed methods has its own advantages and disadvantages when used in Li-ion battery models. The following section gives a brief summary about each of the methods and discusses their merits and usability.
Duhamel's superposition method
The Duhamel's superposition method1 is a robust method available for representing the solid-phase diffusion for constant diffusivities. This method relates the solution of a boundary value problem with time-dependent boundary conditions to the solution of a similar problem with time-independent boundary conditions by a simple relation. More information about the method and equations are presented elsewhere.1, 2
Duhamel's superposition method is a robust method and is valid for any kind of operating condition, such as high rates of discharge, pulse power, etc. One of the major drawbacks of this method is that it cannot be used in DASSL-like solvers, which do not accept equations that are discretized in time and might as well be time-consuming for very stiff problems depending on the time steps taken. In addition, it cannot be used for nonlinear diffusivities.
Diffusion length method
The diffusion length method's approach3 is based on a parabolic profile approximation for the solid phase. The diffusion length method predicts that the surface concentration and volume-averaged concentration inside a particle are linearly dependent on each other, which should be valid only after the diffusion layer builds up to its steady state. Therefore, the method is inadequate at short times or under dynamic operations, such as pulse or current interrupt operations. The prediction based on the diffusion length method is inadequate at short times and is very efficient at long times and low rates.
Polynomial approximation
The polynomial approximation method by Subramanian et al.4 is based on parabolic profile approximation and volume averaging of the solid-phase diffusion equation. This high order polynomial method uses a different approach from the diffusion length method to improve the solution accuracy at short times. The diffusion length method uses the empirical exponential term in the equation and determines the multiplier value by matching surface concentration profiles to the exact solutions. The high order polynomial method uses a higher order polynomial for the concentration profile in the derivation, and one could derive new sets of equations with an even higher order polynomial model, if needed, following the same procedures discussed in the papers.4, 7
This method is very efficient at long time ranges and for low/medium rates and is ideal for adaptive solvers for pseudo-2D models. However, it is inaccurate at short times and for high rates/pulses and hence would not be a suitable method for implementing in models for Hybrid Electric Vehicles (HEVs) and other high rate applications.
PSS method
The PSS approach by Liu5 is very robust and, by having enough equations, this approach can cover the entire spectrum of high/low rates, pulses, etc. This is a form of a finite integral transform technique to eliminate the independent spatial variable 
from the solid-phase diffusion equation. For diffusion problems with a time-dependent pore wall flux, 
appears in the boundary condition and in the calculation of coefficients.
However, this method involves terms/coefficients, which blow up when the number of terms increases adding numerical difficulties for simulation. More details are given in the Results and Discussion section where this method is compared with our proposed approach implemented in this paper.
Penetration depth method
The penetration depth analysis approach was used earlier with empirical fits to the numerical solution for penetration depth near the surface of the particle. The advantage of this method is that it is very accurate at short times/pulses and more accurate and efficient penetration depth solutions can be directly obtained from the PDE as discussed elsewhere.6
The drawback with this approach is the need to be reinitialized every time and that it does not give a good match for varying δ values. Though this method is very accurate and efficient at short times, it is not ideal for adaptive solvers in a pseudo-2D model (increases stiffness).
Finite element method
While the governing equation 1 describes solid-phase concentration along the radius of each spherical particle of the active material, the macroscopic model requires only the concentration at the surface, 
, as a function of the time history of local reaction current density, 
. The PDE is transformed from spherical to planar coordinates and is discretized in the 
direction with 
suitably chosen linear elements. They used five elements with node points placed at 
. Transformed back to spherical coordinates, the discretized system is represented as ordinary differential equations in state space form and then solved.6
The finite element node sizes were probably obtained by trial and error and may not be optimal at long times or different operating conditions. The following section describes two methods that can be used for solid-phase diffusion approximation and explains the derivation of the same.
Galerkin Reformulation of Solid-Phase Diffusion
The reformulation discussed here is based on eigenfunction based Galerkin collocation for constant diffusivity. In constant diffusivity, 
in Eq. 6, 9 is 1. The following section describes the derivation for the eigenfunction based Galerkin collocation, and the equations for approximation are given at the end. Expanding Eq. 6 can be written as
A general solution for the variable 
can be written as a polynomial approximation given as
where 
are the eigenvalues of the given problem (Eq. 22). To get the value of 
, we need to evaluate the time-dependent coefficients in the solution 
, 
, and 
. This trial function is assumed to be the solution for 
as the functions 
are the eigenfunctions satisfying the given problem as given in Ref. 5. In general, a polynomial form is assumed.7 However, having the solution for constant boundary conditions as the trial function helps in simplifying the integrals. The choice is up to the researchers; however, eigenfunctions for problems with constant boundary conditions and linear models are a good basis for nonlinear and time-varying boundary conditions. This idea follows the Duhamel's superposition method1 wherein models with time-varying boundary conditions are obtained from constant boundary condition models using the convolution theorem.
The coefficients are obtained by solving the final equation using the boundary conditions given above and solving for 
. We introduce the average concentration as
The coefficients 
and 
are obtained in terms of 
and simplifying the expressions substituting 
we have
Substitution in 
yields
Substitution of 
in our original Eq. 6 yields
Calculating 
from the above equation, multiplying with differential volume 
, also multiplying by 
and integrating to find 
, after simplification we have
To solve the above equation efficiently, we introduce variable 
such that
Substituting the above relation and integrating the above equation we have
where 
is a constant for integration.
Grouping like terms we may write Eq. 18 as
Simplifying and substituting 
, the above equation can be written as
The above equation can be written as
Substituting 
does not alter the value of 
as 
remains unchanged. Hence, the final set of equations is not affected by the value of integration constant 
.
The final set of equations can be written as
An important advantage of this approach is that this reformulation is very robust, and by having enough 
values, in other words, enough number of Eq. 24, this approach can cover the entire spectrum of high/low rates, pulses, etc., such as the PSS method by Liu. However, in Liu's PSS approach, the 
values vary as 
and the order of 
is as high as 
causing stiffness and numerical instability in the pseudo-2D models using this approach. The 
values mentioned here are the 
values in the final set of equations obtained from the PSS method5 as mentioned in Eq. 9a-c in Ref. 2. The present reformulation overcomes this problem. In this method, the 
values vary as 
; in other words, we have a converging series in 
, which makes this approach equivalent to the PSS model in accuracy but highly efficient in a pseudo-2D environment for computation avoiding stiffness and computational difficulties.
Finite Difference Approach with Unequal Node Spacing
Finite Difference method is one of the most widely used numerical techniques to solve ordinary or partial differential equations. The use of the FD method has been the first choice for solving first principles based lithium-ion battery models. However, for a pseudo-2D model, when dealing with a second dimension 
for discretization, the number of equations increases by many folds and thereby making the simulation of the system slower and complex. About 20 node points in the 
direction increase the number of equations by a great deal, and hence, based on the mixed order finite element approach discussed earlier, where the size of linear elements were unequal instead of fixed equal elements, we used a mixed order FD approach, wherein we use less node points with unequal node spacing. The derivation of FD notations for different approximations for the derivatives is given in the following section.
Taylor series expansions at 
and 
are written as
where 
is the unequal node spacing between the 
node in the domain.
Truncating the series expansion with the required amount of accuracy and solving for the first and second-order derivatives, we can obtain the formula for the central FD for the first and second-order derivatives. We use an order of 
accuracy for all our approximations
Similarly, forward and backward FD relations for the derivatives can be obtained and used for boundary conditions.
Figure 1 presents a general methodology for obtaining an efficient reformulation/representation for the solid-phase equation for nonlinear diffusivities. First, a mixed FD representation is written, say with 
node points. For the optimization, 
is the initial guess with 
as the constraint and the error between the expected full-order numerical solution and the mixed FD method is minimized to a set tolerance. Typically, Jacobian-based methods are sufficient for convergence.8 For difficult nonlinearities, global and robust optimization might be needed for convergence and robustness9, 10 though they are likely to be very slow compared to direct Jacobian based methods.
Figure 1. Schematic of steps involved in mixed FD method for optimized spacing and hence reformulation of solid-phase diffusion.
For the solid-phase diffusivity 
varying as a function of concentration, the Galerkin approach cannot be used. We used FD with unequal node spacing in the 
direction and discretized the diffusion equation.6 The mixed FD form of this equation with constant 
using the above derived FD stencil is given as
where 
is the number of interior node points. A similar expression for varying 
can be derived.
One of the advantages of this method is that, for our case, the concentration gradient is nearer the surface compared to the center, and hence, strategically placing more node points near the surface and less node points at the center would give results with the same accuracy with lesser total number of node points compared to a large number of equally spaced node points. Lesser number of node points in 
leads to lesser number of equations and hence faster simulation for the whole battery model. The placement of these node points are important, and to find the exact position of these node points, we ran an optimization algorithm to find the best 
, 
, 
, etc., and minimized 
and the central processing unit (CPU) time. This method is very accurate for short times/high rates/pulses and is recommended for varying diffusion coefficients. Varying diffusivities are important and is likely to get more attention because of its requirement for addressing stress effects in the Li-ion batteries.11 The approach is robust as an optimization algorithm proposed in this paper will automatically detect instead of having to guess and arrive at the node spacing by trial and error.
Coupling Solid-Phase Diffusion with Full-Order Pseudo-2D Battery Models
In the traditional formulation of solid-phase diffusion, Eq. 1, 2, 3, 4 are coupled to the equations of a full-order pseudo-2D model for lithium-ion batteries, which is described elsewhere.1, 12, 13 For comparison, two efficient methods, the eigenfunction based Galerkin method and the mixed FD method (with five internal nodes) are also coupled to the full-order pseudo-2D model.
Three pseudo-2D codes were written in Fortran and solved with the DASKR differential–algebraic equation solver, which is a root-finding version of DASPK.14 They consisted of the traditional FD in the 
and 
(with 50 node points in the 
direction and 35 in 
direction) model, the Galerkin model, and the mixed FD model for intraparticle diffusion. For all cases in Fig. 6, 
is constant.
Figure 6. Discharge curves at 5C and 10C rate for a pseudo-2D model for Li-ion battery: Comparison of full-order pseudo-2D, Galerkin based, and mixed FD methods for solid-phase diffusion.
Results and Discussion
Figure 2 shows the comparison of the Galerkin method, with traditional FD (full-order) numerical solution for varying 
values and a constant 
along with the PSS method. From this figure, results obtained with the full-order numerical solution (50 node points in 
) can clearly be efficiently obtained at the reduced computational time with no compromise in accuracy. Though this figure compares the results for a single PDE (solid-state diffusion alone), the results obtained help in simulating a pseudo-2D model with efficient approximation/representation for the solid phase.
Figure 2. Comparison of eigenfunction based Galerkin reformulation with full-order numerical solution and PSS by Liu for 
and 
.
Figure 3a shows the values of 
obtained when solving the above equation for the given sin function as input for the current. As described above, it can be seen that the values of 
decrease with increasing 
as well as with time as opposed to the PSS method thus giving a converging solution and also making it easier for stiff equation solvers to converge faster and more efficiently. The 
values from Ref. 2 are plotted in Fig. 3b, and it is clearly observed that they are a diverging series at low precision computations and reach very high values, which may cause the solvers to become unstable. Thus the Galerkin method provides a more efficient way of including the solid-phase diffusion without compromising accuracy for all possible operating conditions with constant diffusion coefficients.
Figure 4 shows the comparison of the mixed FD method with five internal nodes for a constant value of 
and a constant 
compared with the full-order numerical solution. The reformulation agrees accurately with the full-order numerical solution. Again, both at short times and long times, the mixed FD representation matches with the full-order solution. The values of optimized node spacing obtained in this case for different values of 
are 0.2183372643, 0.1779355824, 0.1228253438, 0.1698047152, 0.1499086011, and 0 0.1611884932.
Figure 4. Comparison of mixed FD method with five interior nodes with full-order numerical solution for constant 
and 
, etc.
Proton diffusion into nickel hydroxide electrodes used in Ni–MH batteries is a strong function of solid-phase concentration and decreases approximately 3 orders of magnitude. This varying transport property was captured by using the complex faradaic impedance of the nickel hydroxide active material and reported as Eq. 5 elsewhere.15 This work has been used for accounting the variable diffusion coefficient in Ref. 16 to determine a diffusion coefficient that is a function of the dimensionless flux rate of the material diffusing into the particle. Verbrugge and Koch17 expressed the intercalate diffusion coefficient as an indirect function of solid-phase concentration consisting of fraction occupancy of the intercalating host material and activity coefficient. The significance of taking an account of this variation in intercalating electrodes was demonstrated by Subramanian et al.18 Here, mathematical models are developed to simulate the potentiostatic charge/discharge of a partially graphitic carbon fiber and the galvanostatic discharge of a lithium foil cell under solid diffusion limitations. Evidence that shows the importance of accounting for nonlinear diffusion was shown by Karthikeyan et al.19 for the recently popular 
positive active material in lithium-ion batteries where the thermodynamic expressions along with the activity correction are incorporated into a single particle diffusion model for a Li-ion cell. Hence, the use of nonlinear diffusion, wherein the diffusion coefficient is a function of concentration, is becoming more and more popular in the battery modeling domain, and the mixed order FD method is capable of giving accurate results with nonlinear diffusivities as well. To illustrate this fact, Fig. 5 is presented with the comparison of mixed order FD method with five internal nodes for constant 
and 
varying as a simple function of 
with the full-order numerical solution. For varying 
, 
, mixed FD approach was efficient and accurate at short times.
Figure 5. Comparison of mixed FD method with five interior nodes with full-order numerical solution for 
and 
.
Comparisons between the solid-phase diffusion models for the full-order pseudo-2D model and approximations developed show excellent agreement. The time history of the cell voltage was monitored for several discharge rates. For very high discharge rates of 5C and 10C, Fig. 6 shows the comparison of the full-order Galerkin (with five 
values) and mixed FD (with five internal nodes) methods for constant 
. The computations were terminated when the potential dropped to 2.5 V. Agreement is very good between the traditional FD and the mixed FD methods, indicating that this efficient method could substitute for the traditional method for any discharge rate (low or high). The mixed FD method agrees remarkably with the full-order FD method and can also be used for nonlinear diffusivities and hence can increase the computational efficiency of the whole battery model. The eigenfunction based Galerkin collocation approach included in a full pseudo-2D Li-ion battery model also shows excellent agreement at high rates of discharge indicating that this is a good alternative for cases with constant diffusivities.
Table I shows the CPU time taken for performing the above mentioned simulations and it can clearly be seen that the two proposed methods take lesser simulation time compared to the simulation time for a full order pseudo-2D model. The eigenfunction based Galerkin method takes more time compared to the mixed FD method. The solid-phase reformulations are necessary for a faster simulation of Li-ion battery models, which help in faster estimation of parameters from these models from the experimental data.20
Table I. Comparison of CPU times taken for full order pseudo-2D Galerkin-based and mixed FD methods for obtaining discharge curves in Fig. 6 at 5C and 10C rates.
| Method | CPU time for 5C rate(s) | CPU time for 10C rate(s) |
|---|---|---|
| Full order pseudo-2D | 19.2 | 20 |
| Galerkin based | 4.5 | 4.84 |
| Mixed FD | 1.438 | 1.297 |
Conclusion
The different approximation schemes available for solid-phase approximation for Li-ion batteries were reviewed and certain disadvantages pertaining to some of those methods were discussed. To overcome these small disadvantages, an eigenfunction based Galerkin weighted residual approximation was presented, which provides efficient reformulation for the solid-phase equation for constant diffusivities. The result from this method compares very well with the PSS method and also facilitates easily converging solutions for the 
values. Mixed order FD based on finite volume equations can be derived for varying values of diffusion coefficients as a function of concentration (nonlinear case) and are very efficient for short times and so far seem to be the only option for reformulating nonlinear diffusivities efficiently. The two methods presented here seem to be better compared to the existing solid-phase approximations with absolutely no shortcomings as they are valid for any kind of operating conditions. For the linear diffusivity (
is constant), the eigenfunction based Galerkin collocation is the most efficient method, which does not compromise on accuracy, however, providing excellent computational efficiency. For nonlinear diffusivities, the mixed order FD method can be used to obtain accurate solutions by using lesser optimally spaced node points thereby reducing the total number of equations being solved for the full pseudo-2D Li-ion battery model.
Future work involves efficient reformulation for moving boundary models for phase change materials and coupling the solid-phase models with reformulated model12 to enable simulation at very high rates.
Acknowledgments
The authors are thankful for the partial financial support of this work by the NSF: SGER 0609915 , CBET 0828002 , and CBET 1008692 and the United States government.
Washington University in Saint Louis assisted in meeting the publication costs of this article.
List of Symbols
| reference concentration, 
|
| concentration of lithium ions in the intercalation particle of electrode, 
|
| average concentration of lithium ions in the intercalation particle of electrode, 
|
| dimensionless concentration of lithium ions in the intercalation particle of electrode |
| dimensionless average concentration of lithium ions |
| dimensionless concentration at  node point |
| Li-ion diffusion coefficient in the intercalation particle of electrode, 
|
| diffusion coefficient at reference concentration  , 
|
| node spacing at  node point |
| pore wall flux of Li ion the intercalation particle of electrode, 
|
| volume averaged concentration flux, 
|
| radius of the intercalation particle of electrode, m |
Greek
| positive eigenvalues |