- Email: [email protected]

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Bayesian calibration for electrochemical thermal model of lithium-ion cells Piyush Tagade a, Krishnan S. Hariharan a, *, Suman Basu a, Mohan Kumar Singh Verma a, Subramanya Mayya Kolake a, Taewon Song b, Dukjin Oh b, Taejung Yeo b, Seokgwang Doo b a Computational Simulations Group (SAIT-India), Samsung R&D Institute India-Bangalore, #2870 Phoenix Building, Bagmane Constellation Business Park, Outer Ring Road, Doddanekundi Circle, Marathahalli Post, Bangalore, 560 037, India b Energy Material Lab, SAIT, Samsung Electronics, Gyeonggi-do, 443-803, Republic of Korea

h i g h l i g h t s Bayesian framework for calibration of the P2D-ECT model. Matrix-variate Gaussian process for computationally efﬁcient implementation. P2D-ECT model parameter estimation and quantiﬁcation of model uncertainty. Accurate model prediction across a range of temperatures. Novel insights into low temperature Li-ion cell operation.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 January 2016 Received in revised form 18 April 2016 Accepted 22 April 2016 Available online 29 April 2016

Pseudo-two dimensional electrochemical thermal (P2D-ECT) model contains many parameters that are difﬁcult to evaluate experimentally. Estimation of these model parameters is challenging due to computational cost and the transient model. Due to lack of complete physical understanding, this issue gets aggravated at extreme conditions like low temperature (LT) operations. This paper presents a Bayesian calibration framework for estimation of the P2D-ECT model parameters. The framework uses a matrix variate Gaussian process representation to obtain a computationally tractable formulation for calibration of the transient model. Performance of the framework is investigated for calibration of the P2D-ECT model across a range of temperatures (333 Ke263 K) and operating protocols. In the absence of complete physical understanding, the framework also quantiﬁes structural uncertainty in the calibrated model. This information is used by the framework to test validity of the new physical phenomena before incorporation in the model. This capability is demonstrated by introducing temperature dependence on Bruggeman's coefﬁcient and lithium plating formation at LT. With the incorporation of new physics, the calibrated P2D-ECT model accurately predicts the cell voltage with high conﬁdence. The accurate predictions are used to obtain new insights into the low temperature lithium ion cell behavior. © 2016 Elsevier B.V. All rights reserved.

Keywords: Lithium-ion batteries Electrochemical thermal model Bayesian calibration Parameter estimation Markov chain monte carlo

1. Introduction Majority of the contemporary consumer electronic devices use rechargeable batteries as the primary energy source, while many other industrial applications employ rechargeable batteries as a secondary source of energy. Rechargeable batteries have also found

* Corresponding author. E-mail address: [email protected] (K.S. Hariharan). http://dx.doi.org/10.1016/j.jpowsour.2016.04.106 0378-7753/© 2016 Elsevier B.V. All rights reserved.

a renewed interest in the automotive industry. Continuous improvement of the battery technology is essential for advancements of the electronic devices and other applications. Advancements in the battery technology are particularly important for the automotive industry, where novel electro-chemistries and cell designs can enable the development of cleaner electric (EVs) and hybrid-electric vehicles (HEVs) [1,2]. Accurate and credible physicsbased battery models are required for investigation of these new electro-chemistries and subsequent optimized cell designs. The

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

physics-based battery models (often after model order reduction [3]) can also be used on-board the next generation battery management systems for intelligent cell monitoring and control purposes. Reliable cell design and subsequent development of the model based battery management system require accurate model predictions. The designer or an on-board model based controller is expected to make safety critical decisions conditional on the model credibility [4]. Here, the model credibility reﬂects the user conﬁdence on the model predictions [5]. Thus, establishing accuracy and the credibility assessment [6,7] of the physics-based battery models is a critical prerequisite for the future advancements of the battery technologies, including the development of novel cell designs and the battery management systems. Lithium ion batteries are currently the most widely used rechargeable energy source for consumer electronic devices. Due to their low weight, low self-discharge rate, high energy and high power density, lithium-ion batteries are particularly attractive for automotive applications [8,9]. With the ever-increasing applications, mathematical modeling of the lithium ion batteries has attracted a signiﬁcant interest of research community. Different lithium ion battery models, varying from the low ﬁdelity empirical equivalent circuit models to the high ﬁdelity continuum scale electrochemical models are reported in the literature [10,11]. Several authors use equivalent circuit models for lithium ion battery design, optimization and control tasks [12e15]. Comparison of the most widely used equivalent circuit models is available in the literature [16]. Though computationally efﬁcient, the equivalent circuit models do not have a capacity to capture internal electrochemistry of the battery, and have limited accuracy and predictive capabilities [10]. For high energy automotive applications like EVs and HEVs, internal electro-chemical information can be exploited to extend the driving range while ensuring a safe battery operation. Continuum scale pseudo-two dimensional electrochemical thermal models (P2D-ECT) provide this detailed internal battery information. The P2D-ECT models are based on the porous electrode theory developed by Newman et al. [17,18]. Governing equations for the P2D-ECT model consists of a set of ﬁve coupled PDEs; four PDEs that deﬁne charge and mass balance for solid and electrolyte phases, and a ﬁfth PDE that models thermal balance for the battery [19,20]. Butler-Volmer kinetics is used to obtain closure for the governing equations [10,21]. In the P2D-ECT model, the solid and electrolyte phase charge balance, electrolyte phase mass balance and energy balance equations are solved along the length of the cell. The solid phase mass balance equation is solved for a spherical particle in a so called pseudo-second dimension [10]. Details of the P2D-ECT model can be found in Refs. [10,11] and its reformulations are discussed in Ref. [22] and references therein. The P2D-ECT model contains many parameters that are unknown or poorly known. Moreover, some of the electrochemical behavior of the lithium ion battery, especially at low temperature and high current operations, is still poorly understood [23]. This poorly understood electro-chemistry either remains unaccounted, or heavily approximated, in the P2D-ECT model. The poorly known parameters and the unaccounted/approximated electro-chemistry induce uncertainty in the P2D-ECT model outputs. A detailed model calibration process is essential to reduce these uncertainties. The model calibration process is expected to achieve two key objectives [4]: a) parameter estimation to reduce uncertainty in the poorly known parameters; and b) quantiﬁcation of the model uncertainty due to unaccounted/approximated electro-chemistry. The quantiﬁed model uncertainty assesses the credibility of the model. Current literature on calibration of lithium ion battery models, however, is limited to the parameter estimation. In the context of lithium ion batteries, the parameter estimation

297

methods are primarily focused on the equivalent-circuit models. Several authors use traditional least-square regression methods for equivalent circuit model parameter estimation [24e26]. Although limited, there are notable exceptions where parameter estimation of full P2D-ECT model is reported [27]. Forman et al. [28] use genetic algorithm (GA) for estimation of the P2D-ECT model parameters. Reimers et al. [29] use Levenberg-Marquardt (LM) algorithm in a least-square setting to estimate parameters that minimize the squared error between the predicted and measured cell voltage across a range of temperature and load current. The parameter estimation using methods like GA or MA, however, require several executions of the P2D-ECT model. The high computational cost of the P2D-ECT model, thus, renders these methods computationally very expensive. Several other authors use model order reduction to ensure computational tractability of the parameter estimation methods. Ramadesigan et al. [30,31] use the Gauss-Newton method with a reformulated model to estimate the parameters. Santhagopalan et al. [32] use the LM algorithm with a single particle model for parameter estimation. Masoudi et al. [27] use the homotopy optimization method with their reduced order model for parameter estimation that ensures convergence of the solution to the global optimum. The optimization methods for parameter estimation are limited by two key deﬁciencies: 1) the optimization methods can not quantify uncertainty in the estimated parameters; and 2) the optimization methods can not take into account any prior information about the parameters, though, a physical interpretation can be obtained for the estimated parameters a-posteriory, as demonstrated in Ref. [33]. The Bayesian inference method effectively resolves both these deﬁciencies. Effectiveness of the Bayesian inference for lithium ion battery state estimation is demonstrated by Saha et al. [34]. In the context of P2D-ECT model, Ramadesigan et al. [30,31] compare the Bayesian inference method with the traditional Gauss-Newton optimizer for parameter estimation. Although Ramadesigan et al. [30,31] quantify model parametric uncertainty, they do not consider the effect of model structural uncertainty to asses credibility of the calibrated model. In this paper, a Bayesian framework [4] is presented for calibration of the P2D-ECT model. The framework can simultaneously estimate the parameters and also assess the credibility of the model. The framework is based on a Bayesian inference method developed by Kennedy and O'Hagan [35] for calibration of a computer simulator. The Bayesian framework of Kennedy and O'Hagan [35] is widely used in the varied ﬁelds for calibration of simulators [36,37]. The implementation of the framework for calibration of the P2D-ECT model, however, is challenging due to the high computational cost and transient nature of the model output. In this paper, an extension of the Bayesian framework of Kennedy and O'Hagan [35] is proposed, which efﬁciently resolves both these challenges. Key innovation of the proposed Bayesian framework is a use of matrix variate Gaussian distribution [38e40], which makes calibration of the P2D-ECT model computationally tractable. The Markov Chain Monte Carlo (MCMC) sampling [41,42] is used for numerical implementation of the proposed framework. Research work presented in this paper advances the current state of the art as follow: a) to the authors knowledge, this research work is a ﬁrst reported full calibration of the spatially resolved P2DECT model; b) use of the matrix-variate Gaussian process priors is novel. The proposed algorithm for numerical implementation of the framework extends current state of the art of the Bayesian calibration methods; c) this paper introduces new physics to the existing P2D-ECT model and a method for efﬁciently testing the validity of the new physics. The paper also present new insights into the temperature dependence of the Bruggeman's relation. As compared to the state of the art optimization based P2D-ECT

