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

Authors: Michal Beneš; Radek Fučík; Vladimír Havlena; Vladimír Klement; Miroslav Kolář; Ondřej Polívka; Jakub Solovský; Pavel Strachota
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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 combined with 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.


---

## Body

## 1. 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 Li-ion batteries subject to intense and variable loading over extended time periods. We adopt the full-order pseudo two-dimensional 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 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 Section2, 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.


## 2. Summary of the Mathematical Model

### 2.1. Geometrical Setting

As shown in Figure1, the one-dimensional representation of the lithium-ion cell computational domain Ω=[xa,xd] is divided into three parts such that Ω=Ω1∪Ω2∪Ω3, where Ω1=[xa,xb] is the positive electrode, Ω2=[xb,xc] is the separator, and Ω3=[xc,xd] is the negative electrode. Any single point x0∈Ω corresponds to a cross section through the real three-dimensional cell by the plane x=x0. 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 Iapp and the mass fluxes are calculated per unit cross section area.Figure 1
One-dimensional representation of the lithium-ion cell.In the following text, the quantities corresponding to the solid electrodes and electrolyte are indexed bys and e, respectively, and the quantities defined in Ωk are enumerated by the appropriate subdomain index k∈1,2,3. The values of all quantities depend on time t which lies within the interval J=0,tfinal.


### 2.2. 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(1)ϵk∂ce,k∂t=∂∂xDeτk∂ce,k∂x+1-t+0akjkinΩk×Jfor each k∈1,2,3. The meanings of the symbols in (1) are as follows. ce,kmolm-3 denotes the concentration of Li+ in the electrolyte.ϵk[1] is the material porosity (void fraction).De[m2s-1] is the diffusion coefficient of Li+ in the electrolyte.τk[1] is the tortuosity factor of the porous medium given by the Bruggeman relationship [19] τk=ϵkbruggk, where ϵk[1] is the porosity of the medium and bruggk[1] is the Bruggeman coefficient.t+0[1] is the transference number of Li+ in the electrolyte.jk[molm-2s-1] denotes the (de)intercalation flux of Li+ from the internal surface of the porous material into the electrolyte per unit surface area (j2=0 as no lithium ions are stored in the material of the separator).ak[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 (2)∂ce,1∂xx=xa=∂ce,3∂xx=xd=0.The other boundary conditions ensure continuity of the concentration in the electrolyte at the interdomain boundaries (3)ce,2x=xb=ce,1x=xb,ce,3x=xc=ce,2x=xc,as well as the continuity of the interdomain concentration fluxes (4)Deτ2∂ce,2∂xx=xb+=Deτ1∂ce,1∂xx=xb-,Deτ3∂ce,3∂xx=xc+=Deτ2∂ce,2∂xx=xc-.The initial conditions are given by (5)ce,kt=0=ce,k0,inΩk,∀k∈1,2,3.


### 2.3. 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]:(6)∂cs,k∂t=Dsr2∂∂rr2∂cs,k∂rinΩpart,k×Ωk×Jfor each k∈1,3. The distance from the center of the spherical particle r lies in the interval Ωpart,k=0,Rs,k. Rs,k is the radius of the solid particles in the electrode Ωk.cs,k[molm-3] is the concentration of lithium in the particle.Ds[m2s-1] is the diffusion coefficient of lithium in the particle.The initial concentration distribution in a particle is given by(7)cs,kt=0=cs,k0inΩpart,k×Ωk.At the center of the particle, the boundary condition(8)-Ds∂cs,k∂rr=0=0imposes zero flux of lithium. On the particle surface, the flux is equal to the consumption/production rate of Li+ due to the electrochemical reaction occurring at the solid/liquid interface; that is,(9)-Ds∂cs,k∂rr=Rs,k=jk.


#### 2.3.1. Scale Coupling

For eachx∈Ωk, the number of spherical particles per unit cross section area and in the range x,x+dx is(10)dnk=1-ϵk-ϵf,k4/3πRs,k3dx.The solid matrix occupies a volume fraction 1-ϵk-ϵf,k, because apart from electrolyte, there can be some amount of inert material(filler) with the volume fraction ϵf,k [15]. According to (10), the number of spherical particles per unit volume at x∈Ωk is (11)nk=1-ϵk-ϵf,k4/3πRs,k3.All such particles are assumed to have the same radial distribution of concentration cs,k. In particular, we denote the concentration at the particle surface as(12)cs,surf,kx,t=cs,kRs,k,x,t.The total surface area of these particles per unit volume(13)ak=nk4πRs,k2=31-ϵk-ϵf,kRs,kis 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 Rs,k from (13).


### 2.4. Electrochemical Interactions

Equations (1) and (6) determine the chemical state of the cell provided that the fluxes j1 and j3 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 Section2.4, the equations do not contain partial derivatives with respect to t and r. For any function f:Ωk×J→R, k∈{1,2,3}, we therefore simplify the notation ∂f/∂x to f′.


