Journal of The Electrochemical

Society

You may also like

- Report on the Electrolytic Industries for the
LIONSIMBA: A Matlab Framework Based on a Year 2001
Pankaj Arora and Venkat Srinivasan
Finite Volume Model Suitable for Li-Ion Battery - Analytical Model Based Prediction of
State-of-Charge (SoC) of a Lithium-Ion
Design, Simulation, and Control Cell under Time-Varying
Charge/Discharge Currents
Mohammad Parhizi, Manan Pathak and
To cite this article: Marcello Torchio et al 2016 J. Electrochem. Soc. 163 A1192 Ankur Jain

- State-of-the-art of battery state-of-charge

determination
V Pop, H J Bergveld, P H L Notten et al.

View the article online for updates and enhancements.

This content was downloaded from IP address 201.71.169.193 on 26/10/2023 at 20:40

A1192 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)
0013-4651/2016/163(7)/A1192/14/$33.00 © The Electrochemical Society

LIONSIMBA: A Matlab Framework Based on a Finite Volume

Model Suitable for Li-Ion Battery Design, Simulation, and Control
Marcello Torchio,a Lalo Magni,a R. Bhushan Gopaluni,b Richard D. Braatz,c
and Davide M. Raimondoa,z
a University of Pavia, 27100 Pavia, Italy
b University of British Columbia, Vancouver, BC V6T 1Z3, Canada
c Massachusetts Institute of Technology, Cambridge, Massachusetts 02142, USA

Consumer electronics, wearable and personal health devices, power networks, microgrids, and hybrid electric vehicles (HEVs) are
some of the many applications of lithium-ion batteries. Their optimal design and management are important for safe and profitable
operations. The use of accurate mathematical models can help in achieving the best performance. This article provides a detailed
description of a finite volume method (FVM) for a pseudo-two-dimensional (P2D) Li-ion battery model suitable for the development
of model-based advanced battery management systems. The objectives of this work are to provide: (i) a detailed description of the
model formulation, (ii) a parametrizable Matlab framework for battery design, simulation, and control of Li-ion cells or battery packs,
(iii) a validation of the proposed numerical implementation with respect to the COMSOL MultiPhysics commercial software and
the Newman’s DUALFOIL code, and (iv) some demonstrative simulations involving thermal dynamics, a hybrid charge-discharge
cycle emulating the throttle of an HEV, a model predictive control of state of charge, and a battery pack simulation.
© 2016 The Electrochemical Society. [DOI: 10.1149/2.0291607jes] All rights reserved.

Manuscript submitted November 18, 2015; revised manuscript received March 22, 2016. Published April 12, 2016.

The increasing demand for portable devices (e.g., smartphones) rected to the handling of interface boundary conditions in the three pri-
and hybrid electric vehicles (HEVs) calls for the design and manage- mary sections of the battery - positive and negative electrodes and the
ment of storage devices of high power density and reduced size and separator. Due to possible discontinuities between adjacent sections,
weight. During the many decades of research, different chemistries a mishandling of such conditions may lead to physically inconsistent
of batteries have been developed, such as Nickel Cadmium (NiCd), solutions. Due to its intrinsic properties, the finite volume method has
Nickel Metal Hydride (NiMH), Lead Acid and Lithium ion (Li-ion) been chosen to easily deal with these particular interface conditions.
and Lithium ion Polymer (Li-Poly) (e.g., see Refs. 1–4). Among elec- Finally, based on the proposed finite volume discretization, we provide
trochemical accumulators, Li-ion batteries provide one of the best the Li-ION SIMulation BAttery Toolbox (LIONSIMBA), a set of fully
tradeoff in terms of power density, low weight, cell voltage, and customizable Matlab functions suitable for simulating the dynamic be-
low self-discharge.5 Mathematical models can support the design of havior of Li-ion batteries. These functions are freely downloadable
new batteries as well as the development of new advanced battery from the website http://sisdin.unipv.it/labsisdin/lionsimba.php. This
management systems (ABMS).6–8 According to the literature, math- article describes the features of our software. The user can implement
ematical models for Li-ion battery dynamics fall within two main his/her own custom-defined control algorithm to test different ABMS
categories: Equivalent Circuit Models (ECMs) and Electrochemical strategies, simulate cell behavior, optimize manufacturing parameters
Models (EMs). ECMs use only electrical components to model the or test battery packs composed of series-connected cells. The pack-
dynamic behavior of the battery. ECMs include (i) the Rint model age also allows the ready implementation of algorithms to estimate
where only a resistance and a voltage source are used to model the indexes such as the State of Charge (SOC) and the State of Health
battery, (ii) the RC model (introduced by the company SAFT9 ) where (SOH). The SOC is an important property of batteries that quantifies
capacitor dynamics have been added to the Rint model,10 and (iii) the the amount of remaining charge (e.g., Ref. 15) and can be used to pre-
Thevenin model, which is an extension of the RC model (e.g., see vent damage, ensure safety, and minimize charging time.16 The SOH
Refs. 11, 12 and references therein). In contrast, EMs explicitly rep- index measures the ability of the battery to store and deliver electrical
resent the chemical processes that take place in the battery. While energy; similar to the SOC, estimation-based approaches are used to
ECMs have the advantage of simplicity, EMs are more accurate due predict the value of the SOH (e.g., see Refs. 17–19). The SOH tracks
to their ability to describe detailed physical phenomena.13 The most the long-term changes in a battery and its knowledge can help ABMS
widely used EM in the literature is the porous electrode theory-based to anticipate problems through online fault diagnosis while providing
pseudo-two-dimensional (P2D) model,14 which is described by a set charging profiles to slow down the battery aging. The package comes
of tightly coupled and highly nonlinear partial differential-algebraic with the experimental parameters of the battery reported in Ref. 20.
equations (PDAEs). In order to exploit the model for simulation and An initialization file allows changes in battery and simulator param-
design purposes, the set of PDAEs are reformulated as a set of ordinary eters. The simulator works under Matlab using IDA21 to solve the
differential-algebraic equations (DAEs). The model reformulation is set of resulting DAEs with a good trade-off between accuracy and
very challenging to carry out in a way that is simultaneously compu- computational time.
tationally efficient and numerically stable for a wide range of battery The battery model and its numerical implementation is described
parameters and operating conditions. To the authors’ best knowledge, first. Then the proposed framework is validated with respect to
no publication is available in the literature that provides a detailed step- the results obtained using the COMSOL MultiPhysics commercial
by-step description of the numerical implementation of the P2D model software22 and the Newman’s Fortran code named DUALFOIL.23
or a freely available Matlab framework suitable for simulation, design, Finally, to demonstrate the effectiveness of the proposed software,
and development of ABMS for Li-ion batteries. In this article, start- thermal dynamics, model predictive control of state of charge, hybrid
ing from the P2D model, a computationally efficient and numerically charge-discharge cycles, and a battery pack of series-connected cells
stable finite volume DAE formulation is described in detail in order to are simulated. The toolbox is equipped with all the Matlab source files
facilitate implementation by the reader, while also addressing potential able to reproduce the simulations presented in this article.
pitfalls and relative loopholes. Boundary conditions used to enforce
physical meaningfulness of the system are thoroughly discussed and
their numerical implementation is explained. Particular attention is di-
Battery Model
The P2D model consists of coupled nonlinear PDAEs for the
z
E-mail: [email protected] conservation of mass and charge in the three sections of the
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1193

Table I. Li-ion P2D model governing equations.

Current Collectors, i ∈ {a, z} Boundary Conditions



  2 (t)
Iapp −λa ∂ T∂(x,t) = h(Tref − T (x, t))
x x=0
ρi C p,i ∂ T∂t
(x,t)
= ∂
λi ∂ T∂(x,t) + 
∂x x σeff,i ∂ T (x,t) 
−λz ∂ x  = h(T (x, t) − Tref )
x=L
Positive and Negative
  i ∈ { p, n}
Electrodes, 
∂ce (x,t) 
i ∂ce∂t(x,t) = ∂
∂x Deff,i ∂ce∂(x,t)
x + ai (1 − t+ ) j(x, t) ∂ x x=x̂ ,x̂ =0
avg 0 n
∂cs (x,t)
∂t = −3 j(x,t)
R p,i
R
cs∗ (x, t) − cs (x, t) =
avg j(x,t)
− Dsp,i
eff,i
5

∂s (x,t) 
  σeff,i ∂x  = −Iapp (t)
