An Efficient and Robust Numerical Solution of the Full-Order Multiscale Model of Lithium-Ion Battery

We propose a novel and efficient numerical approach for solving the pseudo two-dimensional multiscale model of the Li-ion cell dynamics based on first principles, describing the ion diffusion through the electrolyte and the porous electrodes, electric potential distribution, and Butler-Volmer kinetics. The numerical solution is obtained by the finite difference discretization of the diffusion equations combinedwith an original iterative scheme for solving the integral formulation of the laws of electrochemical interactions. We demonstrate that our implementation is fast and stable over the expected lifetime of the cell. In contrast to some simplified models, it provides physically consistent results for a wide range of applied currents including high loads. The algorithm forms a solid basis for simulations of cells and battery packs in hybrid electric vehicles, with possible straightforward extensions by aging and heat effects.


Introduction
Modern Li-ion batteries possess advantages making them a popular choice in many different applications.Their light weight, low self-discharge rate, and performance especially matter for power storage in hybrid electric vehicles (HEV).There are ongoing efforts to optimize the control strategy of the HEV powertrain in order to improve not only the range of the vehicle, but also the total useful battery capacity during its lifetime.Mathematical models of Li-ion cells based on first principles provide insight into the dynamics of the battery cycling and the information from the computational simulations can be used during the design of the control algorithms.
The basic reference for simulating Li-ion cell dynamics is the isothermal model proposed by Newman and Tiedemann [1] and Doyle et al. [2].For the model with thermal effects included, we refer the reader, for example, to Cai and White [3], Kumaresan et al. [4], or Gu and Wang [5].The aging effects of Li-ion batteries are discussed, for example, by Ramadass et al. [6,7] or Ning et al. [8,9].
Recently, several approximate techniques were successfully applied to reduce the computational complexity of these models.In [10,11], Subramanian et al. used perturbation techniques.Model reduction and Chebyshev polynomial methods were used by Bhikkaji and Söderström in [12].In Smith et al. [13], the residue grouping method was used.Cai and White [14] developed a reduced-order model by means of orthogonal decomposition.
In our work, we strive to create a robust yet computationally efficient algorithm for predicting the state of Liion batteries subject to intense and variable loading over extended time periods.We adopt the full-order pseudo twodimensional model of Li-ion cell dynamics describing ion diffusion through the electrolyte, charge flow, and the Butler-Volmer kinetics, as summarized in [15].We propose a novel approach to the solution of electrochemical interactions by means of integral reformulation of the governing equations and an iterative scheme for their solution.As a result, our algorithm remains stable for a wide range of applied currents.The simulations are fast enough to cover the 2 Mathematical Problems in Engineering long term behavior of the battery while still resolving both the macroscale diffusion processes across the cell as well as position-dependent microscale dynamics in the porous material of the electrodes (in contrast to the single-particle models; see, e.g., [16,17]).
The paper is structured as follows.In Section 2, we present the summary of the mathematical model.Section 3 is dedicated to the derivation of the integral solution of electrochemical interactions, leading to the formulation of an iterative algorithm.In Section 4, we shortly comment on the implementation of the algorithm so as to be able to explain some further ideas.In Section 5, we first use our model to replicate the study by Subramanian et al. [18] and Dao et al. [15] and discuss the obtained outcomes.Further on, we proceed with the analysis of the effect of the numerical algorithm parameters.We draw conclusions about the applicability of the proposed approach for long term simulations.

Summary of the Mathematical Model
2.1.Geometrical Setting.As shown in Figure 1, the onedimensional representation of the lithium-ion cell computational domain Ω = [  ,   ] is divided into three parts such that Ω = Ω 1 ∪ Ω 2 ∪ Ω 3 , where Ω 1 = [  ,   ] is the positive electrode, Ω 2 = [  ,   ] is the separator, and Ω 3 = [  ,   ] is the negative electrode.Any single point  0 ∈ Ω corresponds to a cross section through the real three-dimensional cell by the plane  =  0 .Neither the area nor the shape of this cross section are known or needed in the 1D model.All respective quantities such as the applied current density  app and the mass fluxes are calculated per unit cross section area.
In the following text, the quantities corresponding to the solid electrodes and electrolyte are indexed by  and , respectively, and the quantities defined in Ω  are enumerated by the appropriate subdomain index  ∈ {1, 2, 3}.The values of all quantities depend on time  which lies within the interval J = (0,  final ).