#### 2.4.1. 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 (14a)σkeffϕs,k′′=akFjkinΩk×J,k∈1,3with the boundary conditions(14b)-σ1effϕs,1′xa=Iapp,(14c)-σ1effϕs,1′xb=0,(14d)-σ3effϕs,3′xc=0,(14e)-σ3effϕs,3′xd=Iapp,(14f)ϕs,3xd=0,where ϕs,k[V] is the solid phase electrical potential,F[Cmol-1] is Faraday’s constant,Iapp[Am-2] is the current density applied to the electrode (Iapp>0 corresponds to charging and Iapp<0 to discharging),Fjk[Am-2] denotes the charge flux in terms of Li+ ions from the internal surface of the porous material into the electrolyte per unit surface area,σkeff[Sm-1] is the effective electronic conductivity defined as (15)σkeff=1-ϵk-ϵf,kσk,σk[Sm-1] is the electronic conductivity.


#### 2.4.2. Electrical Potential in the Electrolyte

Based on [15], the charge continuity equation in the electrolyte is given by (16a)κkeffϕe,k′′=-akFjk+2RTF1-t+0κkeffce,k′ce,k′inΩk×J,k∈1,2,3with the boundary conditions(16b)-κ1effϕe,1′xa=0,(16c)-κ1effϕe,1′xb=-κ2effϕe,2′xb=Iapp,(16d)-κ2effϕe,2′xc=-κ3effϕe,3′xc=Iapp,(16e)-κ3effϕe,3′xd=0,(16f)ϕe,1xb=ϕe,2xb,(16g)ϕe,2xc=ϕe,3xc,where ϕe,k[V] is the electrical potential in the electrolyte,R[Jmol-1K-1] is the universal gas constant,T[K] is the temperature,κkeff[Sm-1] is the effective ionic conductivity of the electrolyte given by κkeff=τkκk that accounts for the tortuous path of the porous medium, where κk[Sm-1].For the reference ionic conductivityκk, we use the empirical correlation reported in [15]:(17)κk=4.1253·10-2+5.007·10-4ce,k-4.7212·10-7ce,k2+1.5094·10-10ce,k3-1.6018·10-14ce,k4,k∈1,2,3.


#### 2.4.3. Butler-Volmer Reaction Kinetics

Lithium is conserved and its fluxes from the solid particles and into the electrolyte are both equal tojk. The charge transfer Fjk is proportional to the mass transfer and is subject to Butler-Volmer reaction kinetics [6, 7, 15, 21] in the form (18a)jk=δkexp⁡αa,kFRTηk-exp⁡-αc,kFRTηk,k=1,3,with(18b)δk=Kkcs,max,k-cs,surf,k1/2cs,surf,k1/2ce,k1/2.The meanings of the symbols in (18a) and (18b) are the following: Kk[mol-1/2m5/2s-1] is the reaction rate coefficient,αa,k and αc,k[1] are the anodic and cathodic transfer coefficients of electrochemical reaction,cs,max,k[molm-3] is the saturated concentration of Li+ ions in the solid phase,ηk[V] is the intercalation overpotential described as(19)ηk=ϕs,k-ϕe,k-Uk.In (19), Uk[V] is the open circuit potential determined by the following empirical correlations: (20)U1=-4.875+5.839θ1-1.507θ13+0.531θ15θ1-1.005,U3=0.15-0.10θ3+0.00778θ3,where θk=cs,surf,k/cs,max,k, k=1,3.


## 2.1. Geometrical Setting

As shown in Figure1, the one-dimensional representation of the lithium-ion cell computational domain Ω=[xa,xd] is divided into three parts such that Ω=Ω1∪Ω2∪Ω3, where Ω1=[xa,xb] is the positive electrode, Ω2=[xb,xc] is the separator, and Ω3=[xc,xd] is the negative electrode. Any single point x0∈Ω corresponds to a cross section through the real three-dimensional cell by the plane x=x0. 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 Iapp and the mass fluxes are calculated per unit cross section area.Figure 1
One-dimensional representation of the lithium-ion cell.In the following text, the quantities corresponding to the solid electrodes and electrolyte are indexed bys and e, respectively, and the quantities defined in Ωk are enumerated by the appropriate subdomain index k∈1,2,3. The values of all quantities depend on time t which lies within the interval J=0,tfinal.


## 2.2. 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(1)ϵk∂ce,k∂t=∂∂xDeτk∂ce,k∂x+1-t+0akjkinΩk×Jfor each k∈1,2,3. The meanings of the symbols in (1) are as follows. ce,kmolm-3 denotes the concentration of Li+ in the electrolyte.ϵk[1] is the material porosity (void fraction).De[m2s-1] is the diffusion coefficient of Li+ in the electrolyte.τk[1] is the tortuosity factor of the porous medium given by the Bruggeman relationship [19] τk=ϵkbruggk, where ϵk[1] is the porosity of the medium and bruggk[1] is the Bruggeman coefficient.t+0[1] is the transference number of Li+ in the electrolyte.jk[molm-2s-1] denotes the (de)intercalation flux of Li+ from the internal surface of the porous material into the electrolyte per unit surface area (j2=0 as no lithium ions are stored in the material of the separator).ak[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 (2)∂ce,1∂xx=xa=∂ce,3∂xx=xd=0.The other boundary conditions ensure continuity of the concentration in the electrolyte at the interdomain boundaries (3)ce,2x=xb=ce,1x=xb,ce,3x=xc=ce,2x=xc,as well as the continuity of the interdomain concentration fluxes (4)Deτ2∂ce,2∂xx=xb+=Deτ1∂ce,1∂xx=xb-,Deτ3∂ce,3∂xx=xc+=Deτ2∂ce,2∂xx=xc-.The initial conditions are given by (5)ce,kt=0=ce,k0,inΩk,∀k∈1,2,3.