∂
σeff,i ∂∂s (x,t) = ai F j(x, t) x=x̂0 ,x̂n
∂x x ∂s (x,t) 
σeff,i ∂x  =0
x=x̂ ,x̂

p s
    ∂e (x,t) 
 =0
ai F j(x, t) = − ∂∂x κeff,i ∂∂e (x,t)
x + ∂
∂x κeff,i ϒ T (x, t) ∂ ln c∂ex(x,t) ∂x x=x̂0
e (x, t)|x=x̂n = 0
 
 
  −λz ∂ T∂(x,t)
x x=x̂ − = −λ p ∂ T∂(x,t)x x=x̂ +
ρi C p,i ∂ T∂t
(x,t)
= ∂
λi ∂ T∂(x,t) + Q ohm + Q rxn + Q rev 

0  0
∂x ∂ T (x,t) 
−λn ∂ T∂(x,t)
x
x  −
= −λ z ∂ x  +
   x=x̂n x=x̂n
max − c∗ (x, t))c∗ (x, t) sinh
j(x, t) = 2keff,i ce (x, t)(cs,i F T (x,t) ηi (x, t)
0.5R
s s
ηi (x, t) = s (x, t) − e (x, t) − Ui
Separator, i = s  
∂ce (x,t)  ∂ce (x,t) 
  − Deff, p ∂ x x=x̂ − = − Deff,s ∂ x x=x̂ +
i ∂ce∂t(x,t) = ∂
Deff,i ∂ce∂(x,t)  p  p
∂x x ∂ce (x,t)  ∂ce (x,t) 
− Deff,s ∂ x x=x̂ − = − Deff,n ∂ x x=x̂ +
 s
 s
∂e (x,t)  ∂e (x,t) 
    −κeff, p ∂x  − = −κeff,s ∂x  +
0 = − ∂∂x κeff,i ∂∂e (x,t)
x + ∂
∂x κeff,i T (x, t)ϒ ∂ ln c∂ex(x,t) ∂e (x,t) 

x=x̂ p
∂e (x,t) 
x=x̂  p

−κeff,s ∂x  − = −κeff,n ∂x  +
x=x̂ x=x̂
 s
 s
∂ T (x,t)  ∂ T (x,t) 
  −λ p ∂ x x=x̂ − = −λs ∂ x x=x̂ +
ρi C p,i ∂ T∂t
(x,t)
= ∂
λi ∂ T∂(x,t) + Q ohm  p  p
∂x x ∂ T (x,t)  ∂ T (x,t) 
−λs ∂ x x=x̂ − = −λn ∂ x x=x̂ +
s s

battery – cathode, separator, and anode – denoted respectively by imated by means of average and surface concentration of the solid
the indexes p, s, and n. The positive and negative current collectors particles,
are denoted by a and z. The index i ∈ S is used to refer to a par-
∂csavg (x, t) j(x, t)
ticular section of the battery, where S := {a, p, s, n, z}. All model = −3 ,
equations are reported in Tables I and II. Variables ce (x, t), csavg (x, t), ∂t Rp
and cs∗ (x, t) ∈ R+ denote the electrolyte concentration, the average R p j(x, t)
concentration in the solid particles, and the surface concentration in cs∗ (x, t) − csavg (x, t) = − .
the solid particles of Li-ions respectively, where time t ∈ R+ and D sp 5
x ∈ R is the spatial direction along which the ions are transported. This reformulation leads to a one-dimensional problem in x by re-
Diffusion inside solid spherical particles with radius R p is described moving the pseudo-second dimension r . Despite the reduced com-
by Fick’s law, putational burden, such approximation could lead to a decrease of
 
∂cs (r, t) 1 ∂ ∂cs (r, t) the prediction accuracy for high rates, short time responses or pulse
= 2 r 2 D sp , [1] currents.26 For these applications, higher-order polynomials or Fick’s
∂t r ∂r ∂r
law are recommended, as discussed in Solid-phase diffusion models
with boundary conditions section.
∂cs (r, t)  ∂cs (r, t)  j(x, t) The electrolyte and solid potential are represented by e (x, t) and
 = 0,  =− s , s (x, t) ∈ R, while T (x, t) and j(x, t) represent the temperature and
∂r r =0 ∂r r =R p Deff the ionic flux. Note that the ionic flux is present only in the positive
where r is the radial direction along which the ions intercalate within and negative electrodes, and not in the separator. The open circuit volt-
the active particles. This model introduces a pseudo-second dimen- age (OCV) is denoted by U while the entropic variation of the OCV
sion (r ). To reduce complexity and computational burden, Refs. 24 is denoted by ∂U ∂T
. The cathode, anode, and separator are composed
and 25 proposed different efficient reformulations for the solid-phase of different materials; for a given section i, different electrolyte diffu-
diffusion equation. As discussed in Ref. 26, according to the particu- sion coefficients Di , solid-phase diffusion coefficients Dis , electrolyte
lar application, different model reformulations can be employed while conductivities κi , porosities i , thermal capacities C p,i , thermal con-
maintaining good accuracy. For low to medium C rates, the diffusion ductivities λi , densities ρi , solid-phase conductivities σi , particle sur-
max
length method27 or the two-term polynomial approximation method face area to volumes ai , maximum solid phase concentrations cs,i ,
are accurate. At high C rates, higher-order polynomial approximations overpotentials ηi , and particle radiuses R p,i can be defined. The terms
or the Pseudo Steady State (PSS)28 approximation can be employed. R and F are the universal gas constant and the Faraday constant, re-
For more details, refer to Ref. 26 and the references therein. spectively, with t+ representing the transference number. The applied
In the two-term polynomial approximation, concentration profiles current density is Iapp (t), and Tref denotes the environment tempera-
inside the particle are assumed to be quadratic in r and 1 is approx- ture. In order to take into account the properties of different materials
A1194 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

Table II. Additional equations.

Open Circuit Potential (Thermal dependence)

∂U p
U p = U p,ref + (T (x, t) − Tref ) ∂ T | T ref

U n = U n,ref + (T (x, t) − Tref ) ∂U

∂T | T ref
n

Entropy change
 
∂U p  0.199521039−0.928373822θ p +1.364550689000003θ 2p −0.6115448939999998θ 3p
∂T  = −0.001
T ref 1−5.661479886999997θ p +11.47636191θ 2p −9.82431213599998θ 3p +3.048755063θ 4p
⎛ ⎞
⎜
0.005269056 + 3.299265709θ n − 91.79325798θ 2n + 1004.911008θ 3n − 5812.278127θ 4n + ⎟
0.001⎝ ⎠
 19329.7549θ 5n − 37147.8947θ 6n + 38379.18127θ 7n − 16515.05308θ 8n
∂U n 
∂ T T = ⎛ ⎞
ref
⎜
1 − 48.09287227θ n + 1017.234804θ 2n − 10481.80419θ 3n + 59431.3θ 4n − ⎟
⎝ ⎠
195881.6488θ 5n + 374577.3152θ 6n − 385821.1607θ 7n + 165705.8597θ 8n
Open Circuit Potential (Reference value)
−4.656+88.669θ 2p −401.119θ 4p +342.909θ 6p −462.471θ 8p +433.434θ 10
U p,ref = p
−1+18.933θ 2p −79.532θ 4p +37.311θ 6p −73.083θ 8p +95.96θ 10
p

U n,ref = 0.7222 + 0.1387θ n + 0.029θ 0.5

n −
0.0172
θn + 0.0019
+ 0.2808e0.9−15θn − 0.7984e0.4465θn −0.4108
θ 1.5
n
∗ (x,t)
cs,
θp = p
max
cs, p
∗ (x,t)
cs,n
θn = max
cs,n

Heat source terms (Anode and Cathode)

 2  2
Q ohm = σ eff,i ∂∂s (x,t) + κeff,i ∂∂e (x,t) (1 − t+ ) ∂ ln c∂ex(x,t) ∂∂e (x,t)
2κeff,i RT(x,t)
x x + F x , i ∈ { p, n}

Q rxn = F ai j(x, t) ηi (x, t), i ∈ { p, n}


∂U i 
Q rev = F ai j(x, t) T (x, t) ∂ T T , i ∈ { p, n}
ref

Heat Source terms (Separator)

 2
Q ohm = κeff,i ∂∂e (x,t) (1 − t+ ) ∂ ln c∂ex(x,t) ∂∂e (x,t)
2κ RT(x,t)
x + eff,iF x , i =s

Various Coefficients
bruggi −4.43− 54 −0.22×10−3 ce (x,t)
Deff,i = i × 10−4 × 10 T (x,t)−229−5×10−3 ce (x,t)

⎛ ⎞2
−10.5 + 0.668 · 10−3 · ce (x, t) + 0.494 · 10−6 ce2 (x, t)+
⎜ ⎟
× 10−4 × ce (x, t) ⎜ ⎟
bruggi −5 −10 2
κeff,i = i ⎝ (0.074 − 1.78 × 10 ce (x, t) − 8.86 × 10 ce (x, t))T (x, t)+ ⎠
(−6.96 × 10−5 + 2.8 × 10−8 ce (x, t))T 2 (x, t)
ki  
E 1 1
− Ra T (x,t) − Tref
keff,i = ki e
Dis  
E 1 1
− aR T (x,t) − Tref
Dseff,i = Dis e

σ eff,i = σi (1 − i −  f,i )
2(1−t+ )R
ϒ := F

used in the battery, effective diffusion and conductivity coefficients side, zero-flux conditions are imposed. Within the battery, interface
are evaluated according to the Bruggeman’s theory, with “eff” suffixes conditions are imposed across the different materials. In order to get a
representing effective values of suchcoefficients. The thickness of the more detailed description of the conductivity and diffusion phenom-
overall battery is L, where L = i li and li , represents the length ena inside the electrolyte, all the related coefficients are determined
of each battery section. Due to physical constraints, it is necessary to as a function of ce and T , as discussed in Ref. 29.
impose (i) zero-flux boundary conditions for the ce diffusion equa- Excessive heat generation may lead to performance degradation
tion at the two ends of the battery, (ii) Newton’s cooling law for the and, in extreme cases, thermal runaway of the cell.30,31 In order
dissipation of heat in the system, and (iii) null-flux conditions for s to address these possible safety issues, thermal dynamics are in-
at the interface between electrodes and the separator as well as the cluded with the set of conservation equations describing the sys-
enforcement of Ohm’s law at the end of the electrodes. Given that tem. The thermal equations include different source terms, which are
only potential differences are measurable, without loss of generality, the ohmic, reversible, and reaction generation rates Q ohm , Q rev , and
e can be set to zero at the end of the anode. Similarly, on the cathode Q rxn , respectively.32 The ohmic generation rate takes into account heat
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1195

generated as a consequence of the motion of Li-ions in the solid/liquid

phase. The reaction generation rate accounts for heat generated due to
ionic flux and over-potentials, and the reversible generation rate takes
into account the heat rise due to the entropy change in the electrodes’
structure. The next section uses the notation x̂0 = la , x̂ p = la + l p ,
x̂s = la + l p + ls , and x̂n = la + l p + ls + ln . For a clearer compre-
hension, bold is used in tables for coefficients whose dependence on
other variables is made explicit in other equations. The nomenclature
of the variables is reported in List of symbols section. The model
equations are from Ref. 14, where for convenience the electrolyte
potential is related to the ionic flux j(x, t) rather than to the applied
current density.33,34 The thermal model is taken from Ref. 32 while all
of the parameters describing the particular chemistry are taken from
Ref. 20.

Numerical Implementation
Most numerical methods for model-based estimation and control
algorithms require the model to be formulated in terms of AEs or
DAEs rather than PDAEs. Different numerical methods can be used
to achieve this objective. The reformulation process from PDAEs to Figure 1. Example of a 2D FVM mesh where the set of neighbor cells C(k)
AEs or DAEs is carried out by discretizing the domains of the in- is represented by the green cells.
dependent variables (e.g., the time domain t and the n-dimensional
spatial domain x ∈ Rn ). The discretization can involve both time and
space, to produce AEs, or only space, to produce DAEs. An example According to the FVM, the spatial domain  is divided into a set
of discretization in time and space is given by the FTCS (Forward- of disjoint control volumes (CVs) k centered in xk ∈ R N , such that
Time Central-Space) approach.35 Other techniques, like the method  = ∪k k and i ∩  j = ∅ , ∀i = j. The average value of the
of lines (MOL),36 discretize only the space domain and leave the time unknown variables for each CV is
as a continuous variable. When this latter approach is used, finite 
1
volume, finite difference, or finite element methods can be employed φ̄k (t) ≈ φ(x, t) d V
to obtain the set of DAEs. Alternative approaches can be used. For G k k
example, orthogonal collocation can be used with an efficient coordi- where G k represents the volume of k . Using this equation, the inte-
nate transformation to solve the set of resulting DAEs.20 In this paper, grals in 3 can be reformulated as
in order to exploit the properties of variable-step solvers, MOL is 
used to reformulate the original set of PDAEs. In particular, the finite φ̄˙ k (t) + (F(φ̄) · n)k, j ≈ s̄k (t) [4]
volume method (FVM) is employed. Due to its ability to conserve j∈C(k)
properties with high accuracy (within numerical roundoff errors), the
FVM has been used in literature to discretize models in a wide range where C(k) is the set of the neighbor cells to the kth CV and (F(φ̄)·n)k, j
of applications, such as heat transfer problems,37 flow and transport is the normal component of the numerical approximation of f (φ) · n,
in porous media,38 or more general applications for hyperbolic prob- directed toward x j starting from xk . An illustrative example of the
lems as discussed in Ref. 39. In particular, the FVM together with the set C(k) is given in Fig. 1. Suitable numerical approximations need
harmonic mean (HM) have been used to deal with possible disconti- to be employed for the term F(φ̄); given that the average values of
nuities across different sections of the cell. To the authors’ knowledge, the unknown variables φ̄ are computed in the FVM, interpolation
no published work addresses in detail the numerical issues related to techniques are employed to recover the value of such unknowns at the
the implementation of the Li-ion cell model and, in particular, the edges of the CVs.40 The approximation of F(φ̄) is discussed in the
handling of boundary conditions that ensure physical meaningfulness next section.
of the obtained solutions. For this reason, all the numerical details are
addressed below. Discretization of the governing equations.—The discretization
method introduced in Finite volume formulation section is exploited
Finite volume formulation.—Consider a general diffusion- to reformulate the set of governing equations summarized in Table I.
convection equation defined on a domain in R N of the form Given that all the unknowns of the Li-ion cell model are functions of
the variables t ∈ R+ and x ∈ R, the development of a 1D FVM model
∂φ
+ ∇(ηφ) = ∇(∇φ) + s [2] is addressed. In order to correctly carry out the discretization process,
∂t a mesh structure is defined by subdividing the spatial domain into
where φ is the unknown variable, η is the velocity,  is a diffusion Na + N p + Ns + Nn + Nz non-overlapping volumes with geometrically
coefficient and s a source term. Both the unknown φ and the source centered nodes (as depicted in Fig. 2). Every CV is associated with a
term s depend on time t and space x ∈ R N . For convenience define
f (φ) := ηφ − ∇φ. Integrating 2 over a spatial domain  ⊂ R N
and applying the divergence theorem produces the integral form of
the conservation law:
  
∂φ
dV + ( f (φ) · n) d S = sd V [3]
 ∂t d 

where d is the boundary of the domain , n is the outward pointing

unit normal on the boundary of the domain, and d V and d S repre-
sent the infinitesimal volume of  and the infinitesimal surface of
the boundary d respectively. Alternatively, this integral equation
could be written directly as an exact conservation equation over any
prescribed spatial domain. Figure 2. One-dimensional finite volume mesh.
A1196 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

Table III. FVM P2D equations.

Current Collectors, i ∈ {a, z} Boundary Conditions

 

 x 1
 k+ 2 2 (t) λi ∂ T∂(x,t)
x  = h(Tref − T̄1 (t))
ρi C p,i ∂ T̄∂tk (t) =
Iapp
(T) 1
λi ∂ T∂(x,t)  +  0
xi x σ eff,i ∂ T (x,t) 
x
k− 12 λi ∂ x  = h(T̄end (t) − Tref )
L
Positive and Negative Electrodes, i ∈ { p, n}
 x 1 
∂c (t)  k+ 2 ∂ce (x,t) 
(M1) i e,k ∂t = x
1
i
Deff,i ∂ce∂(x,t)
x  + ai (1 − t+ ) j̄k (t) ∂ x x̂ = 0
x
k− 21 
0

∂ce (x,t) 
∂x  =0
x̂n
avg
∂ c̄s (t)
(M2) ∂t = −3 Rj̄kp,i
(t)

R
c̄s∗ (t) − c̄s (t) = − Dsp,i j̄k5(t)
avg
(M3)
 x 1 eff,i  
 
(C1) σ eff,i ∂∂s (x,t)
x  2 = ai F j̄k (t) xi
k+
σ eff,i ∂s
∂x  = −Iapp
x 1 x̂0 ,x̂n
k− 2 
∂s (x,t) 
∂x  =0
x̂ p ,x̂s
 x 1  x 1 ∂e (x,t) 
∂e (x,t)  k+ 2   =0
(C2) κeff,i ∂x  − κeff,i T (x, t)ϒ ∂ ln c∂ex(x,t)  k+ 2 = xi ai F j̄k (t) ∂x x̂0
x 1

k− 2
x
x 1
k− 2 
¯ e,end = 0
∂ T (x,t)  k+ 12
(T) ρi C p,i ∂ T̄∂tk (t) = x 1
i
λ i ∂x  + Q̄ source,k
x 1
 k− 2
 
max − c̄∗ (t))c̄∗ (t) sinh 0.5R η̄ (t)
j̄k (t) = 2keff,i c̄e,k (t)(cs,i s,k s,k F T̄k (t) i,k
η̄i,k (t) = ¯ s,k (t) − ¯ e,k (t) − Ū i,k

Separator, i = s
 x
∂c (t) ∂ce (x,t)  k+ 21
(M1) i e,k∂t = x1
i
D eff,i ∂ x 
x 1
 x 1  k− 2
x 1
∂e (x,t)  k+ 2 
(C2) κeff,i ∂x  − κeff,i T (x, t)ϒ ∂ ln c∂ex(x,t)  k+ 2 = 0
x 1 x 1

k− 2
x k− 2
∂ T̄k (t) ∂ T (x,t)  k+ 12
(T) ρi C p,i ∂t = xi λi ∂ x 
1
+ Q̄ ohm,k
x
k− 12

center xk and spans the interval [xk− 1 ; xk+ 1 ]. To facilitate the treatment between two CVs. In order to recover such value, linear interpola-
2 2
of boundary and interface conditions, the edges of each CV are aligned tion techniques are used. The same approach is also applied for ce
with the domain boundaries and internal interfaces. The width of every and κeff .
CV is defined as xi = li /Ni , where i represents a particular section As discussed in Finite volume formulation section, a suitable nu-
of the battery. merical approximation for F(φ̄) is needed. Given that no convective
Once the discretization mesh is structured, the governing equa- terms are present in the set of governing equations, numerical ap-
tions are discretized as summarized in Table III. All the interface proximation is only required for the diffusive terms (e.g., −∇φ). In
conditions used to enforce continuity between adjacent materials are this work, all the diffusive terms are numerically approximated with
discussed in Implementation of Boundary and Interface conditions a first-order scheme:
section. ∂φ(x, t)  φ̄k+1 (t) − φ̄k (t)
Particular attention is required for the thermal dynamics. The re-  ≈
versible and reactive heat sources can be discretized as ∂x x 1
k+
x
2
∂Ui,k
Q̄ rev,k = Fai j̄k (t) T̄k (t) ∂φ(x, t)  φ̄k (t) − φ̄k−1 (t)
∂T  ≈
∂x x 1
k−
x
Q̄ rxn,k = Fai j̄k (t) η̄i,k (t) 2

whereas the derivatives present in the ohmic source are numerically All the values coming from the additional equations in Table II are
approximated as obtained as a function of the average values of the unknowns. Equation
(T) is used to obtain the values of T , while equations (M1), (M2), and
∂s (x, t)  ¯ s,k+1 (t) − ¯ s,k−1 (t) (M3) are used to obtain the values of ce , csavg , and cs∗ respectively.
 ≈
∂x xk 2xi
∂e (x, t)   e,k+1 (t) − 
¯ ¯ e,k−1 (t)
 ≈
∂x xk 2xi
∂ ln ce (x, t)  c̄e,k+1 (t) − c̄e,k−1 (t)
 ≈
∂x xk 2xi c̄e,k (t)
using a central differencing scheme. Finally the term Q̄ source,k :=
Q̄ ohm,k + Q̄ rev,k + Q̄ rxn,k .
Equation (C2) in Table III requires the evaluation of T (x, t),
ce (x, t), and κeff at the edges of the CVs. For example, consider Fig. 3,
where the value of the unknown T̄ has to be evaluated at the interface Figure 3. Interpolation technique to recover edge values of the unknowns.
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1197

whereas
∂ c̄e,k+1 (t) (c̄e,k+2 (t) − c̄e,k+1 (t))
s = Deff,k+ 3
∂t 2 xs2

(c̄e,k+1 (t) − c̄e,k (t))

− Deff,k+ 1
2 xs (x̃)
Figure 4. Electrolyte diffusion process: interface across the cathode and
x +x
separator. in the first volume of the separator, with x̃ = s
2
p
. The
same approach is used to enforce interface conditions where
needed.
The values of s are obtained from (C1) while the values of e are When dealing with battery packs, in particular with series-
calculated through (C2). connected cells, all the aforementioned numerical schemes have to
be replicated for each cell. Moreover, when temperature dynamics
Implementation of boundary and interface conditions.—Bound- are considered, the numerical scheme has to be adapted in order to
ary conditions must be enforced to have a physically meaningful account for continuity fluxes across the cells. Indeed, if two cells are
solution. As shown in Table I, null-flux boundary conditions on the connected in series,
electrolyte diffusion equation ce can be straightforwardly enforced by ∂ T1 (x, t)  ∂ T2 (x, t) 
imposing ∂c e
= 0 at x = x̂0 and x = x̂n . The same procedure can −λz,1  ∗ = −λa,2  ∗
∂x ∂x x=x ∂x x=x
be used to enforce ∂ ∂x
e
= 0 at x = x̂0 , while e = 0 at x = x̂n is must hold at the interface of the current collectors across the two
enforced by setting to zero the value of e at the last CV of the anode. cells (e.g., at x = x ∗ ), where Ti (x, t) is the temperature of the current
Solid-phase potential boundaries are enforced by substituting ∂ ∂x
s
at collector of the ith cell. Finally, Newton’s law of cooling has to be
x = x̂0 and x = x̂n the value of −Iapp /σeff,i . Similarly, at x = x̂ p and enforced respectively at the positive current collector of the first cell
x = x̂s , ∂
∂x
s
is replaced by the value 0. To enforce heat exchange with and at negative current collector of the second cell.
the surrounding environment, the terms ∂∂Tx evaluated at x = 0 and
x = L are substituted with the terms h(Tref − T̄1 ) and −h(T̄end − Tref )
Li-ion Simulation Battery Toolbox (LIONSIMBA)
respectively. The suffixes 1 and end refer to the first and last CV of the
entire mesh. All these conditions have been formulated also for the Different implementations of Li-ion cell simulation have been re-
FVM discretization as shown in Table III. ported in the literature written in such languages as Maple and Fortran
Due to changes in material properties along the length of the (DUALFOIL23 ), and in numerical analysis commercial software such
battery, interface conditions are required to enforce continuity of the as COMSOL Multiphysics22 and Modelica42 which provide a variety
solution. For this reason, the values of different coefficients (e.g., of models to simulate the behavior of a Li-ion cell. Matlab, how-
Deff,i , κeff,i , λi ) need to be evaluated at the interface between two ever, is the software language most commonly used by researchers
different materials. The easiest way would be to use an arithmetic for the development and evaluation of different identification, es-
mean; however, in some cases, this approach cannot accurately handle timation, and control algorithms, as Matlab has by far the largest
the abrupt changes of coefficients that may occur. Instead, the HM is number of toolboxes that implement the widest variety of such al-
employed to evaluate the value at the edges of the CVs. The HM of gorithms. Combined with its Instrument Control Toolbox that has
two generic coefficients (k1 and k2 ) can be expressed as a very extensive suite of protocols for directly communicating and
controlling laboratory equipment, Matlab has the maximum flexi-
k1 k2 bility for evaluation of control algorithms through simulations and
βk2 + (1 − β)k1 experiments.
Based on the aforementioned Li-ion cell model, this work provides
where β represents a weight to account for the difference between the a freely available Matlab based software for the simulation of Li-ion
different CV widths. A common value for β is β = xx 1
2 +x 1
, where cells, Li-ION SIMulation BAttery Toolbox (LIONSIMBA), to serve
x1 and x2 represent the CV widths. This formulation produces as a reference simulation environment for the facile development and
results that are more robust in presence of the abrupt changes of the evaluation of different ABMSs. Due to its native integration with the
coefficients, without requiring an excessively fine grid in the vicinity Matlab environment, the software facilitates the development of new
of the interface.41 algorithms, such as for the identification of Li-ion cell parameters,
Consider Fig. 4 where the interface across the last volume of the State of Charge, and optimal charging. LIONSIMBA is freely down-
cathode and the first volume of the separator is depicted. Remember loadable at: http://sisdin.unipv.it/labsisdin/lionsimba.php.
that, as discussed in Discretization of the governing equations section, The package comes with different Matlab editable .m
the mesh structure has been chosen in order to align the CV edges scripts:d
with the interfaces or physical boundaries of the battery. The value of
Deff,k+ 1 can be obtained using the HM as
2
r electrolyteDiffusionCoefficients.m: computes the electrolyte
diffusion coefficients.
Deff,k Deff,k+1 r electrolyteConductivity.m: computes the electrolyte conduc-
Deff,k+ 1 =
2 βDeff,k+1 + (1 − β)Deff,k tivity coefficients.
r openCircuitPotential.m: used to compute the Open Circuit
x
where β = x p +x
p
. The electrolyte diffusion in the last volume of Potential (OCP).
s r reactionRates.m: computes the reaction rate coefficients for
the cathode is
the ionic flux.
∂ c̄e,k (t) (c̄e,k+1 (t) − c̄e,k (t)) r solidPhaseDiffusionCoefficients.m: computes the solid phase
p = Deff,k+ 1 diffusion coefficients.
∂t 2 x p (x̃)
(c̄e,k (t) − c̄e,k−1 (t))
− Deff,k− 1
2 x 2p
d
This set of scripts refer to version 1.02 of the software; modifications or other scripts can
+ a p (1 − t+ ) j̄k (t) be added in future releases of the software.
A1198 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

All the parameters related to the simulator and to the battery are parameter values have been taken from the real battery data in Ref.
managed through the script Parameters_init.m. The customization 20, and are summarized in Table IV.
of this script allows the user to disable features such as the thermal
dynamics, change the number of CVs of the mesh, enable real-time LIONSIMBA validation.—The P2D model has been experimen-
display of results, and change the battery section lengths, thermal con- tally validated numerous times since,14 this section validates the nu-
ductivities, porosities, and so on. The user can define the operating merical implementation of LIONSIMBA by comparing the results
mode of the charge/discharge cycle by selecting between galvanos- coming from the proposed framework with COMSOL MultiPhysics
tatic, potentiostatic, or variable current profile operations. and DUALFOIL. While COMSOL has been supplied with the same
The script getInputCurrent.m contains an example for the def- model used in our framework, where a heat diffusion PDE is used to
inition of the variable current profile, and can be used to apply a describe the thermal dynamics, DUALFOIL neglects the spatial distri-
customized current profile during the simulation of the Li-ion battery. bution of the temperature and averages the heat generation rates over
A generic nonlinear function can be used for this purpose; extra pa- the cell.45 For this reason, the comparison among the three different
rameters can be used inside this function: current time instant t, initial codes is carried out considering isothermal conditions. The thermal
integration time t0 , final integration time t f , and a structure-containing model is included in the comparison with COMSOL. For isothermal
extra user data. For example, a possible implementation is and thermal enabled scenarios, a 1C discharge cycle is performed,
while the same set of parameters are maintained across the different
t − t0
I (t) = α + ξ, [α, ξ] ∈ R codes.
t f − t0 The cell potentials V (t), electrolyte concentrations ce (x, t), poten-
tials φe (x, t), and surface solid-phase concentrations cs∗ (x, t) for the
More sophisticated control strategies such as model predictive control isothermal battery are nearly identical for LIONSIMBA, COMSOL,
(MPC) can also be implemented in this framework (see the next and DUALFOIL (see Fig. 5). For the thermal enabled scenario, the
section for an example). An additional degree of freedom is set by cell potentials, temperature, and other internal states for LIONSIMBA
the possibility of defining a custom algorithm for the estimation of are nearly identical to COMSOL (see Fig. 6).
SOC and SOH. Within the Parameters_init.m script, the user can set
custom functions to be externally called after each integration step; Solid-phase diffusion models.—As introduced in Battery model
these functions will receive all the integration data of the battery and section, according to the P2D model developed in Ref. 14, diffusion
an extra structure-containing user-defined data. inside the solid particles is described using Fick’s law, where the pres-
A simulation can be initiated by calling from the Matlab command ence of a second-pseudo dimension (r ) can significantly increase the
line: computational burden. According to the particular application under
out = startSimulation(t0, tf, initialStates, I, param) study, different approximations of 1 can be employed without signif-
where icant loss of accuracy. The choice of the solid-phase diffusion model
r t0: represents the initial integration time.
should be cautious, as approximate models can have poor accuracy in
r tf: represents the final integration time.
scenarios comprising high rate of charge/discharge, short time simu-
lations, or pulse currents.26 For this reason, LIONSIMBA allows the
r initialStates: represents the structure of initial states.
r I: represents the value of the applied input current.
user to chose among three different models for solid-phase diffusion:
r param: represents the cell array of parameters structures to be r Fick’s law (including the pseudo-second dimension r ):
used in simulation.  
∂cs (r, t) 1 ∂ ∂cs (r, t)
= 2 r 2 Dseff
The structure initialStates can be used as initial state from which ∂t r ∂r ∂r
to start a simulation. If left empty, LIONSIMBA will automatically with boundary conditions
∂cs (r, t)  ∂cs (r, t) 
compute a set of consistent initial conditions (CICs) starting from
which the simulation will run. If initialStates is used as a parameter, j(x, t)
 =0  =− s
it has to be a set of CICs for the battery model in Table III. In case ∂r r =0 ∂r r =R p Deff
it is not a set of CIC, the numerical integrator will fail to converge r two-parameter polynomial approximation:25
and no results will be provided. The param array, if passed, is used
as the set of parameters for the simulation. If empty, the software ∂csavg (x, t) j(x, t)
will use a set of parameters according to the settings defined by the = −3
∂t Rp
user in the script Parameters_init.m. When designing ABMSs for
battery packs with series-connected cells, a cell-wise balancing must R p j(x, t)
cs∗ (x, t) − csavg (x, t) = −
be ensured during charging.43,44 LIONSIMBA can support the user Dseff 5
in this task by providing a full independent parametrization of each r higher-order polynomial approximation:25
cell of the series. If the param array contains a multiple parameters
structure, the software will perform a simulation of a battery pack ∂csavg (x, t) j(x, t)
composed of several cells connected in series as shown in Battery = −3
∂t Rp
pack of series-connected cells section. Each element of the pack can
be parametrized individually, leading to independent simulations of ∂q(x, t) Ds 45 j(x, t)
each cell. Finally, the out structure will contain the values of all = −30 eff q(x, t) −
∂t R 2p 2 R 2p
the dependent variables and parameters used in the simulations. The
package requires the SUNDIALS21 suite to be installed and correctly j(x, t)R p
configured with Matlab; in particular, the solver IDA is used. cs∗ (x, t) − csavg (x, t) = − + 8R p q(x, t)
35 Dseff
To obtain further help on any single script, the user can type

help <scriptname> The prediction accuracy of each approximate model is assessed by

comparison of the cell potential vs. time profiles for different C rates
from the Matlab command line or refer to the software manual. for Fick’s law with the two approximate models (see Fig. 7). The
The numerical implementation of the LIONSIMBA has been car- two-parameter approximation accurately describes the cell potential
ried out according to the rules outlined in Numerical Implementation for low to medium C rates (1C to 2C, Figs. 7a and 7b) and the
section and the cell considered is a LiCoO2 and LiC6 system. All the higher order polynomial approximation is accurate up to the 5C rate
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1199

Table IV. List of parameters used in simulation.20

Parameter Description Aluminium CC Cathode Separator Anode Carbon CC

ceinit [mol/m3 ] Initial concentration in the electrolyte – 1000 1000 1000 –
avg,init
cs [mol/m3 ] Initial solid-phase concentration – 25751 – 26128 –
csmax [mol/m3 ] Maximum solid-phase concentration – 51554 – 30555
Di [m2 /s] Electrolyte diffusivity – 7.5 × 10−10 7.5 × 10−10 7.5 × 10−10 –
Dis [m2 /s] Solid-phase diffusivity – 10−14 – 3.9 × 10−14 –
ki [m2.5 /(mol0.5 s)] Reaction rate constant – 2.334 × 10−11 – 5.031 × 10−11 –
li [m] Thickness 10−5 8 × 10−5 2.5 × 10−5 8.8 × 10−5 10−5
R p,i [m] Particle radius – 2 × 10−6 – 2 × 10−6 –
ρi [kg/m3 ] Density 2700 2500 1100 2500 8940
C p,i [J/(kg K)] Specific heat 897 700 700 700 385
λi [W/(m K)] Thermal conductivity 237 2.1 0.16 1.7 401
σi [S/m] Solid-phase conductivity 3.55 × 107 100 100 5.96 × 107
i – Porosity – 0.385 0.724 0.485 –
ai [m2 /m3 ] Particle surface area to volume – 885,000 – 723,600 –
Ds
Ea i [J/mol] Solid-phase diffusion activation energy – 5000 – 5000 –
i
E ak [J/mol] Reaction constant activation energy – 5000 – 5000 –
brugg – Bruggeman’s coefficient – 4 4 4 –
F 96485 [C/mol] Faraday’s constant – – – – –
R 8.314472 [J/(mol K)] Universal Gas constant – – – – –
t+ 0.364 Transference number – – – – –
 f,i – Filler fraction – 0.025 – 0.0326 –

(Fig. 7c), with increased error at the 10C rate typical of HEV appli-
cations (Fig. 7d).
The performance of each approximate model are quantified by
root-mean-square error (RMSE) and normalized time index (NTI) in
Table VI, where the RMSE is evaluated by comparing an approximate
model solution with respect to the full model, while the normalized
time index is the ratio between the computational time required by
an approximate model and the time required by the full model with
Fick’s law to simulate different scenarios. The two-term polynomial
approximation has much less than 1% error for the 1C and 2C rates,
with more than 1% error for higher rates. In all of the scenarios,
this approximate model takes ≈80% less time than the full model to
simulate the cell.
The higher-order polynomial approximation has a factor of 4 to
5 lower RMSE for each scenario than for the two-term model, but
an increase in computational time by a factor of two or more due to
the inclusion of another set of PDEs. Although the RMSE of 1.9%
at 10C could be small enough for some applications, the reduction in
computational time is only about 37% compared to solving the full
Fick’s law model.

Simulations
Simulation results were obtained using Matlab R2014b on a Win-
dows [email protected] PC with 8 GB of RAM for the experimental battery
parameters in Ref. 20 with a cutoff voltage of 2.5 V and environ-
mental temperature of 298.15 K. For the proposed chemistry, the 1C
value is ≈30 A/m2 . The effectiveness and ease of use of the proposed
framework are shown in a series of simulations.
In the first scenario (Fig. 8), 1C discharge simulations are compared
for a very wide range of heat exchange coefficient h, with high h being
the most challenging for retaining numerical stability in dynamic
simulations. As expected, decreasing the value of the parameter h
leads to a faster increase of the cell temperature. Moreover, due to the
coupling of all of the governing equations, it is possible to note the
influence of different temperatures on the cell voltage.
In the second scenario (Fig. 9), for a fixed value of h = 1 W/(m2 K),
different discharge cycles are compared at 0.5C, 1C, and 2C. Accord-
Figure 5. Validation of the LIONSIMBA numerical implementation in ing to the different applied currents, the temperature rises in different
isothermal conditions, with the legend given in part a. (a) Cell potential. (b) ways; it is interesting to note the high slope of the temperature during
Electrolyte Li-ion concentration. (c) Solid phase Li-ion concentration. (d) a 2C discharge, mainly due to the electrolyte concentration ce being
Electrolyte potential. driven to zero in the positive electrode by the high discharge rate.
A1200 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

Figure 7. Comparison of the three different solid-phase diffusion equations

implemented in LIONSIMBA. (a) 1C rate comparison. (b) 2C rate comparison.
(c) 5C rate comparison. (d) 10C rate comparison.

In the third scenario, the framework is used to simulate a hybrid

charge-discharge cycle, emulating the throttle of a HEV. During brak-
ing, the battery is charged. Table V resumes the configuration of the
car throttle during simulations. In Fig. 10 it is possible to analyze
the response of a single cell inside an HEV pack under a hybrid
charge-discharge cycle. In this case, the effects of temperature among
Figure 6. Validation of the LIONSIMBA numerical implementation with the different cells have been neglected. The solid potential behavior
thermal dynamics, with the legend given in part a. (a) Cell potential. (b) is primarily due to the different applied C rates, with discontinu-
Electrolyte Li-ions concentration. (c) Solid phase Li-ions concentration. (d) ous changes producing voltage drops. Different slopes of the voltage
Electrolyte potential. (e) Temperature. curve are related to the different C rates applied. Temperature rise is
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1201

Table V. Throttle configuration for hybrid charging-discharging

simulation.

Time (s) C rate Description

0–50 −1 C Moderate speed
50–60 0.5 C Charge
60–210 −0.5C Normal speed
210–410 −1C Moderate speed
410–415 −2C Overtaking
415–615 −1C Moderate speed
615–620 0.5C Charge

Figure 8. 1C discharge cycle run under different heat exchange parameters:

blue line h = 0.01W/(m2 K), dashed orange line h = 1 W/(m2 K) and dot-
dashed yellow line h = 100 W/(m2 K).

recorded in the first 50 seconds of simulations, which are followed

by a slight decrease of the temperature mainly due to the exchange of
heat with the surrounding environment (h = 1 W/(m2 K)) and due to
the lower current density applied. At around 250 s, the temperature
starts to increase due to the 1C rate applied during moderate speed;
the high slope of increase at around 410 s is due to the higher value of
the discharge current which during an overtake reaches the value of
2C. Returning to moderate speed makes the temperature slope more

Figure 9. Full discharge cycle run under different C rates: 2C (dot-dashed

yellow), 1C (dashed orange line), and 0.5C (blue line). Figure 10. Hybrid charging-discharging cycle.
A1202 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

gentle. During the last 10 seconds, temperature decreases due to the

significant change in applied current and due to dissipation of heat
with surrounding ambient. A sketch of the code used for this simula-
tion is presented in Algorithm 1.

Algorithm 1 Car cycling example code.

Input setup:
1: I = {−29.5, 14.75, −14.75, −29.5, −58, −29.5, 14.75}
 Simulation current densities
2: time = {50, 10, 150, 200, 5, 200, 10}
 Duration of each element of Iapplied (in seconds).
3: t0 = 0;  Init all the useful variables
4: tf = 0;
5: initialStates.Y = [ ];
6: initialStates.YP = [ ];
7: Phis_tot = [ ];
8: t_tot = [ ];
9: T_tot = [ ];
Core script:
10: for i = 1:length(I) do
11: tf = tf + time(i);
12: results = startSimulation(t0,tf,initialStates,I(i),[ ]);
13: Phis_tot = [Phis_tot;results.original.Phis];
 Concatenate results
14: T_tot=[T_tot;results.original.Temperature];
15: t0 = time(i);
16: initialStates = results.initialStates;
 Update initial states for the next simulation
17: end for

The simulation of the application of an ABMS is shown in Fig. 11.

In this particular simulation, a model predictive control algorithm8
is adopted to drive the State of Charge (SOC) of the battery to a
given value, while accounting for input and output constraints. The
initial SOC was around 20% and its reference value was set to 85%.
According to LIONSIMBA, the estimation of the SOC can be easily
simulated by defining a custom function. In this particular scenario,
the SOC has been computed as
 ln
1
SOC(t) = max
csavg (x, t) d x
ln cs,n 0

Algorithm 2 High level control code.

Init script:
1: t0 = 0;
2: tf = dt;  Simulations are run over a sampling time periods
3: initialStates.Y = [ ];
4: initialStates.YP = [ ];
5: Condition = 1;
Core script: Figure 11. ABMS control: an MPC algorithm8 is used to drive the charge
6: while Condition do of the battery from 20% to 85% while considering voltage, temperature, and
current constraints.
7: I = ComputeControlLaw(initialStates);
8: results = startSimulation(t0,tf,initialStates,I,[ ]);
9: [...]
Elaborate and concatenate the results The temperature maximum bound was set to 313.5 K, with the
and update the time indexes voltage set to 4.2 V. The ABMS applies a current density which is al-
10: initialStates = results.initialStates; most fixed at 1C value for the entire simulation, while starting to drop
Update initial states for the next simulation as the SOC approaches its final value. The SOC raises almost linearly
11: if SOC reached reference value then during the first 2500 s, then plateauing according to the current drop,
12: Condition = 0 in order to approach smoothly the final stage of charge. Fig. 12 shows
13: end if that, according to the different charging stages, the electrolyte concen-
14: end while tration diffuses in different ways. Starting from ceinit = 1000 mol/m3 ,
the input current induces a drop in electrolyte concentration within
the battery sections due to the diffusion of ions from the cathode to
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1203

Table VI. Comparison of different approximation methods for the diffusion in the solid particles. Root Mean Square Error (RMSE) and the
Normalized Time Index (NTI) are shown.

1C 2C 5C 10C

RMSE NTI RMSE NTI RMSE NTI RMSE NTI

two-parameter 0.082% 20% 0.25% 18% 1.6% 22% 6.6% 24%
higher-order 0.017% 37% 0.053% 44% 0.36% 60% 1.9% 63%

shows the times required by the simulator to simulate the different

scenarios, which are all under 100 s.

Battery pack of series-connected cells.—The results of a battery

pack simulation are shown in Figs. 14 and 15. To emphasize the ability
to independently parametrize each cell, in this scenario the SOC of
cell #1 is set to the 95% of its initial value while the thickness of the
cathode of cell #2 is doubled with respect to its initial value. All the
other parameters are the same for the three cells. The time responses
of the output voltage of the overall pack and the voltage of each cell
are plotted in Fig. 14. The starting voltage of the pack is around 12.1 V
and decreases subjected to a 1C discharge current. In 3346 s, the pack
is completely discharged due to cell #1 first reaching the cutoff voltage
Figure 12. ABMS control – Electrolyte concentration: The behavior of the (set to 2.5 V). The lowered starting SOC determined this behavior.
first and last volume of each section of the battery is depicted, where the The electrolyte and solid-phase surface concentrations as well as the
continuous lines belong to the cathode, the dashed lines to the separator and electrolyte potentials are compared for the three cells in Fig. 15. Cell
the dotted lines to the anode. #2 has a significantly different behavior mainly due to the presence
of a cathode with a thickness two times that of the other two cells.
This variation has effects over the output voltages, as shown in Fig.
14. Besides cell # 1 which is starting from a different SOC value, the
different behaviors of V (t) between cell #2 and cell #3 are driven by
the thickness variation.

Conclusions
This work describes a detailed procedure for the numerical imple-
mentation of the P2D model.14 Two published approximate models for
the solid-phase diffusion are also implemented, in which the pseudo-
second dimension is removed to reduce the computational complexity.
The treatment of boundary conditions is addressed with particular at-
tention to the interface conditions across the different sections of the
Figure 13. Full discharge cycle in an isothermal environment: blue line 0.5C, battery. Following the procedures and rules outlined in Numerical
dashed orange line 1C, and dot-dashed yellow line 2C. implementation section, the reader can implement his/her own ver-
sion of the model in different programming languages. Moreover, a
freely available Matlab framework LIONSIMBA is provided that is
the anode. Approaching the final stage of charging, the concentration suitable for battery design, simulation, and control. The framework
starts to converge back to the initial value of 1000 mol/m3 and, around is extended to account for different solid-phase diffusion models to
5500 s, reaches the steady value. This behavior emphasizes the prop- meet required accuracy. The simulations demonstrate high numeri-
erty of the FVM to conserve properties within numerical roundoff. cal stability for different operating scenarios. The effectiveness and
Algorithm 2 provides a high-level description of how to implement a reliability of LIONSIMBA is verified considering a heterogeneous
closed-loop controller in LIONSIMBA. sequence of applied current coming from an HEV and through the
In Fig. 13, simulations have been run disabling the thermal dynam- assessment of an ABMS strategy, in particular, the model predictive
ics leading to an isothermal environment. This particular configuration control of state of charge.
can be exploited in order to assess the influence of different constant A battery pack composed of series-connected cells can be sim-
temperatures at which the battery can operate. ulated by considering several independent cells with their own pa-
All the results of the proposed simulations can be reproduced by rameters. Due to its integration with the Matlab environment, the
running the example scripts available with LIONSIMBA. Table VII framework facilitates the development and test of different algo-
rithms such as control algorithms, identification procedures, or op-
timization of manufacturing parameters. The computational times of
LIONSIMBA are compared to DUALFOIL and COMSOL in Table
Table VII. Timing comparisons of different simulation scenarios. VIII. For each code, average simulation times are reported for re-
peated simulations of a 1C discharging cycle in isothermal conditions.
C rate h value Simulation Duration Effective Simulation Time LIONSIMBA and DUALFOIL have very similar computation times
1C 0.01 3523 s 72 s for all discretizations, with COMSOL being significaintly slower at
1C 1 3523 s 81 s lower discretizations. For the same numerical algorithm, an imple-
1C 100 3523 s 77 s mentation in a compiled language such as Fortran is inherently much
0.5C 1 7050 s 56 s faster than an interpreted language such as Matlab, indicating that
2C 1 1522 s 85 s the underlying numerical algorithm used by LIONSIMBA is more
A1204 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016)