Diffusion of Li
+ in the Electrolyte.Depending on the mode of operation (charge/discharge), the lithium ions are extracted (deintercalated) from the porous material of one electrode, transferred through the electrolyte across the separator by diffusion, and finally intercalated into the porous material of the other electrode.Based on [15], the governing equations for diffusion read for each  ∈ {1, 2, 3}.The meanings of the symbols in (1) are as follows.
, [mol m −3 ] denotes the concentration of Li + in the electrolyte.
[m 2 s −1 ] is the diffusion coefficient of Li + in the electrolyte.
[1] is the tortuosity factor of the porous medium given by the Bruggeman relationship [19] , where   [1] is the porosity of the medium and brugg  [1] is the Bruggeman coefficient. 0 + [1] is the transference number of Li + in the electrolyte.
[mol m −2 s −1 ] denotes the (de)intercalation flux of Li + from the internal surface of the porous material into the electrolyte per unit surface area ( 2 = 0 as no lithium ions are stored in the material of the separator).  [m −1 ] is the internal surface area of the porous material per unit volume.
Equation ( 1) is accompanied by a number of boundary conditions.First, the lithium ions cannot leave the outer boundary of the cell which implies The other boundary conditions ensure continuity of the concentration in the electrolyte at the interdomain boundaries The initial conditions are given by

Diffusion inside the Porous Electrodes.
The electrodes are made of porous material, that is, a mixture of void space (filled with electrolyte) and a solid continuum.At the microscopic level, the solid matrix is modeled in the form of small spherical particles releasing (deintercalating) or absorbing (intercalating) lithium ions through their surface.Lithium then diffuses through each particle in the radial direction according to the current concentration distribution.The material balance for lithium in a single active solid material particle in the positive or negative electrode is governed by Fick's second law in spherical coordinates [3,15]: for each  ∈ {1,3}.The distance from the center of the spherical particle  lies in the interval Ω part, = (0,  , ).
, is the radius of the solid particles in the electrode Ω  .
, [mol m −3 ] is the concentration of lithium in the particle.
[m 2 s −1 ] is the diffusion coefficient of lithium in the particle.
The initial concentration distribution in a particle is given by At the center of the particle, the boundary condition The solid matrix occupies a volume fraction 1 −   −  , , because apart from electrolyte, there can be some amount of inert material (filler) with the volume fraction  , [15].According to (10), the number of spherical particles per unit volume at  ∈ Ω  is All such particles are assumed to have the same radial distribution of concentration  , .In particular, we denote the concentration at the particle surface as The total surface area of these particles per unit volume is a material property.Its values for some porous materials can be found in literature [20], allowing the calculation of the radius of the respective spherical particles  , from (13).

Electrochemical Interactions.
Equations ( 1) and ( 6) determine the chemical state of the cell provided that the fluxes  1 and  3 are given.The connection between the concentrations of lithium ions and the respective mass fluxes comes from the modeling of the electrochemical interactions in the cell which also allows the calculations of other quantities such as the cell voltage.Whereas the diffusion equations in Sections 2.2 and 2.3 are readily prepared for numerical solution by standard tools (e.g., the method of finite differences), the equations summarized below require nontrivial treatment.By the end of Section 2.4, the equations do not contain partial derivatives with respect to  and .For any function  : Ω  × J → R,  ∈ {1, 2, 3}, we therefore simplify the notation / to   .

Electrical Potential in the Porous
Electrodes.The charge continuity equations in the solid electrodes placed at Ω 1 and Ω 3 are given by Ohm's law [15] as with the boundary conditions where is the solid phase electrical potential, is the current density applied to the electrode ( app > 0 corresponds to charging and  app < 0 to discharging), denotes the charge flux in terms of Li + ions from the internal surface of the porous material into the electrolyte per unit surface area, is the effective electronic conductivity defined as [S m −1 ] is the electronic conductivity.