## 2.3. 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]:(6)∂cs,k∂t=Dsr2∂∂rr2∂cs,k∂rinΩpart,k×Ωk×Jfor each k∈1,3. The distance from the center of the spherical particle r lies in the interval Ωpart,k=0,Rs,k. Rs,k is the radius of the solid particles in the electrode Ωk.cs,k[molm-3] is the concentration of lithium in the particle.Ds[m2s-1] is the diffusion coefficient of lithium in the particle.The initial concentration distribution in a particle is given by(7)cs,kt=0=cs,k0inΩpart,k×Ωk.At the center of the particle, the boundary condition(8)-Ds∂cs,k∂rr=0=0imposes zero flux of lithium. On the particle surface, the flux is equal to the consumption/production rate of Li+ due to the electrochemical reaction occurring at the solid/liquid interface; that is,(9)-Ds∂cs,k∂rr=Rs,k=jk.


### 2.3.1. Scale Coupling

For eachx∈Ωk, the number of spherical particles per unit cross section area and in the range x,x+dx is(10)dnk=1-ϵk-ϵf,k4/3πRs,k3dx.The solid matrix occupies a volume fraction 1-ϵk-ϵf,k, because apart from electrolyte, there can be some amount of inert material(filler) with the volume fraction ϵf,k [15]. According to (10), the number of spherical particles per unit volume at x∈Ωk is (11)nk=1-ϵk-ϵf,k4/3πRs,k3.All such particles are assumed to have the same radial distribution of concentration cs,k. In particular, we denote the concentration at the particle surface as(12)cs,surf,kx,t=cs,kRs,k,x,t.The total surface area of these particles per unit volume(13)ak=nk4πRs,k2=31-ϵk-ϵf,kRs,kis 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 Rs,k from (13).


## 2.3.1. Scale Coupling

For eachx∈Ωk, the number of spherical particles per unit cross section area and in the range x,x+dx is(10)dnk=1-ϵk-ϵf,k4/3πRs,k3dx.The solid matrix occupies a volume fraction 1-ϵk-ϵf,k, because apart from electrolyte, there can be some amount of inert material(filler) with the volume fraction ϵf,k [15]. According to (10), the number of spherical particles per unit volume at x∈Ωk is (11)nk=1-ϵk-ϵf,k4/3πRs,k3.All such particles are assumed to have the same radial distribution of concentration cs,k. In particular, we denote the concentration at the particle surface as(12)cs,surf,kx,t=cs,kRs,k,x,t.The total surface area of these particles per unit volume(13)ak=nk4πRs,k2=31-ϵk-ϵf,kRs,kis 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 Rs,k from (13).


## 2.4. Electrochemical Interactions

Equations (1) and (6) determine the chemical state of the cell provided that the fluxes j1 and j3 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 Section2.4, the equations do not contain partial derivatives with respect to t and r. For any function f:Ωk×J→R, k∈{1,2,3}, we therefore simplify the notation ∂f/∂x to f′.


### 2.4.1. 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 (14a)σkeffϕs,k′′=akFjkinΩk×J,k∈1,3with the boundary conditions(14b)-σ1effϕs,1′xa=Iapp,(14c)-σ1effϕs,1′xb=0,(14d)-σ3effϕs,3′xc=0,(14e)-σ3effϕs,3′xd=Iapp,(14f)ϕs,3xd=0,where ϕs,k[V] is the solid phase electrical potential,F[Cmol-1] is Faraday’s constant,Iapp[Am-2] is the current density applied to the electrode (Iapp>0 corresponds to charging and Iapp<0 to discharging),Fjk[Am-2] denotes the charge flux in terms of Li+ ions from the internal surface of the porous material into the electrolyte per unit surface area,σkeff[Sm-1] is the effective electronic conductivity defined as (15)σkeff=1-ϵk-ϵf,kσk,σk[Sm-1] is the electronic conductivity.


### 2.4.2. Electrical Potential in the Electrolyte