Figure 14. Simulation of a 3-cell pack. The upper curve represents the overall voltage of the 3 series connected Li-ion cells, while the lower plots depict the
voltage of each cell in the pack. The different parametrization of each cell determines different behaviors.

Figure 15. Simulation of a 3-cell pack. The profiles of different internal states inside the three cells. Individual parametrizations lead to different behaviors.

efficient but the higher efficiency is offset by Matlab being an in- reduced by at least a factor of ten by using a multicore CPU using par-
terpreted language: the two effects approximately cancel so that the allel DAE solvers. Modern versions of Matlab have easy-to-implement
overall simulation times for LIONSIMBA and DUALFOIL were very built-in options for distributing calculations among multiple cores on
similar. a single CPU, and among multiple CPUs.
The results in this article demonstrate the promise of the proposed
framework as a reliable, efficient, and freely available Matlab-based
software for the P2D model simulation. Further developments such as
code optimization and distribution of compiled versions can only im- List of Symbols
prove the current performance. Moreover, as the proposed simulations ce (x, t) Electrolyte salt concentration [mol/m3 ]
were written in standard serial mode, the computation time could be csavg (x, t) Solid-phase average concentration [mol/m3 ]
cs∗ (x, t) Solid-phase surface concentration [mol/m3 ]
Deff Effective electrolyte diffusion coefficient [m2 /s]
Table VIII. Timing comparisons among different P2D model Dseff Effective solid-phase diffusion coefficient [m2 /s]
implementations. The number of discretized nodes has been set h Heat exchange coefficient [W/(m 2 K )]
equal for each section of the cell. Iapp (t) Applied current density [A/m2 ]
j(x, t) Ionic flux [mol/(m2 s)]
# of discrete nodes keff Effective reaction rate [m2.5 /(mol0.5 s)]
Q ohm Ohmic heat source term [W/m3 ]
10 20 30 40 50
Q rev Reversible heat source term [W/m3 ]
COMSOL 96 s 114 s 143 s 189 s 244 s Q rxn Reaction heat source term [W/m3 ]
DUALFOIL 28 s 57 s 97 s 137 s 185 s T (x, t) Temperature [K]
LIONSIMBA 28 s 69 s 105 s 134 s 223 s U ref Open Circuit Voltage [V]
 Journal of The Electrochemical Society, 163 (7) A1192-A1205 (2016) A1205