Electrical Potential in the Electrolyte.
Based on [15], the charge continuity equation in the electrolyte is given by with the boundary conditions where is the electrical potential in the electrolyte, is the effective ionic conductivity of the electrolyte given by  eff  =     that accounts for the tortuous path of the porous medium, where   [S m −1 ].
For the reference ionic conductivity   , we use the empirical correlation reported in [15]: with The meanings of the symbols in (18a) and (18b) are the following: is the reaction rate coefficient,  , and  , [1] are the anodic and cathodic transfer coefficients of electrochemical reaction,  ,max, [mol m −3 ] is the saturated concentration of Li + ions in the solid phase,   [V] is the intercalation overpotential described as In (19),   [V] is the open circuit potential determined by the following empirical correlations: where   =  ,surf, / ,max, ,  = 1, 3.

Integral Solution of Electrochemical Interactions
For a given time , the flux   ,  ∈ {1, 3}, can be calculated by the solution of the system of two differential equations (14a) and (16a) and the algebraic equation (18a) for the unknowns  , ,  , ,   , and   , together with the respective boundary conditions (see Section 2.4).It turns out that this set of equations can be transformed into a system of two ordinary differential equations (ODEs) for   and   , where by   we denote the integrated mass transfer function defined by and  ℓ, ,  , represent the left and right boundary coordinates of Ω  , respectively.

Derivation of the Differential Equations.
In order to derive the ODE for   , we first integrate (14a) with respect to  over the interval ( ℓ, , ) ⊂ Ω  ,  ∈ {1, 3}.Using the boundary conditions (14b) and (14d), we obtain Similarly, by integrating (16a) for  ∈ {1,2,3} and using the boundary conditions (16b), (16c), and (16d), we arrive at Subtractions of (22c) from (22a) and (22e) from (22b) allow expressing the derivative of   defined by (19) as The ODE for   is given directly by the Butler-Volmer reaction kinetics (18a) in the form Differentiating (24) with respect to  and plugging in   from (23a) and (23b), respectively, lead to a single second-order ODE for   .The corresponding boundary conditions follow from the definition of   given by ( 21) and evaluation of (22a), (22b) at  =   and  =   together with the boundary conditions given by (14c) and (14e).They read and they complete two well posed problems for the unknowns  1 and  3 .However, in the following, we use (23a) and (23b), (24), and (25a) and (25b) and transform them into integral equations that can be solved iteratively by means of numerical integration.

ODE System in General
Form.Assume a general coefficient form of ODEs ((23a) and (23b)) and (24) in the following compact form: for all  ∈ (  ,   ) with the following boundary conditions: where   =  , ,   =  , ,   =  , ,   =  , , the coefficients  =   are given by (18b), and  =   ,  =   , and  =   correspond to the following domain-defined coefficients (denoted by the indices  = 1 and  = 3): where, in (28a),  = (  ) is the unknown integration constant and in (28b) the boundary condition given by (26d) was already employed.Equations (28a) and (28b) can be combined into a single integral equation in two different ways: Under the assumption of constant temperature ,  becomes constant in (26c).As a result, (26c) and (31) can be combined to obtain where Λ = exp() denotes the term with the unknown parameter  and the coefficients   and   read as If   =   , (32) is easily resolved by In general, (32) is a highly nonlinear equation that needs to be solved numerically.

Iterative Scheme for Solving the Integral Equation.
We propose the following iterative scheme for solving the integral equation given by (30): where   is the solution of (32) which has to be updated at every iteration step for given   ,  = 0, 1, 2, . .., and  ∈ (0, 1] is a relaxation coefficient that serves as a tuning parameter allowing control of the convergence of the iteration scheme.Note that, with values of  close to 1, the iterative scheme given by (35) diverges rapidly because of the exponential functions in (26c).As the initial guess in (35), we choose  0 ≡ 0 in the first call to the iterative solver and the previously calculated value of  in the subsequent calls.The iterative process is terminated when the norm of J+1 is below a given threshold .The values of  and  are discussed further in Section 5.2.

Computational Algorithm for Solving the System of ODEs.
Let us summarize the algorithm that allows computing the electrolyte and the solid phase electrical potentials and the respective charge fluxes.The steps of the computational algorithm are as follows: (1) In the given time step, use the prescribed profiles of  , and  ,surf, and the value of  app to solve (26a), (26b), (26c), (26d), and (26e) using the iteration scheme given by ( 35) to obtain the integrated charge fluxes  1 and  3 .
(8) Compute  ,1 by integrating (22a) and using the boundary condition that follows from (19) as (9) The value of  ,1 (  ) represents the external apparent voltage of the cell.