Based on [15], the charge continuity equation in the electrolyte is given by (16a)κkeffϕe,k′′=-akFjk+2RTF1-t+0κkeffce,k′ce,k′inΩk×J,k∈1,2,3with the boundary conditions(16b)-κ1effϕe,1′xa=0,(16c)-κ1effϕe,1′xb=-κ2effϕe,2′xb=Iapp,(16d)-κ2effϕe,2′xc=-κ3effϕe,3′xc=Iapp,(16e)-κ3effϕe,3′xd=0,(16f)ϕe,1xb=ϕe,2xb,(16g)ϕe,2xc=ϕe,3xc,where ϕe,k[V] is the electrical potential in the electrolyte,R[Jmol-1K-1] is the universal gas constant,T[K] is the temperature,κkeff[Sm-1] is the effective ionic conductivity of the electrolyte given by κkeff=τkκk that accounts for the tortuous path of the porous medium, where κk[Sm-1].For the reference ionic conductivityκk, we use the empirical correlation reported in [15]:(17)κk=4.1253·10-2+5.007·10-4ce,k-4.7212·10-7ce,k2+1.5094·10-10ce,k3-1.6018·10-14ce,k4,k∈1,2,3.


### 2.4.3. Butler-Volmer Reaction Kinetics

Lithium is conserved and its fluxes from the solid particles and into the electrolyte are both equal tojk. The charge transfer Fjk is proportional to the mass transfer and is subject to Butler-Volmer reaction kinetics [6, 7, 15, 21] in the form (18a)jk=δkexp⁡αa,kFRTηk-exp⁡-αc,kFRTηk,k=1,3,with(18b)δk=Kkcs,max,k-cs,surf,k1/2cs,surf,k1/2ce,k1/2.The meanings of the symbols in (18a) and (18b) are the following: Kk[mol-1/2m5/2s-1] is the reaction rate coefficient,αa,k and αc,k[1] are the anodic and cathodic transfer coefficients of electrochemical reaction,cs,max,k[molm-3] is the saturated concentration of Li+ ions in the solid phase,ηk[V] is the intercalation overpotential described as(19)ηk=ϕs,k-ϕe,k-Uk.In (19), Uk[V] is the open circuit potential determined by the following empirical correlations: (20)U1=-4.875+5.839θ1-1.507θ13+0.531θ15θ1-1.005,U3=0.15-0.10θ3+0.00778θ3,where θk=cs,surf,k/cs,max,k, k=1,3.


## 2.4.1. 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 (14a)σkeffϕs,k′′=akFjkinΩk×J,k∈1,3with the boundary conditions(14b)-σ1effϕs,1′xa=Iapp,(14c)-σ1effϕs,1′xb=0,(14d)-σ3effϕs,3′xc=0,(14e)-σ3effϕs,3′xd=Iapp,(14f)ϕs,3xd=0,where ϕs,k[V] is the solid phase electrical potential,F[Cmol-1] is Faraday’s constant,Iapp[Am-2] is the current density applied to the electrode (Iapp>0 corresponds to charging and Iapp<0 to discharging),Fjk[Am-2] denotes the charge flux in terms of Li+ ions from the internal surface of the porous material into the electrolyte per unit surface area,σkeff[Sm-1] is the effective electronic conductivity defined as (15)σkeff=1-ϵk-ϵf,kσk,σk[Sm-1] is the electronic conductivity.


## 2.4.2. Electrical Potential in the Electrolyte

Based on [15], the charge continuity equation in the electrolyte is given by (16a)κkeffϕe,k′′=-akFjk+2RTF1-t+0κkeffce,k′ce,k′inΩk×J,k∈1,2,3with the boundary conditions(16b)-κ1effϕe,1′xa=0,(16c)-κ1effϕe,1′xb=-κ2effϕe,2′xb=Iapp,(16d)-κ2effϕe,2′xc=-κ3effϕe,3′xc=Iapp,(16e)-κ3effϕe,3′xd=0,(16f)ϕe,1xb=ϕe,2xb,(16g)ϕe,2xc=ϕe,3xc,where ϕe,k[V] is the electrical potential in the electrolyte,R[Jmol-1K-1] is the universal gas constant,T[K] is the temperature,κkeff[Sm-1] is the effective ionic conductivity of the electrolyte given by κkeff=τkκk that accounts for the tortuous path of the porous medium, where κk[Sm-1].For the reference ionic conductivityκk, we use the empirical correlation reported in [15]:(17)κk=4.1253·10-2+5.007·10-4ce,k-4.7212·10-7ce,k2+1.5094·10-10ce,k3-1.6018·10-14ce,k4,k∈1,2,3.


## 2.4.3. Butler-Volmer Reaction Kinetics

Lithium is conserved and its fluxes from the solid particles and into the electrolyte are both equal tojk. The charge transfer Fjk is proportional to the mass transfer and is subject to Butler-Volmer reaction kinetics [6, 7, 15, 21] in the form (18a)jk=δkexp⁡αa,kFRTηk-exp⁡-αc,kFRTηk,k=1,3,with(18b)δk=Kkcs,max,k-cs,surf,k1/2cs,surf,k1/2ce,k1/2.The meanings of the symbols in (18a) and (18b) are the following: Kk[mol-1/2m5/2s-1] is the reaction rate coefficient,αa,k and αc,k[1] are the anodic and cathodic transfer coefficients of electrochemical reaction,cs,max,k[molm-3] is the saturated concentration of Li+ ions in the solid phase,ηk[V] is the intercalation overpotential described as(19)ηk=ϕs,k-ϕe,k-Uk.In (19), Uk[V] is the open circuit potential determined by the following empirical correlations: (20)U1=-4.875+5.839θ1-1.507θ13+0.531θ15θ1-1.005,U3=0.15-0.10θ3+0.00778θ3,where θk=cs,surf,k/cs,max,k, k=1,3.