Greek 20. P. W. C. Northrop, V. Ramadesigan, S. De, and V. R. Subramanian, “Coordinate

∂U
|
∂ T T ref
Open Circuit Potential Entropic Variation [V/K] transformation, orthogonal collocation, model reformulation and simulation of
electrochemical-thermal behavior of lithium-ion battery stacks,” Journal of The Elec-
κeff Effective electrolyte conductivity [S/m] trochemical Society, 158(12), A1461 (2011).
e (x, t) Electrolyte potential [V] 21. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker,
s (x, t) Solid potential [V] and C. S. Woodward, “SUNDIALS: Suite of nonlinear and differential/algebraic
σ eff Effective solid-phase conductivity [S/m] equation solvers,” ACM Transactions on Mathematical Software (TOMS), 31(3), 363
(2005).
22. L. Cai and R. E. White, “Mathematical modeling of a lithium ion battery with thermal
effects in COMSOL Inc. Multiphysics (MP) software,” Journal of Power Sources,
References 196(14), 5985 (2011).
23. FORTRAN Programs for the Simulation of Electrochemical Systems. [Online]. Avail-
1. S. K. Dhar, S. R. Ovshinsky, P. R. Gifford, D. A. Corrigan, M. A. Fetcenko, and able: http://www.cchem.berkeley.edu/jsngrp/fortran.html.
S. Venkatesan, “Nickel/metal hydride technology for consumer and electric vehicle 24. V. R. Subramanian, V. D. Diwakar, and D. Tapriyal, “Efficient macro-micro scale
batteries – A review and update,” Journal of Power Sources, 65(1), 1 (1997). coupled modeling of batteries,” Journal of The Electrochemical Society, 152(10),
2. D. Linden, Handbook of Batteries and Fuel Cells, New York: McGraw-Hill Book A2002, (2005).
Company (1984). 25. V. Ramadesigan, V. Boovaragavan, J. C. Pirkle, and V. R. Subramanian, “Efficient
3. P. Ruetschi, “Review on the lead-acid battery science and technology,” Journal of reformulation of solid-phase diffusion in physics-based lithium-ion battery models,”
Power Sources, 2(1), 3 (1977). Journal of The Electrochemical Society, 157(7), A854 (2010).
4. W. Zhang, “A review of the electrochemical performance of alloy anodes for lithium- 26. Q. Zhang and R. E. White, “Comparison of approximate solution methods for the solid
ion batteries,” Journal of Power Sources, 196(1), 13 (2011). phase diffusion equation in a porous electrode model,” Journal of Power Sources,
5. P. Van den Bossche, F. Vergels, J. Van Mierlo, J. Matheys, and W. Van Autenboer, 165(2), 880 (2007).
“SUBAT: An assessment of sustainable battery technology,” Journal of Power 27. C. Y. Wang, W. B. Gu, and B. Y. Liaw, “Micro-macroscopic coupled modeling of bat-
Sources, 162(2), 913 (2006). teries and fuel cells: I. Model development,” Journal of the Electrochemical Society,
6. S. Santhanagopalan, Q. Guo, P. Ramadass, and R. E. White, “Review of models 145(10), 3407 (1998).
for predicting the cycling performance of lithium ion batteries,” Journal of Power 28. S. Liu, “An analytical solution to Li/Li+ insertion into a porous electrode,” Solid
Sources, 156(2), 620 (2006). State Ionics, 177(1), 53 (2006).
7. V. Ramadesigan, P. W. C. Northrop, S. De, S. Santhanagopalan, R. D. Braatz, and 29. L. O. Valøen and J. N. Reimers, “Transport properties of LiPF6 -based Li-ion battery
V. R. Subramanian, “Modeling and simulation of lithium-ion batteries from a sys- electrolytes,” Journal of The Electrochemical Society, 152(5), A882 (2005).
tems engineering perspective,” Journal of The Electrochemical Society, 159(3), R31 30. D. Bernardi, E. Pawlikowski, and J. Newman, “A general energy balance for battery
(2012). systems,” Journal of the Electrochemical Society, 132(1), 5 (1985).
8. M. Torchio, N. A. Wolff, D. M. Raimondo, L. Magni, U. Krewer, R. B. Gopaluni, 31. T. M. Bandhauer, S. Garimella, and T. F. Fuller, “A critical review of thermal is-
J. A. Paulson, and R. D. Braatz, “Real-time model predictive control for the op- sues in lithium-ion batteries,” Journal of The Electrochemical Society, 158(3), R1
timal charging of a lithium-ion battery,” in Proceedings of the American Control (2011).
Conference, 4536 (2015). 32. K. Kumaresan, G. Sikha, and R. E. White, “Thermal model for a li-ion cell,” Journal
9. V. H. Johnson, A. A. Pesaran, T. Sack, and S. America, Temperature-dependent bat- of the Electrochemical Society, 155(2), A164 (2008).
tery models for high-power lithium-ion batteries, Golden, Colorado: National Re- 33. K. Smith and C.-Y. Wang, “Solid-state diffusion limitations on pulse operation of a
newable Energy Laboratory (2001). lithium ion cell for hybrid electric vehicles,” Journal of Power Sources, 161(1), 628
10. B. Y. Liaw, G. Nagasubramanian, R. G. Jungst, and D. H. Doughty, “Modeling of (2006).
lithium ion cells – A simple equivalent-circuit model approach,” Solid State Ionics, 34. D. M. Bernardi and J.-Y. Go, “Analysis of pulse and relaxation behavior in lithium-ion
175(1), 835 (2004). batteries,” Journal of Power Sources, 196(1), 412 (2011).
11. H. He, R. Xiong, and J. Fan, “Evaluation of lithium-ion battery equivalent circuit 35. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, New York: Springer-
models for state of charge estimation by an experimental approach,” Energies, 4(4), Verlag, vol. 12, (2013).
582 (2011). 36. W. E. Schiesser, The Numerical Method of Lines, Academic Press, San Diego,
12. X. Hu, S. Li, and H. Peng, “A comparative study of equivalent circuit models for (1991).
Li-ion batteries,” Journal of Power Sources, 198, 359 (2012). 37. J. C. Chai and S. V. Patankar, “Finite-volume method for radiation heat transfer,”
13. S. K. Rahimian, S. Rayman, and R. E. White, “Comparison of single particle and Advances in Numerical Heat Transfer, 2, 109 (2000).
equivalent circuit analog models for a lithium-ion cell,” Journal of Power Sources, 38. P. Jenny, S. H. Lee, and H. A. Tchelepi, “Adaptive multiscale finite-volume method
196(20), 8450 (2011). for multiphase flow and transport in porous media,” Multiscale Modeling & Simula-
14. C. M. Doyle, “Design and Simulation of Lithium Rechargeable Batteries,” Ph.D. tion, 3(1), 50 (2005).
dissertation, University of California, Berkeley, August 1995. 39. R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge, UK:
15. R. B. Gopaluni and R. D. Braatz, “State of charge estimation in Li-ion batteries using Cambridge University Press, vol. 31 (2002).
an isothermal pseudo two-dimensional model, ” in Proceedings of the 10th IFAC 40. J. H. M. ten Thije Boonkkamp and M. J. H. Anthonissen, “The finite volume-
International Symposium on Dynamics and Control of Process Systems, 135 (2013). complete flux scheme for advection-diffusion-reaction equations,” Journal of Sci-
16. K. L. Man, K. Wan, T. O. Ting, C. Chen, T. Krilavicius, J. Chang, and S. H. Poon, entific Computing, 46(1), 47 (2011).
“Towards a hybrid approach to SOC estimation for a smart battery management 41. S. Patankar, Numerical Heat Transfer and Fluid Flow, Boca Raton, Florida: CRC
system (BMS) and battery supported cyber-physical systems (CPS),” in Proceedings Press (1980).
of the 2nd Baltic Congress on Future Internet Communications, 113 (2012). 42. M. Tiller, Introduction to Physical Modeling with Modelica, Boston: Kluwer Aca-
17. J. Remmlinger, M. Buchholz, M. Meiler, P. Bernreuter, and K. Dietmayer, “State- demic Publishers, (2012).
of-health monitoring of lithium-ion batteries in electric vehicles by on-board internal 43. S. W. Moore and P. J. Schneider, “A review of cell equalization methods for lithium
resistance estimation,” Journal of Power Sources, 196(12), 5357 (2011). ion and lithium polymer battery systems,” Society of Automative Engineers, Tech.
18. Z. Guo, X. Qiu, G. Hou, B. Y. Liaw, and C. Zhang, “State of health estimation for Rep., 2001.
lithium ion batteries based on charging curves,” Journal of Power Sources, 249, 457 44. W. F. Bentley, “Cell balancing considerations for lithium-ion battery systems, ” in
(2014). Proceedings of the Twelfth Annual Battery Conference on Applications and Advances,
19. Y.-H. Chiang, W.-Y. Sean, and J.-C. Ke, “Online estimation of internal resistance and 223 (1997).
open-circuit voltage of lithium-ion batteries in electric vehicles,” Journal of Power 45. L. Rao and J. Newman, “Heat-generation rate and general energy balance for insertion
Sources, 196(8), 3921 (2011). battery systems,” Journal of the Electrochemical Society, 144(8), 2697 (1997).