Implementation of the Numerical Algorithm
There are two diffusion processes on different time and spatial scales in the model.They are both solved by the implicit Euler scheme of the finite difference method [22] with generally different time steps.In Ω, (1) is solved on a grid of  uniformly spaced nodes  0 , . . .,  −1 and the positions of the individual grid nodes determine their correspondence to the domains Ω  ,  ∈ {1, 2, 3}.This grid is also used for numerical integration in the algorithm described in Section 3.For each node   ∈ Ω  ,  ∈ {1, 3}, another mesh of  nodes exists that discretizes the interval [0,  , ], where ( 6) is solved.The algorithm uses multiple time scales to update the individual quantities.It carries out the diffusion in the solid particles for each   with a constant time step Δ.After a given number of time steps, the integral solver is called, which updates the values of   for each   .For the time period between the updates of   , the macroscopic diffusion in the electrolyte governed by ( 1) and the microscopic diffusion in the spherical particles governed by ( 6) are completely independent.This allows a suitable multiple of Δ to be used as the time step for the diffusion in the electrolyte.

Simulation Results
In this section, we demonstrate the capabilities of our algorithm and investigate its behavior depending on the settings of the numerical solver parameters.For all simulations, we use the model parameters taken from [15,18], as summarized in Table 1.

Single Discharge Cycle.
We compare the results of our algorithm to those described in [15,18].In Figure 2(a), the evolution of cell voltage during one discharge cycle with low discharge currents (0.5C, 1C) is shown.Our results coincide almost completely with the simulations using the full-order finite difference model by Subramanian et al. [18].There is also a fair agreement with the results of the simplified and reduced model proposed by Dao et al. [15].The notable exception is the initial part of the evolution where Dao et al. observe no voltage drop because of loading.The situation changes for higher discharge currents (2C, 3C, and 4C), as plotted in Figure 2(b).In this case, the results of the simplified model are completely different from ours.Again, the simplified model does not account for immediate voltage drop due to an increased load.On the other hand, the voltage readings at  = 0 based on our model form an almost perfectly linear dependence on the applied current, as can be seen in Figure 3.This is in correspondence with the usual representation of battery cells in DC electrical circuits where a constant internal resistance is considered.
The above arguments indicate that the reduced model ceases to be valid for higher currents while our proposed model continues to provide expected and consistent results (see also the model comparison in [17]).Unfortunately, the results of the full-order model from [18] are not available for currents above 1C, although the authors state that they are able to simulate such situations.The insufficient information contained in the reduced model can also be demonstrated on Li + concentration profiles in the electrolyte.Significant differences occur even for low currents, as can be seen in Figure 4.In Figure 5, the concentration profile comparison with another implementation of the full-order model [3] is provided.Qualitatively, the agreement is satisfactory.However, complete match of the curves could not be achieved due to the lack of information about the model parameters in [3].