298

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

model parameter estimation methods, the proposed approach provide following additional information: a) a formal framework to investigate validity of the new physics before incorporation into the model. b) probability distribution of the model uncertainty. Variance of this model uncertainty quantiﬁes credibility of the model. c) updated probability distribution of the estimated parameters. The relative update of the probability distribution quantiﬁes sensitivity of the parameters. These functionalities are demonstrated using numerical results presented in this paper. Though the effectiveness of the method is demonstrated in this paper for a particular Lithium typology, the method is generic in nature and can be applied to different typologies without any change. Details of the P2D-ECT model are provided in section 2. Section 3 brieﬂy describes the proposed Bayesian framework for calibration of the P2D-ECT model. Detailed formulation of the framework is provided in the supplementary details. Effectiveness of the framework is demonstrated using numerical results presented in Section 4. The section also presents physical insights into the low temperature behavior of the lithium ion cell and temperature dependent Bruggeman's relation. The paper is summarized and concluded in Section 5. 2. Pseudo-two dimensional electrochemical thermal model The P2D-ECT model used in this paper is based on the porouselectrode theory of Newman et al. [17,18]. This model couples electrochemical description of the lithium transport, reaction kinetics and the thermodynamics at the electrode level with the thermal model that describes cell level temperature transients. The lithium ion cell consists of two porous electrodes that are separated by an ionically conducting but electronically insulating separator. The separator and the electrodes are immersed in an appropriately conducting electrolyte. The porous electrodes are composed of active solid materials that store intercalated lithium. During discharge, lithium deintercalates from an anode active material, releasing an electron in the outer circuit and a lithium ion into the electrolyte. These lithium ions diffuse through the separator towards the cathode region. To ensure electro-neutrality, the electron instantaneously reacts with a lithium ion at the cathodeelectrolyte interface to form a lithium molecule that intercalates into the cathode active material. The intercalated lithium subsequently diffuses inside the active material lattice. The similar transport takes place from cathode to anode during the charging. This behavior is macroscopically modeled using the P2D-ECT model. The schematic of the lithium ion cell and its approximate representation used in the P2D-ECT model is described in Fig. 1. The P2D-ECT model divides the lithium-ion cell into the anode, separator and the cathode region. The one dimensional celldynamics is described schematically along the through-plane direction (x-direction). Along the length of the electrodes, active material is approximated by spherical particles located at each discretized node. The evolution of the solid phase potential, fs, at the particle surface is given by (all the symbols are summarized in the nomenclature),

v vx

seff

vfs vx

¼ as FjLi ;

(1)

Fig. 1. Schematic of the Lithium-ion cell.

Meff ¼ MεBRUG ;

(2)

where ε is porosity and BRUG is known as the Bruggeman's coefﬁcient. In the P2D-ECT model, different Bruggeman's coefﬁcient values are used for the solid and the electrolyte phases. Note that other forms of the Bruggman relation are also used in the literature (see Ref. [43] and references therein), however, those forms are not explored in this paper. The electrolyte phase potential, fe, is similarly governed by the electrolyte phase charge balance equation

v vf vlogce ¼ as FjLi : keff e þ keff D vx vx vx

(3)

The transport of the lithium ions in the electrolyte is given by the diffusion equation

εe

vce v eff vce De ¼ þ as ð1 tþ ÞjLi : vx vt vx

(4)

The charge balance Equations (1)e(3) and the diffusion Equation (4) are subjected to the appropriate Neumann boundary conditions at the current collector ends. The diffusion of the lithium inside the active particle is modeled using the spherical diffusion equation

vcs 1 v vcs r 2 Ds ¼ 2 ; vt vr r vr

(5)

eff

where s is the effective conductivity. The P2D-ECT model uses effective transport properties to account for the tortuous path of the porous material. An empirical Bruggeman's relation is used to obtain these effective transport properties [10]. For a given transport property M, the Bruggeman's relation is given as

with the Neumann boundary conditions

vcs vcs ¼ 0; Ds ¼ jLi : vr r¼0 vr r¼Rs

(6)

The closure for the governing equations is obtained using the

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

Butler-Volmer reaction kinetics [10].

jLi ¼ i0

aa F aF exp h exp c h ; RT RT

(7)

where the exchange current density, i0, is given by

a i0 ¼ k cs;max cs a ðcs Þac ðce Þaa :

(9)

where U is the open circuit potential and the Rf is the SEI ﬁlm resistance. The electrochemical model is coupled with the thermal model using the energy balance equation

rCp

vT ¼ q_ gen þ q_ conv : vt

(10)

The heat generation rate consists of heat generation due to reaction, reversible heat generation due to the entropy changes and the ohmic heat generation due to the lithium ion transport. Refer [44] and the references therein for a detailed description of the heat generation terms. The temperature obtained from (10) is used in the Arrhenius's relation to model the temperature dependence of the transport properties. The output of the P2D-ECT model is the terminal voltage of the lithium-ion cell for a given charge-discharge cycle. The terminal cell voltage, Vcell, is given by the solid phase potential difference between cathode and anode current collector,

Vcell ¼ fs;c fs;a :

developed in this paper for calibration of the P2D-ECT model. The Bayesian calibration method is ﬁrst introduced for a generic computer simulator with a scalar output, which is then extended for calibration of the transient simulators. The detailed formulation and the numerical implementation of the algorithm is provided in the Supplementary material.

(8)

The Butler-Volmer reaction (7) is driven by the over-potential h given by

h ¼ fs fe U Rf jLi ;

299

(11)

The governing equations and the associated boundary conditions for the P2D-ECT model are summarized in the Table 1. The governing equations are numerically implemented using a fully implicit ﬁnite volume method. 3. Bayesian calibration of the P2D-ECT model This section brieﬂy describes the Bayesian formulation

3.1. Bayesian calibration of computer simulators The Bayesian calibration method is introduced for a simulator, T(x,q), which accepts control inputs x and parameters q as inputs and predicts a scalar output y. The parameters are often unknown/ poorly-known, inducing uncertainty into the simulator predictions. Let qt denote the true but unknown value of the parameters. The model is structurally uncertain due to the approximated or unaccounted physics. For example, tortuous path of the porous electrode is approximated in the P2D-ECT model using the Bruggeman's relation. Due to the structural uncertainty, the simulator prediction deviates from the true system response even if the true parameters, qt, are completely known. This deviation is modeled using a discrepancy function d(x). Thus, the relation between the true system response, x(x), and the simulator prediction is given by

xðxÞ ¼ Tðx; qt Þ þ dðxÞ:

(12)

Similarly, relation between the system response and the experimental observations is given by