## 3. Integral Solution of Electrochemical Interactions

For a given timet, the flux jk, k∈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 ϕe,k, ϕs,k, ηk, and jk, 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 ηk and Jk, where by Jk we denote the integrated mass transfer function defined by(21)Jkx=∫xl,kxjkξdξ,∀x∈xl,k,xr,k,k∈1,3,and xl,k, xr,k represent the left and right boundary coordinates of Ωk, respectively.


### 3.1. Derivation of the Differential Equations

In order to derive the ODE forηk, we first integrate (14a) with respect to x over the interval (xl,k,x)⊂Ωk, k∈{1,3}. Using the boundary conditions (14b) and (14d), we obtain (22a)ϕs,1′=a1Fσ1effJ1-Iappσ1eff,∀x∈Ω1,(22b)ϕs,3′=a3Fσ3effJ3,∀x∈Ω3.Similarly, by integrating (16a) for k∈{1,2,3} and using the boundary conditions (16b), (16c), and (16d), we arrive at(22c)ϕe,1′=-a1Fκ1effJ1+2RTF1-t+0ce,1′ce,1,∀x∈Ω1,(22d)ϕe,2′=2RTF1-t+0ce,2′ce,2-Iappκ2eff,∀x∈Ω2,(22e)ϕe,3′=-a3Fκ3effJ3+2RTF1-t+0ce,3′ce,3-Iappκ3eff,∀x∈Ω3.Subtractions of (22c) from (22a) and (22e) from (22b) allow expressing the derivative of ηk defined by (19) as (23a)η1′=a1FJ1σ1eff+κ1effσ1effκ1eff-Iappσ1eff-2RTF1-t+0ce,1′ce,1-U1′,(23b)η3′=a3FJ3σ3eff+κ3effσ3effκ3eff+Iappκ3eff-2RTF1-t+0ce,3′ce,3-U3′.The ODE forJk is given directly by the Butler-Volmer reaction kinetics (18a) in the form(24)Jk′=δkexp⁡αa,kFRTηk-exp⁡-αc,kFRTηk,k=1,3.Differentiating (24) with respect to x and plugging in ηk from (23a) and (23b), respectively, lead to a single second-order ODE for Jk. The corresponding boundary conditions follow from the definition of Jk given by (21) and evaluation of (22a), (22b) at x=xb and x=xc together with the boundary conditions given by (14c) and (14e). They read (25a)J1xa=0,J1xb=Iappa1F,(25b)J3xc=0,J3xd=-Iappa3F,and they complete two well posed problems for the unknowns J1 and J3.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.


### 3.2. ODE System in General Form

Assume a general coefficient form of ODEs ((23a) and (23b)) and (24) in the following compact form: (26a)η′x=αxJx+βx,(26b)J′x=Bηx,x,where(26c)Bηx,x=δxexp⁡αaγxηx-exp⁡-αcγxηx,for all x∈(xl,xr) with the following boundary conditions:(26d)Jxl=0,(26e)Jxr=ε,where xl=xl,k, xr=xr,k, αa=αa,k, αc=αc,k, the coefficients δ=δk are given by (18b), and α=αk, β=βk, and ε=εk correspond to the following domain-defined coefficients (denoted by the indices k=1 and k=3): (27)α1=a1Fσ1+κ1effσ1κ1eff,α3=a3Fσ3+κ3effσ3κ3eff,β1=-Iappσ1-2RTF1-t+0ce,1′ce,1-U1′,β3=Iappκ3eff-2RTF1-t+0ce,3′ce,3-U3′,γ1=FRT,γ3=FRT,ε1=Iappa1F,ε3=-Iappa3F.


### 3.3. Derivation of the Integral Equation

Integrating (26a), (26b), (26c), (26d), and (26e) from xl to x∈[xl,xr], we obtain (28a)ηx=λ+∫xlxαζJζ+βζdζ,(28b)Jx=∫xlxBηζ,ζdζ,where, in (28a), λ=η(xl) 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:(29)ηx=λ+∫xlxαζ∫xlζBηξ,ξdξ+βζdζ,or(30)Jx=∫xlxBλ+∫xlζαξJξ+βξdξ,ζdζ.


### 3.4. Equation forλ

The integral equation given by (30) allows determining the value of the unknown integration constant λ since the remaining boundary condition given by (26e) implies(31)ε=∫xlxrBλ+∫xlζαξJξ+βξdξ,ζdζ.Under the assumption of constant temperatureT, γ becomes constant in (26c). As a result, (26c) and (31) can be combined to obtain(32)ε=ΛαaIa-Λ-αcIc,where Λ=exp(γλ) denotes the term with the unknown parameter λ and the coefficients Ia and Ic read as (33a)Ia=∫xlxrδζexp⁡αaγ∫xlζαξJξ+βξdξdζ,(33b)Ic=∫xlxrδζexp⁡-αcγ∫xlζαξJξ+βξdξdζ.If αa=αc, (32) is easily resolved by(34a)Λαa=Λαc=ε+ε2+4IaIc2Ia;that is,(34b)λ=1γln⁡ε+ε2+4IaIc2Ia.In general, (32) is a highly nonlinear equation that needs to be solved numerically.