Properties of the Numerical Solver.
We are interested in the influence of the numerical solver parameters on the accuracy and speed of the simulation.Several tests were performed, involving the parameters explained in Section 4. As the test vehicle, the voltage curve of a single discharge cycle with various values of the applied current was used.
(i) We tested several grid resolutions in the spatial domain Ω, ranging from  = 20 to  = 800.
For each value of , uniform grids from  = 20 to  = 400 were utilized for diffusion inside the spherical particles.For the 1C applied current and  fixed, the solution is almost independent of the value  of  provided that  ≥ 50.For different values of , there are slight differences mostly in the final part of the voltage curve, as demonstrated in Figure 6.As  increases, the differences between the subsequent cases become smaller, which suggests convergence of the numerical method.No rigorous convergence tests (such as measuring the experimental order of convergence [23]) were performed, though.With higher applied currents, the value of  also affects the shape of the voltage curve.For the 2C applied current, the differences in the results become negligible for any combination of  and  satisfying  > 200,  > 200.
The simulation of the whole 1C discharge cycle with  =  = 100, Δ = 10 −1 s,  = 0.05, and  = 10 −13 can be performed in less than 10 seconds on a single core of an Intel i7-6700K @ 4 GHz CPU, which is 350x faster than the real discharge process on average.Our computational tests indicate that decreasing  and  below approximately 50 nodes (and losing accuracy as seen in Figure 6) is not necessarily beneficial for the computational time.(ii) As explained in Section 4, the integral solver need not be called in every time step.We performed some computations with a fixed (and unnecessarily small) time step Δ = 10 −3 s.The integral solver was set to be called in the intervals ⋅Δ, for several different values of  ∈ N. In addition, the diffusion in electrolyte was also solved using the time step  ⋅ Δ.All the results were virtually identical.The numbers of iterations of the integral solver for all these cases are summarized in Figure 7.As expected, the number of iterations in both Ω 1 and Ω 3 is higher when the solver is called less frequently because the previous solution is used as the initial guess for the next iteration.In fact, calling the solver in long time steps brings little computational time savings, as demonstrated in Figure 8.
(iii) Next, the influence of the relaxation parameter  was investigated.The maximum value of  sufficient for the convergence of the integral solver depends on the applied current and other settings (e.g., the approximate bounds are  < 10 −1 for 1C discharge and  < 2 × 10 −3 for 4C discharge).However, once such value is found, it is undesirable to further decrease  as it only leads to prolonged computational times.Provided that the integral solver converges, its accuracy is only controlled by the value of  and is independent of the setting of .
For long term simulations with other included effects such as battery aging, the basic version of the algorithm has to be stable over extended time periods.We performed a stability test by simulating over 1000 constant current charge/discharge cycles in the prescribed voltage range.The results in Figure 9 testify that the algorithm not only exhibits excellent stability, but also allows long simulations as its computational time is many hundreds of times shorter than the real duration of the simulated processes.

Conclusion
We have developed an efficient numerical algorithm for the solution of the full-order version of the well known model of Li-ion cell dynamics, as used by [15,18].For low to moderate applied currents, the obtained simulation results are in good agreement with the studies performed in both [15,18].For higher discharge rates, our algorithm proves to maintain physically consistent behavior in contrast to the simplified model by [15].Moreover, its implementation is fast and stable enough to enable cycling simulations over the expected lifetime of the cell.These properties justify the choice of the full-order model for our ongoing efforts to simulate the behavior of both individual cells and battery packs installed in hybrid electric vehicles.The proposed numerical algorithm is a convenient basis for such efforts, as it allows straightforward generalizations in order to incorporate heat effects and aging     phenomena in the scope of the models found, for example, in [3,[6][7][8][9].

Figure 1 :
Figure 1: One-dimensional representation of the lithium-ion cell.

Figure 2 :
Figure 2: Cell voltage evolution during one deep discharge cycle.(a) Discharge currents 0.5C and 1C; comparison of our results (thick solid lines in the background) with Dao et al. ([15], dashed lines) and Subramanian et al. ([18], dotted lines).(b) Discharge currents 2C, 3C, and 4C; comparison of our results (thick lines) with Dao et al. ([15], thin lines).The same line styles correspond to the same discharge current.

Figure 3 :
Figure 3: The almost perfectly linear dependence of the cell voltage on the applied current.The solid line connects the subsequent data points; the dashed line connects the first and the last point.

Figure 4 :Figure 5 :
Figure 4: Comparison of the spatial profiles of   at selected time levels.Thick lines represent our results, thin lines are the results from [15].The same line styles correspond to the same time .Discharge current 1C.

Figure 6 :
Figure 6:  The voltage curves for one 1C discharge cycle in simulations using different mesh resolutions for the discretization of the spherical particles.The marked region of the plot is zoomed in to better illustrate the differences. = 100, Δ = 10 −1 s,  = 0.05, and  = 10 −13 .Integral solver called in every time step.

Figure 9 :
Figure9: Algorithm stability test during 1C constant current discharge/charge cycling within the voltage range from 3.4 V to 4.5 V.A total of 1043 cycles in the physical time period of 4.41×10 6 s (51 days) have been simulated.The voltage curve is plotted by the solid line for the first four cycles and by the thick dashed line for the last four cycles.In addition, the cycle numbers and cycle durations are indicated.Numerical parameters:  = 100,  = 50, Δ = 10 −1 s,  = 0.05, and  = 10 −14 .The integral solver was called in every time step.The computation took less than 3 hours on a single core of an Intel i7-6700K @ 4 GHz CPU.