ye ðxÞ ¼ xðxÞ þ ε;

(13)

where ε is the measurement uncertainty. The parametric, model structural and measurement uncertainties are quantiﬁed using appropriate probability distributions. The Bayesian framework uses the Bayes theorem to update the probability distributions of the parameters and the discrepancy function by assimilating the experimental observations. For Eqs.12 and 13, the Bayes theorem gives [35].

pðqt ; dðxÞjye ðxÞÞfpðye ðxÞjqt ; dðxÞÞpðqt ; dðxÞÞ:

Here p(qt,d(x)) is known as the prior, p(ye(x)jqt,d(x)) is the likelihood, and p(qt,d(x)jye(x)) is known as the posterior probability distribution. The rest of the formulation presented in this paper

Table 1 Compilation of governing equations for the pseudo-2 dimensional electrochemical thermal model. Governing equations Solid phase charge balance s seff vf ¼ as FjLi vx

v vx

Boundary conditions s seff vf vx x¼0;L ¼ I s seff vf vx x¼La ;La þLs ¼ 0

Electrolyte phase charge balance v eff vfe þ keff v log ce ¼ as FjLi vx k vx vx D Electrolyte phase mass balance eff vce v e þ as ð1 tþ ÞjLi εe vc vt ¼ vx De vx Solid phase mass balance v r 2 D vcs ¼ r12 vr s vr

vcs vt

e keff vf vx

x¼0;L

¼0

eff vce vx x¼0;L

De

vcs vr r¼0

¼0

¼0

Li s Ds vc vr r¼R ¼ j s

Butler-Volmer Reaction Kinetics aF h exp aRTc F h jLi ¼ i0 exp aRT Thermal Balance _ _ rCp vT vt ¼ qgen þ qconv

(14)

l vT vt

x¼0;L

¼ hðT Tamb Þ

300

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

primarily focuses on an appropriate speciﬁcation of the prior and the likelihood to obtain a numerically and computationally tractable form of the posterior distribution.

3.2. Bayesian calibration of a transient model Two key challenges in the calibration of the P2D-ECT model, high computational cost and the transient output, are efﬁciently resolved in the proposed formulation by using a matrix variate Gaussian process (MVGP) prior. Similar to the scalar and the multivariate Gaussian process, the MVGP is described using a matrix variate Gaussian distribution. For an arbitrary matrix Q 2ℛ nq the matrix variate Gaussian distribution, M N n;q ðm; S; VÞ is given by

n=2 q=2 1 pðQ ÞfS V exp tr S1 ðQ mÞT V 1 ðQ mÞ ; 2

m2ℛ nq

where is a mean, S is a row-covariance matrix and V is a column-covariance matrix of the distribution. Readers are referred to Gupta and Varga [45] for more details of the matrix variate Gaussian distribution. This distribution is related to the multivariate Gaussian distribution by n;q ðm; S; VÞ

¼N

nq ðvecðmÞ; S5VÞ

(17)

where hð,Þ2ℛ m is a vector of known regression functions and b2ℛ mq is a matrix of regression coefﬁcients. A set of regression functions deﬁned for all the design points is given by H 1 ¼ ½hT ðxi ; qi Þ; i ¼ 1; …:; N. m2 ¼ 0 is used for the discrepancy function, which speciﬁes unbiased structural uncertainty in the model. A second set of regression functions at experimental control input settings is deﬁned as H 2 ¼ ½hT ðxj ; qt Þ; j ¼ 1; …:; n. Squared exponential functions are used to model row and column covariance matrices. For example, (i,j)th element of the row and covariance matrices of D1 are given by

T

V 1 ði; jÞ ¼ exp ðxi ; qi Þ xj ; qj L1 ðxi ; qi Þ xj ; qj

(18) 2

S S1 ði; jÞ ¼ s21 exp l1 ti tj ;

(15)

MN

m1 ð,Þ ¼ hT ð,Þb;

(16)

where vec(m) is a vectorization of a matrix m and 5 denotes a Kronecker product [45]. Deﬁnitions of these operators are summarized in the Appendix. The formulation presented in this paper extensively uses the equivalence between the matrix variate and the multi-variate Gaussian distributions. The MVGP representation in the present formulation is facilitated by treating the transient simulator response as a multivariate output. Thus for a given control input and parameters, the transient response is given by y ¼ T(x,q), where y¼[tj;j ¼ 1,…,q] is a simulator output at discrete time instances tj. The formulation uses a design of experiment (DOE) [46] method to select a statistically optimized set of N design points1 (xi,qi). In this paper, Latin square sampling [47] is used for the DOE. Although not reported in this paper, authors compared Latin square sampling with the Latin hypercube sampling, where the Latin square sampling gave superior performance. Let d1 ¼ [(xi,qi);i ¼ 1,….,N] denote a set of DOE selected design points. The simulator is executed at the design points d1 to obtain a set of simulator outputs, D1 ¼ ½yðxi ; qi Þ; i ¼ 1; …:; N2ℛ Nq . In the MVGP representation, D1 is considered as a sample from a matrixvariate Gaussian distribution D1 M N N;q ðm1 ; S1 ; V1 Þ. Let experimental observations are obtained at the control input settings d2 ¼ [(xj,qt);j ¼ 1,…,n] and D2 ¼ ½ye ðxj ; qt Þ2ℛ nq similarly deﬁne a set of experimental observations. Note here that the control input settings are concatenated with the true, but unknown, parameters qt. The uncertainty in the experimental observations is modeled using independent zero-mean Gaussian distributions with standard deviations se, thus D2 M N n;q ðye ; se I q ; I n Þ. The discrepancy function is similarly deﬁned at d2, thus, a set Dd ¼ ½dðxj ; qt Þ2ℛ nq is obtained such that Dd M N n;q ðm2 ; S2 ; V2 Þ. The MVGP representation is completed by specifying the modeling choices for the mean and the covariance functions. In the present formulation, a linear regression function is used for m1. Thus

(19)

where L1 is a diagonal matrix of row-correlation length parameters and lS 1 is a column-correlation length parameter. The parameters S fb; L1 ; l1 ; s21 g are the hyper-parameters (i.e. parameters of the probability distribution) of M N N;q ðm1 ; S1 ; V1 Þ. The distribution Dd M N n;q ðmd ; Sd ; Vd Þ is also similarly deﬁned using the hyper2 parameters fL2 ; lS 2 ; s2 g. For notational convenience, a combined matrix of simulation outputs and experimental data is deﬁned as D ¼ [D1;D2]. Other relevant matrices are similarly deﬁned, for e.g., H ¼ [H1;H2]. A set of covariance hyper-parameters is deﬁned as S 2 2 f ¼ ½L1 ; lS 1 ; s1 ; L2 ; l2 ; s2 . Using these deﬁnitions, the probability distribution of D for a given set of true parameters qt, hyperparameters f and the regression coefﬁcients b is

1=2 1 vecðD1 H 1 bÞ T 1 pðDjqt ; f; bÞfV exp V 2 vecðD2 H 2 bÞ

vecðD1 H 1 bÞ vecðD2 H 2 bÞ

(20)

:

(21)

Here the covariance matrix V has a following form

V¼

S1 5V 1 ðd1 ;d2 Þ S1 5V 1 ðd1 ;d1 Þ ; 2 S1 5V 1 ðd2 ;d1 Þ S1 5V 1 ðd2 ;d2 Þþ se I q 5I n þS2 5V 2 ðd2 ;d2 Þ (22)

where V1(d1,d2) is a cross-covariance between the simulation and experimental observations. Equation (21) is used in the Bayes theorem (14) with an appropriate speciﬁcation of prior probabilities to obtain the posterior probability distribution. The noninformative prior, p(b)f1, is used for the regression coefﬁcients b. In the ﬁnal step, the regression coefﬁcients b are marginalized to obtain a computationally tractable form of the Bayes theorem,

1=2 1=2 pðqt ; fjDÞfpðqt ÞpðfÞV W 0 0 0

1T

11 b b vec D1 H 1 b vec D1 H 1 b B 1B C 1 B CC

A V @

AA: [email protected] @ 2 vec D H b b b vec D2 H 2 b 2 2 (23)

1

Here and in the following, a design point refer to a set of control input and parameter values.