### 3.5. Iterative Scheme for Solving the Integral Equation

We propose the following iterative scheme for solving the integral equation given by (30):(35)J~n+1x=∫xlxBλn+∫xlζαξJnξ+βξdξ,ζdζ,(36)Jn+1x=1-ωJnx+ωJ~n+1x,where λn is the solution of (32) which has to be updated at every iteration step for given Jn, n=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 J0≡0 in the first call to the iterative solver and the previously calculated value of J in the subsequent calls. The iterative process is terminated when the norm of J~n+1 is below a given threshold ϑ. The values of ω and ϑ are discussed further in Section 5.2.


### 3.6. 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 ofce,k and cs,surf,k and the value of Iapp to solve (26a), (26b), (26c), (26d), and (26e) using the iteration scheme given by (35) to obtain the integrated charge fluxes J1 and J3.(2)
Compute the intercalation overpotentialsη1 and η3 using (28a).(3)
Compute the fluxesj1=J1′ and j3=J3′ from (26b); that is, j1(x)=B(η1(x),x), and j3(x)=B(η3(x),x), respectively.(4)
Computeϕs,3 by integrating (22b) and using the boundary condition given by (14f).(5)
Computeϕe,3 by integrating (22e) and using the boundary condition that follows from (19) as (37)ϕe,3xc=ϕs,3xc-η3xc-U3xc.(6)
Computeϕe,2 by integrating (22d) and using the boundary condition given by (16g).(7)
Computeϕe,1 by integrating (22c) and using the boundary condition given by (16f).(8)
Computeϕs,1 by integrating (22a) and using the boundary condition that follows from (19) as(38)ϕs,1xb=ϕe,1xb-η1xb-U1xb.(9)
The value ofϕs,1(xa) represents the external apparent voltage of the cell.


## 3.1. Derivation of the Differential Equations

In order to derive the ODE forηk, we first integrate (14a) with respect to x over the interval (xl,k,x)⊂Ωk, k∈{1,3}. Using the boundary conditions (14b) and (14d), we obtain (22a)ϕs,1′=a1Fσ1effJ1-Iappσ1eff,∀x∈Ω1,(22b)ϕs,3′=a3Fσ3effJ3,∀x∈Ω3.Similarly, by integrating (16a) for k∈{1,2,3} and using the boundary conditions (16b), (16c), and (16d), we arrive at(22c)ϕe,1′=-a1Fκ1effJ1+2RTF1-t+0ce,1′ce,1,∀x∈Ω1,(22d)ϕe,2′=2RTF1-t+0ce,2′ce,2-Iappκ2eff,∀x∈Ω2,(22e)ϕe,3′=-a3Fκ3effJ3+2RTF1-t+0ce,3′ce,3-Iappκ3eff,∀x∈Ω3.Subtractions of (22c) from (22a) and (22e) from (22b) allow expressing the derivative of ηk defined by (19) as (23a)η1′=a1FJ1σ1eff+κ1effσ1effκ1eff-Iappσ1eff-2RTF1-t+0ce,1′ce,1-U1′,(23b)η3′=a3FJ3σ3eff+κ3effσ3effκ3eff+Iappκ3eff-2RTF1-t+0ce,3′ce,3-U3′.The ODE forJk is given directly by the Butler-Volmer reaction kinetics (18a) in the form(24)Jk′=δkexp⁡αa,kFRTηk-exp⁡-αc,kFRTηk,k=1,3.Differentiating (24) with respect to x and plugging in ηk from (23a) and (23b), respectively, lead to a single second-order ODE for Jk. The corresponding boundary conditions follow from the definition of Jk given by (21) and evaluation of (22a), (22b) at x=xb and x=xc together with the boundary conditions given by (14c) and (14e). They read (25a)J1xa=0,J1xb=Iappa1F,(25b)J3xc=0,J3xd=-Iappa3F,and they complete two well posed problems for the unknowns J1 and J3.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.


## 3.2. ODE System in General Form

Assume a general coefficient form of ODEs ((23a) and (23b)) and (24) in the following compact form: (26a)η′x=αxJx+βx,(26b)J′x=Bηx,x,where(26c)Bηx,x=δxexp⁡αaγxηx-exp⁡-αcγxηx,for all x∈(xl,xr) with the following boundary conditions:(26d)Jxl=0,(26e)Jxr=ε,where xl=xl,k, xr=xr,k, αa=αa,k, αc=αc,k, the coefficients δ=δk are given by (18b), and α=αk, β=βk, and ε=εk correspond to the following domain-defined coefficients (denoted by the indices k=1 and k=3): (27)α1=a1Fσ1+κ1effσ1κ1eff,α3=a3Fσ3+κ3effσ3κ3eff,β1=-Iappσ1-2RTF1-t+0ce,1′ce,1-U1′,β3=Iappκ3eff-2RTF1-t+0ce,3′ce,3-U3′,γ1=FRT,γ3=FRT,ε1=Iappa1F,ε3=-Iappa3F.