b ¼ W½I 5HT V 1 vecðDÞ and W ¼ ½½I 5HT V 1 ½I 5H. Here b q q q

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

The formulation (23) is completed by specifying the prior distributions for the parameters, p(qt), and the hyper-parameters p(f). In this paper, independent Gaussian distributions with known mean and variance are used to specify p(qt). The hyper-parameters of M N N;q ðm1 ; S1 ; V1 Þ have negligible effect on the estimated parameters [35], thus, the hyper-parameters s1 and l1 are ﬁxed apriory. Gamma prior is used for l2 and inverse-Gamma prior [48] is used for s22. Note that although speciﬁc choices for the prior distribution are used in this paper, the formulation works without any change for other priors. For all the test cases presented in this paper, load current and the ambient temperature are used as the control inputs, while the cell voltage is considered as the system response. Thus, each row of the matrix D contains a cell voltage transient for a given C-rate, ambient temperature and parameter values, while each column contains cell voltage for a given time instance. The output of the calibration process is the posterior probability distribution p(qt,fjD) given by Eq. (23). The MCMC sampling method is used to sample from this posterior probability distribution. See Ref. [42] and references therein for a detailed description of the MCMC algorithm. The estimated values of the parameters, qbt , are given by

qbt ¼ argmax qt

Z pðqt ; fjDÞdf;

(24)

which is a maximum a-posteriory (MAP) estimate of the parameters. Note that this MAP estimate is easily obtained from the MCMC samples. The posterior distribution of s22 assesses the credibility of the calibrated simulator. As demonstrated in this paper, this posterior distribution can also be used for testing validity of the new physics before inclusion in the model. The Bayesian formulation for calibration of the P2D-ECT model is summarized in Fig. 2. Readers are referred to the supplementary material for details of the formulation, and its numerically efﬁcient implementation.

301

4. Results and discussion Key result of this paper is an estimated set of parameters for the P2D-ECT model and the quantiﬁed model structural uncertainty. The framework is applied for estimation of the Li-ion cell parameters. The cell is completely speciﬁed using a set of parameters consisting of cell geometry deﬁnitions, electrochemical properties and the transport properties. In this paper, a subset of the most critical parameters is inferred using the Bayesian framework. This subset can be determined using sensitivity analysis, prior knowledge and the information obtained from the experts. Although demonstrated for particular parameters, the method can be used for estimation of other parameters without any change. 4.1. Methodology validation In the ﬁrst set of results, accuracy of the Bayesian framework is demonstrated using a synthetic data set. This data set is generated by designing a synthetic test bed, which is deﬁned using the P2DECT model with the completely known parameters and control input settings. This P2D-ECT model is run for room temperature (298 K) with various load currents from a fully charged state to obtain a set of synthetic cell voltage data. This cell voltage data is used for the model calibration. 0.5C discharge data at the room temperature is used for this test case. For the test case, anode and cathode solid phase VOFs and Bruggeman coefﬁcients are assumed unknown. Total 49 design points are generated using the Latin square method. The P2D-ECT model is executed at these design points and the resultant set D1 is used in the Bayesian framework for calibration. The design space, true, initial and the estimated parameter values are summarized in the Table 3. Fig. 3 show results of the calibration process. In Fig. 3a), the true cell voltage is compared against the model predictions obtained using the initial and estimated parameters. From the ﬁgure, it can

Fig. 2. Flowchart of the Bayesian framework.

302

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

Table 2 List of the parameters used in the model. S.no

Parameter

Cathode

Anode

Separator

1 2 3 4 5 6

Length (m) Particle radius (m) No. of particles Max. conc. (mol m3) React. rate const. (ms 1(mol m3)0.5) Sp. heat (J Kg1 K1)

68 106 1.42 106 17 49,495 5.7 108 1 103

92 106 16.9 106 23 28,200 7.77 109 8 102

16 106

Table 3 Estimated parameters for methodology validation. Parameter

Anode VOF Cathode VOF Anode Brugg. Cathode Brugg.

Design space Low

High

0.56 0.56 1.2 1.2

0.84 0.84 4.5 4.5

True

Initial

Calibrated

0.678 0.7 3.0 3.0

0.63 0.63 1.5 1.5

0.6816 0.7061 2.9240 2.9450

4 1150 8 102

be noted that the model prediction with the initial parameters deviates signiﬁcantly from the true response. The Bayesian framework accurately estimates the parameters, and the resultant cell voltage prediction matches closely with the true system response. For brevity, the comparison is shown only for 0.5C discharge; however, the results are comparable at other C-rates. Fig. 3b) shows the movement of the Markov Chain in the space of electrode VOFs and the cathode Bruggeman's coefﬁcient for ﬁrst thousand samples inside the design space. The ﬁgure also shows the true parameter values. The chain is initialized from a random starting point, and moves rapidly towards the true parameter values. Although the chain probe states over a complete design space, states near the true parameters are accepted with higher probability. States

(a)

(b)

(c)

(d)

Fig. 3. Validation of the proposed methodology. Frame (a) compares the true cell voltage with the predictions obtained using the initial and estimated parameters. Frame (b) shows the movement of the Markov chain. Frame (c) compares the prior and posterior probability distributions of the parameters and frame (d) compares the prior and posterior probability distributions of s22 .

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

accepted after the burnout period are samples from the posterior probability distribution. Fig. 3c) shows the prior and the posterior probability distributions (PDF) of the parameters. The distributions are scaled using the true value for better visualization. The MAP estimates of the posterior distributions are close to the true values (indicated in the ﬁgure using vertical lines). The posterior PDFs are narrow compared to the prior, signifying the increased conﬁdence on the estimated parameters. The change in the PDF of VOFs is more pronounced as compared to the Bruggeman's coefﬁcients. The narrower PDF indicates that the cell voltage contains more information about the VOFs. This implies that the cell voltage prediction is more sensitive to VOFs as compared to the Bruggeman's coefﬁcient. The posterior distribution of the Bruggeman's coefﬁcients is non-Gaussian with the mode near the true values (Fig. 3(c) bottom panel). As can be observed from the ﬁgure, the MAP is a more accurate estimator compared to the mean for such non-Gaussian distributions. The prior and posterior distribution of s22 for cell voltage and the cathode potential is shown in Fig. 3(d). The prior PDF speciﬁes nonzero probability to a wide range of s22 values. The higher value of s22 indicates wider range of probable values for the discrepancy function d(x). This wider range speciﬁes lower credibility to the P2D-ECT model. As the cell voltage data is assimilated, the posterior PDF of s22 moves toward left of the prior. This posterior PDF gives high probability to the lower values of s22 , indicating increased credibility of the model.

4.2. Computational cost of the framework Detailed computational complexity analysis of the framework is presented in the Supplementary material. This subsection brieﬂy investigates the computational cost for implementation of the framework. The computational cost analysis is performed using MATLAB R2013a on a desktop computer with Intel Core i7 3.40 GHz processor and 4.00 GB RAM. Total computational cost of the framework is divided into ofﬂine and online components. The ofﬂine computational cost is primarily incurred by N executions of the P2D-ECT model to generate the D1 matrix. Note that this computational cost depend on the complexity of the P2D-ECT model, and incurred only once for a given calibration. Table 4 summarizes the online computational cost incurred by the proposed framework. The computational cost is presented for three test cases. For all the test cases, single data point is assimilated for calibration. In the ﬁrst test case, effect of the number of simulation runs, N, on the computational time is investigated. For this test case, cell voltage data at q ¼ 10 time instances are assimilated and total 20,000 samples are collected using the MCMC method. The second test case investigates the effect of number of data time instances on the computational cost. Total N ¼ 75

Table 4 Computational cost of the Bayesian framework. Computational cost analysis Simulation data points

Time instances

MCMC samples

No.

Time (sec)

No.

Time (sec)

No.

Time (sec)

25 50 75 100 121

94.64 274.56 564.82 923.61 1266.42

5 7 10 15 20

168.04 296.88 564.82 1459.47 2849.04

5000 10,000 20000 35,000 50000

162.88 313.06 564.82 1017.28 1467.05

303