## 3.3. Derivation of the Integral Equation

Integrating (26a), (26b), (26c), (26d), and (26e) from xl to x∈[xl,xr], we obtain (28a)ηx=λ+∫xlxαζJζ+βζdζ,(28b)Jx=∫xlxBηζ,ζdζ,where, in (28a), λ=η(xl) 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:(29)ηx=λ+∫xlxαζ∫xlζBηξ,ξdξ+βζdζ,or(30)Jx=∫xlxBλ+∫xlζαξJξ+βξdξ,ζdζ.


## 3.4. Equation forλ

The integral equation given by (30) allows determining the value of the unknown integration constant λ since the remaining boundary condition given by (26e) implies(31)ε=∫xlxrBλ+∫xlζαξJξ+βξdξ,ζdζ.Under the assumption of constant temperatureT, γ becomes constant in (26c). As a result, (26c) and (31) can be combined to obtain(32)ε=ΛαaIa-Λ-αcIc,where Λ=exp(γλ) denotes the term with the unknown parameter λ and the coefficients Ia and Ic read as (33a)Ia=∫xlxrδζexp⁡αaγ∫xlζαξJξ+βξdξdζ,(33b)Ic=∫xlxrδζexp⁡-αcγ∫xlζαξJξ+βξdξdζ.If αa=αc, (32) is easily resolved by(34a)Λαa=Λαc=ε+ε2+4IaIc2Ia;that is,(34b)λ=1γln⁡ε+ε2+4IaIc2Ia.In general, (32) is a highly nonlinear equation that needs to be solved numerically.


## 3.5. Iterative Scheme for Solving the Integral Equation

We propose the following iterative scheme for solving the integral equation given by (30):(35)J~n+1x=∫xlxBλn+∫xlζαξJnξ+βξdξ,ζdζ,(36)Jn+1x=1-ωJnx+ωJ~n+1x,where λn is the solution of (32) which has to be updated at every iteration step for given Jn, n=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 J0≡0 in the first call to the iterative solver and the previously calculated value of J in the subsequent calls. The iterative process is terminated when the norm of J~n+1 is below a given threshold ϑ. The values of ω and ϑ are discussed further in Section 5.2.


## 3.6. 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 ofce,k and cs,surf,k and the value of Iapp to solve (26a), (26b), (26c), (26d), and (26e) using the iteration scheme given by (35) to obtain the integrated charge fluxes J1 and J3.(2)
Compute the intercalation overpotentialsη1 and η3 using (28a).(3)
Compute the fluxesj1=J1′ and j3=J3′ from (26b); that is, j1(x)=B(η1(x),x), and j3(x)=B(η3(x),x), respectively.(4)
Computeϕs,3 by integrating (22b) and using the boundary condition given by (14f).(5)
Computeϕe,3 by integrating (22e) and using the boundary condition that follows from (19) as (37)ϕe,3xc=ϕs,3xc-η3xc-U3xc.(6)
Computeϕe,2 by integrating (22d) and using the boundary condition given by (16g).(7)
Computeϕe,1 by integrating (22c) and using the boundary condition given by (16f).(8)
Computeϕs,1 by integrating (22a) and using the boundary condition that follows from (19) as(38)ϕs,1xb=ϕe,1xb-η1xb-U1xb.(9)
The value ofϕs,1(xa) represents the external apparent voltage of the cell.


## 4. 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 N uniformly spaced nodes x0,…,xN-1 and the positions of the individual grid nodes determine their correspondence to the domains Ωk, k∈1,2,3. This grid is also used for numerical integration in the algorithm described in Section 3. For each node xi∈Ωk, k∈1,3, another mesh of M nodes exists that discretizes the interval 0,Rs,k, 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 xi with a constant time step Δt. After a given number of time steps, the integral solver is called, which updates the values of jk for each xi. For the time period between the updates of jk, 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 Δt to be used as the time step for the diffusion in the electrolyte.


## 5. 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.Table 1
Model parameters based on literature values [15, 18].
(a)
Global parameters
Parameter
Unit
Value
t+0
[
1
]
0.363
F
[
C
m
o
l
-
1
]
96487
T
[
K
]
298.15
xa
[
μ
m
]
0
xb
[
μ
m
]
80
xc
[
μ
m
]
105
xd
[
μ
m
]
193
De
[
m
2
s
-
1
]
7.5
×
1
0
-
10
(b)
Domain-specific parameters
Parameter
Unit
Positive electrodeΩ1
SeparatorΩ2
Negative electrodeΩ3
k
1
2
3
σk
[
S
m
-
1
]
100
100
ϵf,k
[
1
]
0.025
0.0326
ϵk
[
1
]
0.385
0.724
0.485
bruggk
[
1
]
4
4
4
αa,k
[
1
]
0.5
0.5
αc,k
[
1
]
0.5
0.5
cs,max,k
[
k
m
o
l
m
-
3
]
51.554
30.555
Rs,k
[
μ
m
]
2
2
Kk
[
m
o
l
-
0.5
m
2.5
s
-
1
]
2.3444
×
1
0
-
11
5.0307
×
1
0
-
11
Ds,k
[
m
2
s
-
1
]
1
×
1
0
-
14
3.9
×
1
0
-
14
cs,k0
[
k
m
o
l
m
-
3
]
0.4955
⋅
c
s
,
m
a
x
,
1
0.8551
⋅
c
s
,
m
a
x
,
3
ce,k0
[
k
m
o
l
m
-
3
]
1
1
1