simulation data points are used for this test case and 20000 samples are collected. The computational cost increases quadratically with the increase in number of simulation data points and the number of data time instances used for the calibration. Effect of the number of MCMC samples on the computational cost is investigated in the third test case using N ¼ 75 and q ¼ 10. Computational cost of the proposed framework increases linearly with the total number of MCMC samples. Around 10 data time instances sufﬁces for practical implementation, where inference can be obtained within half an hour of the computational time. In comparison, the P2D-ECT model parameter estimation using the GA algorithm reported in Ref. [28] takes about three weeks on a ﬁve quad core processor. 4.3. Experimental data The results presented in the previous subsections have demonstrated accuracy of the Bayesian framework to efﬁciently calibrate the P2D-ECT model using the synthetic test bed data. In this and the subsequent subsections, experimental data obtained for a Li-ion cell is used to calibrate the P2D-ECT model. The Li-ion cell considered in this study is a second generation, 18,650 size commercial cell with a maximum capacity of 3.3 Ah. The cell consist of a graphite anode and Lithium Nickel Cobalt Aluminum Oxide (NCA) cathode. This NCA cathode has a layered structure. For the cell, the anode thickness is 92 106 m and the cathode thickness is 68 106 m. The electrodes are separated by a 16 106 m thick separator. The salt used in the cell is LiPF6 mixed in a nonaqueous mixture of ethylene carbonate (EC) and dimethyl carbonate (DMC). Maximum Li concentration in the salt is 1150 molm3. All the relevant parameters of the cell are given in Table 2. These parameters were obtained using experiments with high conﬁdence. The cell was tested at room temperature (T ¼ 298 K), high temperatures (T ¼ 318 K, T ¼ 333 K) and low temperatures (T ¼ 283 K, T ¼ 273 K, T ¼ 263 K) in a controlled instrumented laboratory setup. The cell was charged using a CC-CV protocol at different C-rates (0.1C, 0.3C, 0.5C, 0.7C and 1.0C) and discharged using a CC protocol at 0.5C. The CC charging phase was continued till 4.2 V and followed by a CV phase to ensure an equilibrated cell. The CC discharge was continued till 2.75 V. This charge-discharge protocol was followed at all the temperatures. For these tests the load current, cell voltage and electrode potentials were monitored at non-uniform frequency. The cell was also tested for US06 and UDDS protocols at the room temperature, however, electrode potentials were not monitored for these tests. In the test cases presented in this paper, the cell voltage measurements obtained for 0.5C discharge are assimilated for calibration. 4.4. Calibration at room temperature Total seven parameters are estimated for this test case. The estimated parameters are: the solid phase VOFs, solid phase Bruggeman's coefﬁcients, SEI resistance and the solid phase diffusion coefﬁcients for both the electrodes. The prior uncertainty in these parameters is speciﬁed using independent Gaussian distributions. The discharge C-rate is used as a control input. The measurement uncertainty is speciﬁed using se ¼ 0.005 V. This measurement uncertainty corresponds to the typical on-board voltage sensor accuracy. For this test case, 121 design points are generated using the Latin square sampling method. Total 50,000 samples are generated from the posterior distribution using the MCMC sampling after the burnout period of 20000 samples. The calibrated model is subsequently used for cell voltage prediction for charge, discharge and the drive-cycle protocols.

304

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

The calibrated parameters are summarized in the Table 5. Key observations from the estimated parameters are: 1) the cathode VOF is lower than the anode; and b) the estimated Bruggeman's coefﬁcients are higher than the popular choice of 1.5 used in the Liion battery simulations [10,3]. The lower VOF results in the faster ﬁlling up of the cathode during discharge (and corresponding faster emptying during charge) as compared to the anode. This indicates that the cathode electrode is limiting for this NCA/C cell. This observation is also corroborated by the experiments, as discussed in the following subsection. The framework estimates higher value for both the electrode Bruggeman's coefﬁcients. The higher Bruggeman's coefﬁcient indicates the higher electrode tortuosity. High tortuosity of the lithium-ion battery electrodes and the resultant deviation from the Bruggeman's relation is already established in the literature by several other authors (see, for eg. [43,49], and the references therein). These authors have used particle scale modeling to investigate tortuosity of the electrodes [43,49] and reported the Bruggeman's coefﬁcients signiﬁcantly higher than 1.5. The similar conclusion is obtained in this paper using the Bayesian framework at a signiﬁcantly lower computational cost. Though not reported in this paper for brevity, authors conducted a separate case study to estimate the electrolyte Bruggeman's coefﬁcient. The framework, however, estimated the value close to 1.5 for the electrolyte Bruggeman's coefﬁcient with high conﬁdence. Fig. 4 compares the cell voltage predicted using the base and calibrated P2D-ECT model with the experimental data. The comparison for different charge protocols is shown in the frame (a), while, the comparison for 0.5C discharge is shown in the frame (b). From the ﬁgure, it can be observed that the predictions obtained using the calibrated simulator matches closely with the experimental data. Inset of frame (b) compares predicted electrode potential with the experimental data. The electrode potential data was obtained using the three electrode experiments performed at the room temperature. Experimental data shows that the cathode is a limiting electrode for the present NCA/C cell. This behavior is appropriately captured by the calibrated simulator, whereas the base simulator wrongly predicts anode as the limiting electrode. In Fig. 5, the calibrated simulator prediction is compared with the drive cycle test data. The drive cycle tests are designed to simulate typical real world operating scenarios. In this paper, a combination of US06 and UDDS drive protocol is investigated. The US06 test protocol considered in this paper consists of an initial 0.9C CC charging followed by an urban drive cycle proﬁle. Note that the cell is not fully charged during the US06 drive cycle. The UDDS proﬁle is initialized from the ﬁnal SOC of the US06 drive cycle. Top panel of Fig. 5 shows the current proﬁle used for the drive cycle. The cell voltage comparison is shown in the middle panel. Comparison near the end of discharge is also shown in the insets. The bottom panel of the ﬁgure shows the absolute error between the predicted and the experimental cell voltage. The calibrated simulator successfully tracks the complete drive cycle proﬁle. The average cell voltage prediction accuracy is 99.98%. The maximum cell voltage prediction error during the US06 drive cycle protocol is 8.41% while the maximum prediction error during the UDDS protocol is 9.387%.

(a)

(b) Fig. 4. Calibration at room temperature. Frame (a) compares cell voltage for different charge proﬁles, and frame (b) compares cell voltage for 0.5C discharge proﬁle.

4.5. Calibration across temperatures Main aim of the calibration process is to obtain a unique set of parameters that work for all the operating conditions. Estimation of this unique set of parameters, however, is often very difﬁcult for a complex system like the lithium ion cell. Thus, it is essential to quantify the conﬁdence on the estimated parameters and the resultant calibrated simulator. In the proposed Bayesian framework, this conﬁdence is quantiﬁed using the posterior PDF of the variance s22 of the model structural uncertainty. In this paper, effect of the model uncertainty is investigated for calibration of the P2DECT model across a temperature range of 263 Ke333 K. Fig. 6(a) compares the prior and posterior PDFs of s22 for all the test cases considered in this subsection. In the ﬁrst test case, 0.5C discharge data at six different temperatures is used for the calibration. All the parameters are estimated using this discharge data. The resultant posterior PDF of s22 is shown in Fig. 6(a) (Post. T1). As

Table 5 Estimated parameters at room temperature. MAP estimate of parameters

Base MAP Estimate

Anode

Cathode

BRUG

BRUG

Diffus.

Diffus.

SEI

VOF

VOF

(An)

(Ca)

(An)

(Ca)

Res.

0.678 0.7535

0.7 0.6472

1.5 4.5588

1.5 3.2578

12

1.6e 3.66e13

13

2.6e 5.91e14

160e04 125e04

Current, A

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

can be observed from the ﬁgure, posterior PDF moves to the right of the prior PDF, giving higher probability to the larger values of s22 as compared to the prior. Thus, the Bayesian framework estimates the lower conﬁdence on the calibrated simulator after the data assimilation. The predicted cell voltage matches poorly with the experimental data for this test case. For brevity, this cell voltage comparison is not shown in this paper.

5 −10 Experiment

Vcell, V

4.5

Calibrated Simulation

3.5

Error

2.5 −4

10

−9

10

0

5000

Time, s