### 5.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.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.
(a)
(b)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 t=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.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.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].Figure 4
Comparison of the spatial profiles ofce 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 t. Discharge current 1C.Figure 5
Comparison of the spatial profiles ofce at the end of 1C discharge cycle (t=3500s). The solid line represents our results; the thick dashed line is the result from [3].


### 5.2. 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 Section4. 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 N=20 to N=800. For each value of N, uniform grids from M=20 to M=400 were utilized for diffusion inside the spherical particles. For the 1C applied current and M fixed, the solution is almost independent of the value of N provided that N≥50. For different values of M, there are slight differences mostly in the final part of the voltage curve, as demonstrated in Figure 6. As M 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 N also affects the shape of the voltage curve. For the 2C applied current, the differences in the results become negligible for any combination of M and N satisfying M>200, N>200.The simulation of the whole 1C discharge cycle with M=N=100, Δt=10-1s, ω=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 M and N below approximately 50 nodes (and losing accuracy as seen in Figure 6) is not necessarily beneficial for the computational time.(ii)
As explained in Section4, the integral solver need not be called in every time step. We performed some computations with a fixed (and unnecessarily small) time step Δt=10-3s. The integral solver was set to be called in the intervals m·Δt, for several different values of m∈N. In addition, the diffusion in electrolyte was also solved using the time step m·Δt. 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 ω.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. N=100, Δt=10-1s, ω=0.05, and ϑ=10-13. Integral solver called in every time step.Figure 7
The number of iterations of the integral solver depending on the frequency of its calls. Discharge current 1C: N=100, M=50, Δt=10-3s, ω=0.05, and ϑ=10-13.
(a)
(b)Figure 8
Computational times of one discharge cycle depending on the frequency of the calls of the integral solver. Discharge current 1C; time step Δt=10-3s. Computed on an Intel i5-2500 @ 3.3 GHz CPU.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 Figure9 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.Figure 9
Algorithm stability test during 1C constant current discharge/charge cycling within the voltage range from 3.4V to 4.5V. A total of 1043 cycles in the physical time period of 4.41×106s (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: N=100, M=50, Δt=10-1s, ω=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.


## 5.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.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.
(a)
(b)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 t=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.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.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].Figure 4
Comparison of the spatial profiles ofce 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 t. Discharge current 1C.Figure 5
Comparison of the spatial profiles ofce at the end of 1C discharge cycle (t=3500s). The solid line represents our results; the thick dashed line is the result from [3].


## 5.2. 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 Section4. 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 N=20 to N=800. For each value of N, uniform grids from M=20 to M=400 were utilized for diffusion inside the spherical particles. For the 1C applied current and M fixed, the solution is almost independent of the value of N provided that N≥50. For different values of M, there are slight differences mostly in the final part of the voltage curve, as demonstrated in Figure 6. As M 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 N also affects the shape of the voltage curve. For the 2C applied current, the differences in the results become negligible for any combination of M and N satisfying M>200, N>200.The simulation of the whole 1C discharge cycle with M=N=100, Δt=10-1s, ω=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 M and N below approximately 50 nodes (and losing accuracy as seen in Figure 6) is not necessarily beneficial for the computational time.(ii)
As explained in Section4, the integral solver need not be called in every time step. We performed some computations with a fixed (and unnecessarily small) time step Δt=10-3s. The integral solver was set to be called in the intervals m·Δt, for several different values of m∈N. In addition, the diffusion in electrolyte was also solved using the time step m·Δt. 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 ω.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. N=100, Δt=10-1s, ω=0.05, and ϑ=10-13. Integral solver called in every time step.Figure 7
The number of iterations of the integral solver depending on the frequency of its calls. Discharge current 1C: N=100, M=50, Δt=10-3s, ω=0.05, and ϑ=10-13.
(a)
(b)Figure 8
Computational times of one discharge cycle depending on the frequency of the calls of the integral solver. Discharge current 1C; time step Δt=10-3s. Computed on an Intel i5-2500 @ 3.3 GHz CPU.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 Figure9 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.Figure 9
Algorithm stability test during 1C constant current discharge/charge cycling within the voltage range from 3.4V to 4.5V. A total of 1043 cycles in the physical time period of 4.41×106s (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: N=100, M=50, Δt=10-1s, ω=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.


## 6. 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–9].



---

*Source: 3530975-2018-02-07.xml*