10000

15000

Fig. 5. Cell voltage comparison US06-UDDS drive cycle. Top panel of the ﬁgure shows applied load current, middle panel compares the predicted cell voltage with the measurements and the bottom panel shows the absolute error in the cell voltage prediction.

14

Prior Post. T1 Post. T2 Post. T3 Post. T4

12

P.D.F.

10

6

(25)

where the side reaction over-potential is given by

4

hs ¼ fs fe Rf jLi s :

2 0 0

0.2

0.4

2 2

0.6

0.8

1

σ

(a) 7

Anode Brug. Coef. Cathode Brug. Coef.

6 5 BRUG

4.5.1. Incorporation of new physics The reduced conﬁdence on the calibrated simulator necessitates addition of new temperature dependent physics in the P2D-ECT model. In this paper, two physical phenomena are investigated for temperature dependent behavior of the P2D-ECT model. First, the paper introduces a temperature dependence for the Bruggeman's coefﬁcient. To obtain this temperature dependence, anode and cathode Bruggeman's coefﬁcient are estimated for each temperature while maintaining other parameters at the values obtained from the room temperature calibration. Second, the lithium plating is incorporated in the existing P2DECT model. During low temperature charging, some of the lithium ions get deposited on the anode particles instead of intercalation, resulting in the formation of the lithium plate [50,51]. During discharge, the lithium plate dissolves ﬁrst followed by the deintercalation from the anode. During charging, the lithium plating is modeled as an anode side reaction with an appropriate Butler-Volmer reaction kinetics [51]. Direction of this side reaction is reversed during discharge to model the lithium plate dissolution. This side reaction is given by

aa F ac F exp exp ; jLi ¼ i h h 0;s s RT s RT s

8

4 3 2 1 260

305

270

280

290 300 310 Temperature, K

320

330

340

(b) Fig. 6. Frame (a) compares the prior and posterior probability distribution of s22 . The test cases described are: T1-calibration at all temperatures without new hypotheses; T2-calibration at T ¼ 333 K with temperature dependent Bruggeman's coefﬁcient; T3calibration at T ¼ 263 K considering temperature dependent Bruggeman's coefﬁcient while neglecting the Lithium plating; T4-calibration at T ¼ 263 K considering temperature dependent Bruggeman's coefﬁcient and Lithium plating. Frame (b) shows temperature dependence of the Bruggeman's coefﬁcient.

(26)

Note that the OCP of the lithium plate is zero, which is used in the deﬁnition of the hs. Due to the lower OCP, the dissolution is favored during discharge. The solid phase mass conservation is ensured during the dissolution by modifying the anode particle interface boundary conditions as

vcs ¼ 0: vr r¼Rs

(27)

Once the lithium plate is completely dissolved, the anode behavior is modeled using the main Butler-Volmer reaction given by the Eq. (7). A sigmoid function is used in the formulation to ensure a smooth transition between the dissolution and deintercalation reactions. The model incorporated with the new physics is calibrated using the cell voltage data across a range of temperatures. The conﬁdence on the calibrated simulator increases signiﬁcantly at high temperatures when the temperature dependence of Bruggeman's coefﬁcient is considered, as can be seen from Fig. 6(a) for 333 K (Post. T2). The conﬁdence at low temperature, however, is still low (Post. T3 Fig. 6(a) for 263 K). After incorporating the lithium plating, the conﬁdence on the calibrated simulator increases considerably at the low temperature as can be observed from Fig. 6(a) (Post. T4). The estimated temperature dependence of the Bruggeman's coefﬁcient is shown in the Fig. 6(b). The anode Bruggeman's coefﬁcient is higher than the cathode across the temperature range. The estimated anode Bruggeman's coefﬁcient shows the sigmoidal behavior with the values of BRUG ¼ 4.5 at high temperatures and BRUG ¼ 6.0 at low the temperatures with the transition near

306

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

T ¼ 283 K. The estimated cathode Bruggeman's coefﬁcient can be closely approximated using a quadratic function of the temperature. From Fig. 6(a), it can be concluded that the incorporation of the new physics improves credibility of the P2D-ECT model. Thus, the P2D-ECT model with the incorporated physics is expected to predict the cell voltage with high conﬁdence and accuracy. In Fig. 7, predicted cell voltage for 0.5C discharge protocol is compared against the experimental data across a temperature range of 333K263 K. The effect of the lithium plating can be observed near beginning of the discharge at 273 K and 263 K. The cell voltage prediction accuracy is improved across all the temperatures using the proposed framework. For all the test cases presented in the paper, the cell voltage prediction accuracy is greater than 98%. 4.6. Physical insights into the temperature dependent lithium ion cell behavior The calibrated P2D-ECT model predicts the cell voltage with high accuracy and credibility across a range of temperature. This allows inference of the internal unobservable cell behavior with high conﬁdence. In this paper, the internal cell behavior is investigated for 0.5C discharge across a temperature range of 333K263 K. The predicted internal cell behavior is summarized in Fig. 8. Fig. 8(a) shows the evolution of electrolyte phase concentration at the middle of the separator. At low temperature, the lithium ions initially accumulates in the separator region resulting in the steep concentration rise. This initial concentration rise is accentuated by the lithium plate dissolution at 273 K and 263 K. Subsequently, the electrolyte phase concentration drops rapidly at the low temperature. Inset of the ﬁgure shows distribution of the electrolyte phase concentration across the cell at 2500 s, which is near middle of the discharge. As can be observed from the inset, the concentration gradients are much steeper at lower temperatures. The steeper concentration gradient implies that the low temperature cell behavior is controlled by the electrolyte diffusion [52,53]. Fig. 8(b) shows the Butler-Volmer reaction rate at time 1000 s (top panel) and 3500 s (bottom panel). At higher temperatures, the reaction rate is approximately uniform across the electrode thickness. At the low temperature, however, the anode reaction rate shows quadratic behavior. At 1000 s, the reaction rate is maximum at the anode-separator interface for all the temperatures, while the reaction rate becomes negligible towards anode current collector at the low temperature (263 K). This implies majority of the deintercalation takes place for particles near the anode-separator

Experiment

Vcell, V

5 4

T=333K

3 0

5000

,V cell

V

5

5

4

T=318K

0

4

5000

0

T=298K

0

5000

5 T=273K

3 5000 Time, s

4 3

5 T=283K

3 0

Simulation (Calib)

3

5 4

Simulation (Base)

4

T=263K

3 5000 Time, s

0

5000 Time, s

Fig. 7. Comparison of cell voltage predictions with experimental observations for a temperature range of 333Ke263 K.

interface. Due to low electrolyte diffusion, the deintercalation is negligible at the low temperature near the anode current collector. With time, the location of maximum reaction rate moves towards anode current collector, as can be observed from the bottom panel of the ﬁgure. Fig. 8(c) shows the evolution of solid phase surface concentration at the anode (top panel) and the cathode (bottom panel) current collector end. In the ﬁgure, the solid phase concentration is normalized using the respective maximum electrode solid phase concentrations. As can be noted from the ﬁgure, the cathode particle reaches its maximum value for all the temperatures. Thus, cathode is the limiting electrode for the NCA/C cell investigated in this paper, which is also experimentally corroborated using electrode potentials at the room temperature (see inset of Fig. 4(b)). Note that at low temperature (273 K and 263 K), the deintercalation from the anode particles start only after the lithium plate is dissolved, as can be observed from the top panel of the ﬁgure. At 263 K, the cathode particle is full near end of discharge, however, the anode particle is not fully empty. This observation corroborate the ﬁndings in Fig. 8(a)e(b). As the reaction rate is negligible for this anode particle (Fig. 8(b)), the deintercalation is low implying the particle retain majority of its lithium. At 263 K, the cathode particle source majority of its lithium for intercalation from the electrolyte salt, resulting in the depleted electrolyte phase concentration observed in Fig. 8(a). Fig. 8(d) shows the anode solid phase potential for all the temperatures. Effect of the lithium plating can be observed at the low temperatures (273 K and 263 K), where the anode solid potential shows two distinct phases. The ﬁrst phase corresponds to the plate dissolution, while the second phase corresponds to the anode deintercalation. Since the lithium plate OCP is zero (see subsection 4.5.1), low solid phase potential sufﬁcient to overcome the electrolyte resistance is enough for plate dissolution. Note that the electrolyte resistance is high at 263 K due to low diffusivity. Thus the solid phase anode potential is comparatively higher at 263 K. 4.7. Relevance of the estimated Bruggeman's coefﬁcient Bruggeman's relation given by Eq. (2) with the coefﬁcient BRUG ¼ 1.5 is widely used in the lithium ion cell models to approximate tortuosity of the porous electrodes [10]. The Bruggeman's relation is derived using the differential effective medium theory [54] to obtain the effective transport properties for the porous media. The derivation inherently assumes spherical particles and non-interacting particulate phase. These assumptions, however, are not strictly satisﬁed by the lithium ion battery microstructure [55]. The inadequacies of the Bruggeman's relation for the lithium ion cell models are investigated by Ferguson and Bazant [56], where they explore strict upper and lower bounds for effective properties of the porous electrodes. As rigorously demonstrated in Refs. [56], the Bruggeman's relation is close to the upper bound for the effective properties which represents a highly conducting isotropic material interspersed by a disconnected particulate phase. The porous electrode microstructure, however, maintain a network of interconnected particles to form percolative pathways for conduction. Effect of these percolative pathways on the Bruggeman's relationship is investigated by Vadakkepatt et al. [57], where they introduce the term connectivity as an inverse of the tortuosity. For electrolyte phase conduction and diffusion, the effect of percolation is negligible [56], which corresponds to the high value of connectivity. Due to high connectivity, Bruggeman's relationship with BRUG ¼ 1.5 approximates the electrolyte phase ionic conductivity and diffusivity with high accuracy [57]. The Bayesian framework presented

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309 −5

1800

x 10

s

−1

10

−3

−2

1600 1400 T=333 T=318 T=298 T=283 T=273 T=263

1200 1000 800 0

2000

4000 Time, s

6000

Reaction Rate, mol m

Electrolyte Conc., mol m

307

8000

5

Time=1000 s

0 0

10

30

40

x 10 4 2 0 0

Time=3500 s

10

20 Particle Index

30

40

6000

8000

(b)

(a) 1

0.8 Anode

0.5 0 0

2000

4000 Time, s

6000

8000

1 0.5 0 0

Cathode

2000

4000 Time, s

6000

8000

(c)

Anode Solid Potential, V

Solid Phase Conc., mol m−3

20 Particle Index

−5

0.6 0.4 0.2 0 0

2000

4000 Time, s

(d)

Fig. 8. Internal cell behavior across a temperature range of 333Ke263 K. Frame (a) shows electrolyte phase concentration; frame (b) shows Butler-Volmer reaction rate; frame (c) shows the normalized solid phase surface concentration; frame (d) shows anode solid phase potential.

in this paper correctly estimates the electrolyte phase Bruggeman's coefﬁcient close to 1.5. Effect of percolation, however, is signiﬁcant for solid phase conduction, where a network of connected particles is essential for electron transport [58]. The non-conducting phase can percolate in the solid phase blocking the conducting pathways [56], resulting in the loss of connectivity. Appropriate modeling of the low connectivity requires value of the Bruggeman's coefﬁcient signiﬁcantly higher than 1.5 [57]. Several authors have used particle scale modeling to report high values of solid phase Bruggeman's coefﬁcient [43,49]. The Bayesian framework presented in this paper correctly estimates the high values for the solid phase Bruggeman's coefﬁcient. The connectivity of the solid phase decreases with the increase in the volume fraction [57]. Thus, the higher solid phase volume fraction implies higher Bruggeman's coefﬁcient. For the NCA/C cell investigated in this paper, the estimated anode volume fraction is higher than the cathode, corroborating the higher estimate of the anode Bruggeman's coefﬁcient by the proposed Bayesian framework. Although compensated using higher coefﬁcient, effect of percolation renders Bruggeman's relation inadequate for modeling of the solid phase effective conductivity [56]. To accurately capture the solid phase effective conductivity, accounting of the percolation effect is essential. The percolation is dependent on a critical volume fraction, known as the percolation threshold, below which the percolation effect is insigniﬁcant. Above the percolation threshold, the effective conductivity can be modeled using a simple power law [56], however the value of the exponent varies signiﬁcantly due to the tunneling effect [58]. The effective solid phase conductivity, thus, is critically dependent on the percolation threshold. This threshold is affected by the operating temperature [59], introducing temperature dependence on the percolation. This

temperature dependence of the percolation is captured by the temperature dependent Bruggeman's coefﬁcient introduced in this paper. 5. Concluding remarks A Bayesian framework is proposed in this paper for calibration of the P2D-ECT model. The framework exploit properties of the Matrix-variate Gaussian distribution and the Latin square sampling to obtain a computationally and numerically tractable formulation. Functionality of the framework is demonstrated for calibration of a NCA/C lithium ion cell model across a wide range of operating conditions and current protocols. The calibrated model predictions show excellent match with the experimental data, as demonstrated in this paper for cell voltage predictions. Output of the proposed framework is a calibrated model with quantiﬁed model uncertainty. In this paper, posterior variance of the model uncertainty is shown to quantify credibility of the calibrated simulator. Posterior PDF of this variance is exploited to test new physical phenomena of temperature dependent Bruggeman's coefﬁcients and lithium plating at low temperatures. When a combination of these two physics is considered, the model uncertainty is reduced for all the temperatures and operating conditions. The cell voltage prediction accuracy is greater than 98% for all the temperatures. The proposed framework is generic in nature and can be used for different lithium typologies and assimilation of other measurable quantities like cell temperature. The framework can also be used without any change if a higher ﬁdelity ECT model is made available. In the future, authors will demonstrate effectiveness of the proposed framework for different lithium typologies like Lithium Iron Phosphate (LFP), Lithium Nickel Cobalt Manganese

308

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309

Oxide (NCM) etc. In the future, authors will also explore different formulations for model uncertainty and its applicability for online prognostics of lithium ion battery. The authors will also investigate in detail the temperature dependence of the effective electronic conductivity of the porous electrode.

List of Symbols

Appendix A. Deﬁnition of the vectorization and Kronecker product operators For an arbitrary matrix

2

a11 A¼4 … am1

a12 … am2

… … …

3 a1n … 52ℛ mn ; amn

(A.1)

the vectorization operator is deﬁned as Bayesian D,Di H,hi p(,) T Vi x y

b

ε

li f q d(,) m x s2i Si

Framework data matrix regression functions probability density function system model row-covariance matrix of ith distribution control input observations regression coefﬁcients measurement uncertainty Li correlation length parameters of ith distribution hyper-parameters parameters discrepancy function mean of the distribution true system response variance of ith distribution column-covariance matrix of ith distribution

Electrochemical Thermal Model a speciﬁc surface area of the porous active material, m1 BRUG Bruggeman's coefﬁcient c concentration of lithium, mol m3 cs surface concentration of lithium, mol m3 D diffusion coefﬁcient of lithium, m2 s1 F Faraday's Constant, C I current density, A m2 L thickness of the electrode/separator/cell, m k reaction rate, m s1 j Butler-Volmer reaction, mol m2 s1 r radial coordinate, m R radius of the active particle within the electrode, m Rg Ideal gas constant, J (mol K)1 tþ transference number of the lithium ions t time, s T temperature, K U open circuit voltage, V V cell voltage, V VOF Volume fraction x cartesian coordinate ε porosity k electrical conductivity of the electrolyte phase, S m1 s electrical conductivity of the solid phase, S m1 f potential, V in,ip interfacial variables, subscript max maximum value of the variable, subscript s side reaction variable, subscript

Appendix B. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.jpowsour.2016.04.106.

vecðAÞ ¼ ½ a11 a21 … am1 … a1n a2n … amn T 2ℛ mn1 : (A.2) For arbitrary matrices A2ℛ mn and B2ℛ pq

2

a11 A¼4 … am1

a12 … am2

… … …

3 2 b11 a1n … 5; B ¼ 4 … bp1 amn

b12 … bp2

… … …

3 b1q … 5; bpq (A.3)

the Kronecker product A5B is deﬁned as

2

a11 B A5B ¼ 4 … am1 B

a12 B … am2 B

… … …

3 a1n B … 52ℛ mpnq : amn B

(A.4)

References [1] B. Kennedy, D. Patterson, S. Camilleri, J. Power Sources 90 (2000) 156e162. [2] H. Xiaosong, M. Nikolce, L. Johannesson, B. Egardt, Appl. Energy 111 (2013) 1001e1009. [3] V.S. Kumar, J. Power Sources 222 (2013) 426e441. [4] T. Trucano, L. Swiler, T. Igusa, W. Oberkampf, M. Pilch, Rel. Eng. Syst. Saf. 91 (2006) 1331e1357. [5] S. Blattnig, L.L. Green, J. Luckring, J. Morrison, R. Tripathi, Towards a credibility assessment of models and simulations, in: 2008 American Institute of Aeronautics and Astronautics, USA, 2008. [6] W. Oberkampf, T. Trucano, C. Hirsch, Appl. Mech. Rev. 57 (2004) 345e384. [7] W. Oberkampf, M. Barone, J. Comp. Phys. 217 (2006) 5e36. [8] C. Fellner, J. Newman, J. Power Sources 85 (2000) 229e236. [9] S.M. Rezvanizaniani, Z. Liu, Y. Chen, J. Lee, J. Power Sources 256 (2014) 110e124. [10] G.L. Plett, Battery Management Systems, 2014. [11] A. Seaman, T.S. Dao, J. McPhee, J. Power Sources 256 (2014) 410e423. [12] V. Johnson, J. Power Sources 110 (2002) 321e329. [13] R. Kroeze, P. Krein, in: IEEE Power Electronics Specialist Conference, 2008. [14] M. Chen, G. Rincon-Mora, IEEE Trans. Energy Conv. 21 (2006) 504e511. [15] X. Hu, S. Li, H. Peng, F. Sun, J. Power Sources 239 (2013) 449e457. [16] X. Hu, S. Li, H. Peng, J. Power Sources 198 (2012) 359e367. [17] M. Doyle, T.F. Fuller, J. Newman, J. Electrochem. Soc. 140 (1993) 1526e1533. [18] T. Fuller, M. Doyle, J. Newman, J. Electrochem. Soc. 141 (1994) 1e10. [19] L. Rao, J. Newman, J. Electrochem. Soc. 144 (1997) 2697. [20] M. Guo, G. Sikha, White, J. Electrochem. Soc. 158 (2011) A122eA132. [21] K. Smith, C. Rahn, C. Wang, Energy Conv. Manag. 48 (2007) 2565e2578. [22] P. Northrop, V. Ramadesigan, S. De, V. Subramaniam, J. Electrochem. Soc. 158 (2011) A1461eA1477. [23] J. Reimers, J. Electrochem. Soc. 160 (2013) A811. [24] Y. Hu, S. Yurkovich, Y. Guezennec, B. Yurkovich, Control Engg. Pract. 17 (2009) 1190e1201. [25] X. Tang, X. Mao, J. Lin, B. Koch, Li-ion battery parameter estimation for state of charge, in: 2011 American Control Conference, San Antonio, Texas, USA, 2011. [26] H. Eichi, M. Chow, Adaptive parameter identiﬁcation and state-of-charge estimation of lithium-ion batteries, in: 38th Annual Conference of the IEEE Industrial Electronics Society, Montreal, Canada, Oct. 25-28, 2012, p. 2012. [27] R. Masoudi, T. Uchida, J. McPhee, J. Power Sources 291 (2015) 215e224. [28] J. Forman, S. Moura, J. Stein, H. Fathy, Genetic parameter identiﬁcation of the Doyle-Fuller-Newman model from experimental cycling of a LiFePO4 battery, in: 2011 American Control Conference, San Antonio, Texas, USA, 2011. [29] J. Reimers, M. Shoesmith, Y. Lin, L. Valoen, J. Electrochem. Soc. 160 (2013) A1870eA1884. [30] V. Ramadesigan, V. Boovaragavan, M. Arabandi, K. Chen, H. Tsukamoto, B. R.D, V. Subramanian, ECS Trans. 19 (2009) 11e19. [31] V. Ramadesigan, K. Chen, N. Burns, V. Boovaragavan, B. R.D, V. Subramanian,

P. Tagade et al. / Journal of Power Sources 320 (2016) 296e309 J. Electrochem. Soc. 158 (2011) A1048eA1054. [32] S. Santhanagopalan, Q. Zhang, K. Kumaresan, R. White, J. Electrochem. Soc. 155 (2008) A345eA353. [33] T. Huria, G. Ludovici, G. Lutzemberger, J. Power Sources 249 (2014) 92e102. [34] B. Saha, K. Goebel, S. Poll, J. Christophersen, IEEE Trans. Inst. Meas. 58 (2009) 291e296. [35] M. Kennedy, A. OHagan, J. R. Stat. Soc. B 63 (2001) 435e464. [36] D. Higdon, M. Kennedy, J.C. Cavendish, J.A. Cafeo, R.D. Ryne, SIAM J. Sci. Comp 26 (2) (2004) 448e466. [37] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, J. Tu, Technometrics 49(2). [38] S. Conti, A. OHagan, J. Stat. Plan. Infer. 140 (2010) 640e651. [39] O. Stegle, C. Lippert, J.M. Mooij, N.D. Lawrence, K.M. Borgwardt, Efﬁcient inference in matrix-variate gaussian models with iid observation noise, in: Advances in Neural Information Processing Systems, 2011. [40] P.M. Tagade, B.M. Jeong, H.L. Choi, Build. Enviorn. 70 (2013) 232e244. [41] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, J. Chem. Phy. 21 (1953) 1087e1092. [42] S. Chib, E. Greenberg, Americ. Statist. 49 (1995) 327e335. [43] I. Thorat, D. David, E. Stephenson, N. Nathan, A. Zacharias, K. Zaghib, J. Harb, D. Wheeler, J. Power Sources 188 (2009) 592e600. [44] P. Gambhire, N. Ganesan, S. Basu, K. Hariharan, S. Kolake, T. Song, D. Oh, T. Yeo, S. Doo, J. Power Sources 290 (2015) 87e101. [45] A. Gupta, T. Varga, J. Multivar. Anal. 41 (1992) 80e88. [46] D.C. Montgomery, Design and Analysis of Experiments, Wiley, New York,

309

1984. [47] M.D. McKay, R.J. Beckman, W.J. Conover, Technometrics 21 (2) (1979) 239e245. [48] J. D. Cook, Inverse gamma distribution, online: http://www.johndcook.com/ inverse.gamma.pdf, Tech. Rep. [49] A. Gupta, J. Seo, X. Zhang, W. Du, A. Sastry, W. Shyya, J. Electrochem. Soc. 158 (2011) A487eA497. [50] P. Arora, M. Doyle, R. White, J. Electrochem. Soc. 146 (1999) 3543e3553. [51] R.D. Perkins, A. Randall, X. Zhang, G. Plett, J. Power Sources 209 (2012) 318e325. [52] S. Basu, R.S. Patil, S. Ramachandran, K.S. Hariharan, S.M. Kolake, T. Song, D. Oh, T. Yeo, S. Doo, J. Power Sources 283 (2015) 132e150. [53] P. Tagade, K.S. Hariharan, P. Gambhire, S.M. Kolake, T. Song, D. Oh, T. Yeo, S. Doo, J. Power Sources 306 (2016) 274e288. [54] D. Stroud, Superlattices Microstruct. 23 (1998) 567e573. [55] D.W. Chung, M. Ebner, D.R. Ely, V. Wood, R.E. Garca, Modeling and simul, Mater. Sci. Eng. 21 (2013) 074009. [56] T.R. Ferguson, M.Z. Bazant, J. Electrochem. Soc. 159 (2012) A1967eA1985. [57] A. Vadakkepatt, B. Trembacki, S. Mathur, J. Murthy, J. Electrochem. Soc. 163 (2016) A119eA130. [58] A. Awarke, S. Lauer, S. Pischinger, M. Wittler, J. Power Sources 196 (2011) 405e411. [59] R. Parthasarathy, X.-M. Lin, K. Elteto, T.F. Rosenbaum, H. Jaeger, Phys. Rev. Lett. 92 (2004) 076801.