Pr<R

Sbbxu^S

Contents lists available at ScienceDirect

Journal of Power Sources

journal homepage: www.elsevier.com/locate/jpowsour

On-board monitoring of 2-D spatially-resolved temperatures in cylindrical lithium-ion batteries: Part I. Low-order thermal modelling

CrossMark

Robert R. Richardson, Shi Zhao, David A. Howey*

Department of Engineering Science, University of Oxford, Oxford, UK

HIGHLIGHTS

> Derivation and validation of low-order 2-D thermal model for cylindrical cells.

> Efficient numerical implementation based on Chebyshev spectral-Galerkin method.

Includes anisotropic heat conduction

and inhomogeneous convection

boundary conditions.

Applicable to various battery cooling

configurations, such as side or end

cooling.

Suitable for use in a state-estimator with surface temperature or EIS measurements.

GRAPHICAL ABSTRACT

ARTICLE INFO

Article history: Received 22 February 2016 Received in revised form 6 June 2016 Accepted 24 June 2016

Keywords: Lithium-ion battery Thermal model Spectral methods Low-order modelling Temperature estimation

ABSTRACT

Estimating the temperature distribution within Li-ion batteries during operation is critical for safety and control purposes. Although existing control-oriented thermal models - such as thermal equivalent circuits (TEC) - are computationally efficient, they only predict average temperatures, and are unable to predict the spatially resolved temperature distribution throughout the cell. We present a low-order 2D thermal model of a cylindrical battery based on a Chebyshev spectral-Galerkin (SG) method, capable of predicting the full temperature distribution with a similar efficiency to a TEC. The model accounts for transient heat generation, anisotropic heat conduction, and non-homogeneous convection boundary conditions. The accuracy of the model is validated through comparison with finite element simulations, which show that the 2-D temperature field (r, z) of a large format (64 mm diameter) cell can be accurately modelled with as few as 4 states. Furthermore, the performance of the model for a range of Biot numbers is investigated via frequency analysis. For larger cells or highly transient thermal dynamics, the model order can be increased for improved accuracy. The incorporation of this model in a state estimation scheme with experimental validation against thermocouple measurements is presented in the companion contribution (Part II).

© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

DOI of original article: http://dx.doi.org/10.1016/jjpowsour.2016.06.104.

* Corresponding author. E-mail addresses: robert.richardson@eng.ox.ac.uk (R.R. Richardson), shi.zhao@ eng.ox.ac.uk (S. Zhao), david.howey@eng.ox.ac.uk (D.A. Howey).

Lithium-ion batteries generate heat due to electrochemical processes, which results in internal temperature gradients during operation. In a typical usage scenario, such as a standard vehicle

http://dx.doi.org/10.1016/jjpowsour.2016.06.103

0378-7753/© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

drive cycle, cells may experience temperature differences between surface and core of 20 °C or more [1]; and during a rapid overheating event this discrepancy can be as large as 40—50 °C [2]. High battery temperatures could trigger thermal runaway resulting in fires, venting and electrolyte leakage. While such incidents are rare [3], consequences include costly recalls and potential endanger-ment of human life. Consequently, transient thermal modelling of batteries during operation is an essential requirement for battery management systems to ensure safe and optimal performance.

In this study, we present and validate a low-order thermal model of a cylindrical battery cell, capable of capturing 2-D thermal dynamics. The model is based on the spectral-Galerkin method, achieving high accuracy with minimal computational requirements, making it suitable for online applications. The remainder of this paper is organised as follows. In Section 2, the alternative numerical methods for low-order thermal modelling are discussed. In Section 3, a gentle introduction to Galerkin-spectral methods is provided by means of a toy problem - a 1-D heat equation in Cartesian coordinates. In Section 4, the full 2-D thermal model is presented; and in Section 5, the results of the model are validated through comparison with high fidelity Finite Element (FE) simulations. Matlab code to simulate the presented model is available online.1

2. Low-order thermal modelling

Lumped parameter thermal equivalent circuit (TEC) models are perhaps the most popular approach for efficient thermal modelling. These methods have been used extensively for low-order and control-oriented modelling of battery cells and packs [1,4—10]. Their main advantage is that they are simple to implement. However, they are unable to predict the full temperature field throughout the domain of interest; their outputs merely consist of nodal values representing average temperatures. This is a particular limitation when comparison with temperature measurements at discrete locations is necessary for state or parameter estimation. Moreover, since the parameters of the model have no physical meaning, they require parameterization using experimental data for any given set of parameters and operating conditions. On the other hand, TEC models with parameters that can be directly calculated from physical properties have been used extensively in thermal modelling of electric machines and other applications [11—13]. However, these are known to have poor performance for large Biot numbers, whilst increasing the number of elements to improve accuracy comes at the expense of increased computational complexity [14,15].

Physics based models instead solve the underlying diffusion partial differential equation (PDE) governing the heat transfer. They are generally applicable to a broad range of problems, and can predict the full temperature field throughout the domain of interest. Several studies have presented 2-D or 3-D thermal simulations of battery cells [16—19]. However, the solution is typically obtained using computationally intensive numerical methods such as Finite Difference Methods (FDM) or Finite Element Methods (FEM), and so their potential for application in control systems is limited. Analytical solutions have also been developed [20,21], however these are inappropriate for on-line applications since time domain solutions rely on computationally intensive integral transforms.

To reduce the computational burden of physics based models, low order approaches have been proposed using techniques such as balanced truncation [22,23] and polynomial approximation (PA) [24—27]. However, these methods are only suited to 1-D problems

involving infinite or semi-infinite domains, or symmetric boundary conditions [28]. This may be acceptable for cases in which thermal gradients arise predominantly in one direction, such as in air cooling of small form factor (e.g. 18650 or 26650) cylindrical cells. However, large form factor cells have a greater propensity for thermal gradients in multiple directions. Recent studies have used 45 Ah cylindrical cells with a diameter of 64 mm and height of 198 mm [29]. Such dimensions give rise to much larger Biot numbers and hence more significant thermal gradients. Moreover, certain cooling configurations, such as end cooling via cooling plates, are more likely to give rise to axial variations. These systems also introduce additional complexities, such as the potential for each surface of the cell to be in contact with a different cooling fluid with a different free stream temperature and/or convection coefficient. Fig. 1a shows a case with forced convection via a liquid coolant at the base of the cell and natural convection to the air at the other surfaces. Fig. 1b shows another scenario in which unequal cooling of multiple cells can result in 2-D thermal dynamics as heat is transferred axially through the tabs from one cell to a neighbouring cell. Consequently, there is a clear need for low-order thermal models capable of capturing 2-D thermal dynamics.

Spectral methods are an alternative numerical method for solving PDEs, in which the spatial discretization is carried out using global rather than local approximating functions [30,31]. In general, FDM and FEM methods are suitable for complex geometries, whereas spectral methods provide greater computational efficiency at the expense of model flexibility and the assumption of a smooth solution. Spectral methods are a type of weighted residual method - a group of approximation techniques in which the solution errors are minimized in a certain way - and they are classified according to the minimization technique employed. The most common techniques are spectral-collocation and spectral-Galerkin. In a collocation method, the solution is obtained by interpolating an approximating function at a set of domain nodes, whereas in a Galerkin method, the solution is obtained by forcing the residual of an integral multiplied by a test function to zero.

Spectral methods have been used in previous studies for low-order battery modelling [32—35], and recently they have even been applied to 2-dimensional problems [36]. However, to the authors' knowledge, no study has applied a Galerkin method to 2-D battery thermal problems. This is perhaps due to the complexity of accounting for non-homogeneous boundary conditions, in particular convection boundary conditions, using a Galerkin method. One of the advantages of the Galerkin method is that the basis functions implicitly satisfy the boundary conditions, and so it possible to have very low order models with satisfactory accuracy.

In this paper, we present a 2-D model of a cylindrical cell based on a Chebyshev spectral-Galerkin (SG) method. The model accounts for transient heat generation, anisotropic heat conduction,

h forced

1 www.github.com/robert-richardson/Spectral-Thermal-Model-2D.

Fig. 1. Battery cooling configurations resulting in 2-D thermal dynamics: (a) Plate cooling, (b) unequal cooling with intercell heat transfer.

and non-homogeneous convection boundary conditions, with different convection coefficients and external temperatures at each surface. This generality makes it suitable for simulating various battery cooling configurations, such as those discussed previously. The treatment of the non-homogeneous boundary conditions is achieved by an efficient boundary lifting algorithm, adapted from a recently developed boundary lifting algorithm for elliptic problems [37]. Since the underlying basis functions implicitly satisfy the boundary conditions, there is no need for additional equations to impose the boundary conditions. As a result accurate models can be generated using as a few as two basis functions in each direction, and thus 22 = 4 states. The accuracy can be improved arbitrarily by increasing the number of basis functions in either dimension.

3. Toy problem

Firstly, we introduce SG methods by means of a simple example - a 1-D heat equation in Cartesian coordinates on the domain be [-1,1], with non-homogeneous boundary conditions.

3.1. Governing equation

functions, however, for this example, Chebyshev functions were used since they can be conveniently defined such that they adhere to any of a broad class of boundary conditions [30]. Suitable basis functions are found as follows. We denote by Cn the nth degree Chebyshev polynomial of the first kind and then let {fn(b)}N=0 be a basis function such that

fn = + ZnCn+1 + hnCi+2 (V), (8)

where zn and hn are defined according to the formula in Lemma 4.3 of [30] such that they adhere to the boundary conditions,

Zn = {4(n + 1)(a+b_ + a-b+)}/DETn, (9)

hn = { _ 2a_a+ + (n2 + (n + 1)2) (a+b_ _ a_b+)

+ 2b_b+n2(n + 1)2}/DETn, (10)

The governing equation is given by:

VT , d2T

pcPVt - k " q = 0

where t is time, be [_1,1] is the position coordinate and k is the thermal conductivity. The convection boundary conditions are given by

a.T + b.-^ T = er at x = 1 + +9x r

a-T + b-T = e< at x = —1

a+ = h/k, b+ = 1, a— = —h/k, b— = 1,

er = a+Ttx, el = a—Ttx,

where h is the convection coefficient and Tœ is the external temperature. For simplicity the convection coefficient and external temperature at each boundary are assumed to be the same.

3.2. Finite sum approximation

The starting point of the SG method is to approximate the solution T of eq. (1) by a finite sum

T = J2 anfni xl + TM,

n=0 ^ '

DETn = 2a+a- + ((n + 1)2 + (n + 2)2) (a-b+ - a+b-)

- 2b-b+(n + 1 )2(n + 2)2. (11)

From this point on, fn are considered known functions.

3.4. Residual equation

Substituting eq. (7) for T into eq. (1) leads to the residual:

VT . V T R := PCp^z — k-^ — qs0. vt dx2

The principle of the SG method is to force an integral of the residual to zero by requiring

pCpVf — k ^br — q I ndx = 0,

where n is a test Junction.

Substituting eq. (7) for T and fn for n into eq. (13), we have

J f0fndx, —, J fNfndx

d dta0

1 1 / f...../ f

— q J fndx = 0. —1

where an are unknown solution coefficients, and fn are the trial (or basis) functions which must satisfy the Robin boundary conditions of eqs. (2) and (3).

3.3. Basis functions

3.5. State space equation

With n = 0,1,...,N, the above N + 1 equations representing the spatial discretization can be written in compact form as Ex = Ax + Bu, (15)

In principle any Jacobi polynomial may be used for the basis

x = [a0, ..., aN]', u = [q, TM]',

E(i,j) = PCp J fifjdx,

, 1 92fj _ A(i,j) = k / -jdx,

B(i, 1)= / fidx,

(20) (21)

B(i, 2) = 0.

Letting the outputs of the system be the temperature at the left and right boundaries, T(b = -1) and T(b = 1), we have

y = Cx + Du,

y = [Tx=-1, Tx=^T,

x = -1 , -, $N\ x = -1

x = 1 , -, fN X = 1

0 1 0 1

As we shall see, the actual 2-D problem of interest presented in the following section is more complex than the above problem, for the following reasons: (i) it involves cylindrical rather than Cartesian coordinates, (ii) the physical domain is not necessarily in the range [-1,1] and must be scaled accordingly, and (iii) it is 2-D, which renders the homogenization step of the non-homogeneous boundary conditions non-trivial. However, the solution processes in both cases are equivalent.

4. Thermal model

4.1. Overview

We now describe the full 2-D SG model. A step-by-step procedure is given as follows. In Section 4.2, the 2-D thermal problem is defined. In Section 4.3, the model is scaled from the physical coordinates to a coordinate system suitable for implementation of the SG algorithm. In Section 4.4, the scaled model is decomposed into a homogeneous problem and a (time-invariant) boundary lifting function. Section 4.5 presents the solution to the homogeneous problem using the SG method and Section 4.6 presents the calculation of the boundary lifting function. Finally, Section 4.7 presents the overall solution to the original problem which combines the solutions obtained from Sections 4.5 and 4.6. The overall procedure

ultimately generates a state-space model; the resulting model can then be solved efficiently online.

4.2. Model definition

The model consists of the transient energy conservation equation in cylindrical coordinates. Heat generation is assumed to be uniform in space but time-dependant, but as we show in part II of this paper, the error introduced by this assumption is modest. The multi-layer structure of the battery is treated as a homogeneous solid with anisotropic thermal conductivity in the radial and axial directions. The temperature variation in the azimuthal (f) direction is neglected. Convective heat transfer is assumed to occur at the outside surfaces, and the properties of the heat transfer fluid (i.e. the heat transfer coefficient and the fluid free-stream temperature) may be different for each surface (see Fig. 2).

The resulting model is governed by the following 2-D boundary value problem [28]:

dT d2T kr dT d2T _

pCp~dt - krd* - T a? - kz= q

where t is time and r and z are the position coordinates in the radial and axial directions respectively. The functions T(r,z,t) and q(t) are the temperature distribution and volumetric heat generation rate, respectively. The parameters p and cp are the density and specific heat capacity respectively, and kr and kz are the anisotropic thermal conductivities in the r and z directions. The boundary conditions are given by:

dT hr (_ _ .

ar=-k;(T - Tt^at r=rout

dT hl , ,

vr=kr(T - ^ at r = rin

-h- (T - Tr,„) at z = H

dT h dz

f(T - at z = 0

(28c) (28d)

Çoo, h,

Tr,„ hr

Fig. 2. Schematic of cylindrical cell geometry for the thermal model, showing different convection coefficients and/or external temperatures at each surface.

where (rs>œ; s = t,b,r and /} are the free-stream temperatures of the heat transfer fluid at the top, bottom, right and left surfaces,2 and {hs; s = t>b>r and /} are the corresponding convection coefficients.

4.3. Change of sca/e

In order to exploit the orthogonality properties of the polynomial basis functions it is necessary to scale from the original (physical) domain, (r,z) where r 2 [rin> rout], z 2 [0, H], to the spectral domain, (r,z) where r2[-1, 1], ze[-l, 1]. We define the spectral coordinates by

C-T + d-—T--dZ

where,

ei at Z = -1,

a+ = hr/kr, b+ = a,

a— = —h,/kr, b— = a,

c+ = ht/kz, d+ = ß,

c— = —hb/kz, d— = ß,

et = a+Tt,x, eb = a—T

er = c+Tr,œ, ei -- = c—Ti

b,œ ;

4.4. Homogenization of the boundary conditions

, r - rin rout — rir

- 1 s.t.

Z (r = rin) = — 1 Z (r = rout ) = 1

~ 2z _ Z = - — 1s.t.

Z(z = 0) = —1 z(z = H) = 1 •

Rearranging (29a) and (29b), the original coordinates are given

1 + Z + arin

Z +1 ß ■

(30a) (30b)

where a and ß are scaling factors defined by a = 2/(rout - rin) and ß = 2/H.

Substituting (30a) and (30b) into the original eq. (27), we obtain the governing equation in the spectral domain:

dT 2, d2T

pcn--a kr—=7

^ P dt r dZ2

dt dZ2 1 + arin + r d r

with the boundary conditions:

a+T + b+T = et at Z= 1

+ +d r

dT —

a-T + b-—T--d r

eb at Z = —1

c+T + d+— T = er at Z= 1

+ +d z r

(32a) (32b) (32c)

Next, the problem with non-homogeneous boundary conditions is transformed into a problem with homogeneous boundary conditions using a boundary lifting algorithm, as follows: We set

T = T + Te,

where Te(r, z) is an arbitrary function satisfying the original boundary conditions (see later), and ~(r, z, t) is an auxiliary function satisfying the modified problem

- 2| d2~

pcP--a kr—--

r P dt r dZ2 1

— ß2k^ v z2

v r2 1 + arin + r v r subject to the homogeneous boundary conditions

a±T + b±T d T =

± ± d r

c±T + d±T—T =

± ± dz

@ 2, d Te

q — I — a2kr—2e

1 + rin + r d r

^e — ß2kzd2Te

(36a) (36b)

The boundary lifting is similar to that in Ref. [37], except that here we apply it to cylindrical coordinates and neglect the corner component of the lifting (since we are assuming constant external temperatures).

4.5. Chebyshev-Galerkin approximation

Now we would like to solve the modified problem (35) using the Galerkin method. The first step is to multiply the equation by a test function, n, and integrate over the entire domain,

1 + arin + z

dT d2T a2kr

pcP--a kr—--„ „

H p dt r dz2 1 + arin + r d r

dT — ß2k2^ v z2

1 + arin + z

2 Note that although we have prescribed constant (w.r.t the spatial variable) external temperatures at each side, the model is also capable of handling spatially dependant temperatures (e.g. Tr„(z) or Tt>TC(r)). However, since the external temperatures would rarely be known with such a high fidelity in a real application, we chose constant values for the sake of simplicity.

where (f, n) denotes the integral of f weighted by n, throughout the domain r, z,

(f, n)= J J f (V, b) , b) drdz.

Note the inclusion of the r terms as defined in (30a) on each side of eq. (38) to account for cylindrical coordinates.

The second step is to approximate the solution ~(b, b, t) with a finite number of functions

T = J2 J2xkjfk[b, a±, b± ) fb (b, c±, d± k=0 j=0

where xkj are unknown solution coefficients, and f? and f? are the basis functions which must satisfy the homogeneous boundary conditions (36a) and (36b) respectively. For simplicity, the same number of basis functions, N, are chosen for the radial and axial directions but it should be noted that different values could be chosen if one dimension required greater resolution than the other.

Suitable basis functions can be found as follows. First we find basis functions in the radial and axial directions separately. Each of these is obtained in a similar manner to that of the toy problem (eqs. (8)—(11)): we denote by Ck the kth degree Chebyshev poly? ^ N

nomial of the first kind and then let ff ?(?, a±, b±)}^0 be a basis in the radial direction such that

fb = + CkCk+ï(^j + ^+2 (V),

where Zk(a±, b±) and r[(a±, b±) are defined as before such that they adhere to the boundary conditionsNin the radial direction.

Similarly, we let {fb(b,c±,d±)}. 0 be a basis in the axial direction such that

fj = Cj(b) + jCj+1 (b) + h/Cj+2 (V),

where Z] (c±, d±) and rf (c±, d±) adhere to the boundary conditions in the axial direction.

As mentioned previously, the choice of basis function is not limited to Chebyshev polynomials: in general, various Jacobi polynomials may work equally well. In this case, we verified that almost identical results were obtained using Legendre polynomials.

follows. For the right side, we have

c+Te + ^VTe~er at z = +1 Vz

where ~ denotes weakly satisfying conditions. Substituting (43) in (44) gives

J2(dk(c+z + d+) + c+z2 + 2d+b) ) fb

+ £ bd)11 + b2 d k0

c+fb + d+—ir J dz

V B b f

Substituting for b = 1 and noting that the term c+f + d+f

is equal to zero (since the basis functions by definition satisfy the homogeneous boundary conditions), eq. (45) simplifies to

N , , ~

J2{d{(c++ d+)+dk(c + + 2d+)) fk ~er.

The weak boundary condition can now be replaced by the appropriate integral equation,

(41) '(d[(c+ + d+) + dk(c+ + 2d+)Vf, f[\ = 1er, ff\ ,

k=0 \ / b \ / b

where (f, g)^ denotes the integral, 1

(f, g)b =

1 + r + arir a

f b gr db.

Through a similar process for the left, top and bottom sides, we find that

Xt{dk(-c- + d-) + dk(c- - 2d-)) f ei,

4.6. Boundary lifting function

We now wish to determine the boundary lifting function, Te(b,b). Recall that Te(b,b) must satisfy the original (non-homogeneous) boundary conditions (32a), such that when subtracted from the actual temperature, T (eq. (34)), the resulting auxiliary function, satisfies the homogeneous boundary conditions. Hence, the aim of this section is to derive a function which is constant with respect to time and satisfies the non-homogeneous boundary conditions. Since this cannot be achieved exactly, the boundary conditions must instead be satisfied in a weak sense (i.e. the solution converges as more basis functions are included). To begin, we adopt a similar approach to [37] by assuming a form of the solution as follows:

j2Uz+dkz2) fk k0

£ j + d)Vb2) f

' (dj11 (a+ + b+) + djV(a+ + 2b+)) f fb ^ b = { et , fb

XT(dj11(_a_ + b_) + djV(a_ _ 2b_)) f ff^ =(eb, ff^ ,

where (f g)b denotes the integral, 1

(f ,g)b = J f (b)g(b) dz. 1

The conditions at the vertical and horizontal sides are defined as

Equations (47) and (49) and equations (50) and (51) form two linear systems which can be each solved for the corresponding ds vectors:

I _ k4er ~ k2el fprl \-1srl

M4 - k2k3(

Ml - k3er

M4 - k2k3

J4et - J2 eb f

J1J4 - J2J3 V

Prl srl

dn^ = j4et ~ J2ebiptb\ _1stb

(53a) (53b) (53c)

EX = Ax + Bu

x =(x00; x10, ■■■xN0, x01, x11, ■■•, xN1,x0N,xNN ) , (58)

and u = [q(t),1]T. The system matrices are defined as follows. First, let us denote the column vector

r fZ, f1 f0, fN f0, f0 ff, f1 ff , fNff, f0 fN, fN fN

d^ = J1eb - J3et fptb\-1stb

J1J4 - JJ V ' '

where fdg = (dg, df, ..., d|)J ; s = I, II, III and IV} are vectors of unknown expansion coefficient, fsg = (sg, sf, ..., sN )T ; s = rl, and tb} are vectors of known source terms, and

k1 = c+ + d+, j1 = a+ + b+,

k2 = c+ + 2d+, j2 = a+ + 2b+,

k3 = d-- c-, J3 = b-- a_,

k4 = c-- 2d-; J4 = a-- 2b-.

The P matrix for the right and left edges is defined as {Pri = pfk; i, k = 0,1, ..., Ng, where

prk =/ f fkdb,

and that for top and bottom sides as {Ptb = ptb; i,j = 0,1,...,Ng, where

ptb = I ff f dz.

Lastly, the source terms are defined as

fss = sf; i = 0,1,...,N; s = rl and tbg, where

-rl = r I1 + z + arin Affdr,

stb _ 5i

ff dz.

Thus, an explicit expression for the unknown expansion coefficients, ds is obtained, and so the boundary lifting function, Te, can be considered known at this stage.

4.7. Solution algorithm

Finally, we present the solution algorithm to the modified problem. Equation (35) can be expressed in state space form, as follows:

and take n = ji, where ji is the ith element of the vector J; then:

-, /1 + r + arin 1 1 A

E(i;J) = Pcp I-a-~jj, j I,

where E(i, j) denotes the element in the ith row and jth column of the matrix E,

A(i, j):

1 + r + arin

-, V2jJ 2 92 jJ o^kr—Xj + ß2kf^j

+ ^r-vrj, ji

for i, J = 0,1,...,N, and, B(i, 1)

l + T^, ji A,

B(i, 2)

1 + r + arin

o2kr + ß2kf ^

1 9Te ,

+ ak^~9rr, ji

for i = 0,1,...,N.

The complete solution to the original non-homogeneous problem is then given by:

T r, z, t = T r, z, t + Te r, z

We choose as outputs the temperatures at the bottom-centre, left-centre, top-centre and right-centre locations, i.e.

Fig. 3. Schematic of cell showing the model outputs displayed in Figs. 4 and 5.

T1 = T(? = -1,? = 0), T2 = T(? = 0 ,? = -1), T3 = T(? = 1,? = 0) and T4 = T(? = 0 , ? = 1) (see Fig. 3). Thus, the outputs are given by:

y = Cx + Te

where y = (r^,T3,r4)T, Te = (Te,1,Te,2Je,3Je,4)T' and,

C(1 , j) , C(2 , j), C(3 J)

C(4 , j),

= jjiî =—1,z = 0j, = jj (r = 0,x = —1), = jj (r = 1,x = 0

= jj(r = 0,x = 1

for j = 0,1, ...,N.

Note that the mean temperature, T, may also included as an output (T is used in Part II of this paper for computing the overall cell electrochemical impedance). Hence, an additional row is appended to the C matrix,

H Zout

C(5 -j)=1 r2 2 r2 / i rjj(r,Z)drdZ, n 'out 'in i J

for j = 0,1, ...,N, which in the scaled coordinates becomes 1 1

C(5 , j) = Tß

1 + r + arit

1 —1

-jj^r, z^j d?dz.

An additional element must also be included in the boundary lifting function,

H Z out

Te,5 = 1 2 2 2 i i rTe(r,z)drdz,

H rlut — rl 0 /

0 ' in

which in the scaled domain becomes 1 1

Te,5 = aß

1 + r + arit

Te xr xz dxrdxz.

With these equations, the mean temperature is computed as the fifth output.

The frequency domain response of the above linear system, H(s), is calculated by

H(s) = C(sI — A)—1B

where s = ju is the Laplace variable and I is the identity matrix.

The above algorithm was implemented both numerically (using clenshaw-curtis quadrature) and analytically using the Matlab Symbolic Maths Toolbox. The numerical implementation allows the state matrices to be generated more efficiently than the symbolic approach, although the resulting state space model is identical (and therefore equally efficient) in each case.

5. Results and discussion

To validate the SG model, the results were compared with high fidelity FEM simulations, implemented using the Matlab Partial Differential Equation Toolbox. To ensure the accuracy of the FEM solution, a fine mesh consisting of 3760 elements was used. The

Table 1

Thermophysical properties for model validation.

Parameter Symbol Value

Inner radius rin 4 mm

Outer radius rout 32 mm

Height H 198 mm

Density r 2,118 kg m—3

Specific heat capacity CP 765 J kg—1 K—1

Radial thermal conductivity kr 0.66 W mK—1

Axial thermal conductivity kz 66 W mK—1

time step for both the SG and FEM models was set at 1 s.

The thermo-physical parameters chosen for the model validation are shown in Table 1. The dimensions were chosen to match those of the large format lithium-ion cell employed in Ref. [29], and the remaining thermal parameters were chosen based on typical properties of lithium iron phosphate cells ([22,25,38]).

5.1. Time domain

Time domain simulations were carried out using two different cooling scenarios, as shown in Table 2. Case 1 represents forced convection air cooling, with equal temperatures and convection coefficients at each of the external sides. Case 2 represents forced convection liquid cooling at the left end of the cell (for instance, via a cooling plate [39]), and mild forced convection to the ambient air at the remaining faces. Thus, the temperature at the left face is set to 3 °C and a large convection coefficient typical of forced cooling via water or glycol is applied, whereas the remaining faces are exposed to a small convection coefficient at 18 °C. This case was chosen to highlight the ability of the model to account for different external temperatures and/or convection coefficients at each side. Note that the convection coefficient at the bottom side (the inner radius of the jelly roll) was set to zero in both cases since negligible cooling occurs here in a typical thermal management system; however a non-zero value could easily be applied if it were required.

To ensure highly transient conditions with large internal temperature gradients, pulsed power load profiles with large heat generation rates but relatively short durations were applied. The load profiles and resulting temperature distributions for the two cases are compared against the corresponding FEM solutions in Figs. 4 and 5. The locations of the model outputs are shown in Fig. 3.

The results show that the SG method with Ns = 4 states (i.e. 2 basis functions in each of the radial and axial directions) is capable of accurately capturing the temperatures in both cases. Slightly greater accuracy is achieved using Ns = 9 states (3 basis functions in each direction). In Case 1, the temperatures, T1 and T3 (i.e. the first and third model outputs from eq. (65)) are plotted (see Fig. 4a) since the temperature gradient is primarily in the radial direction in this case. In Case 2, the temperatures T2 and T4 were plotted since a large gradient occurs in the axial direction in this case. Plots of T(r)z=H/2 are shown in part (b) of these figures. In each case, the first two sub-plots show the distribution at instances of highly transient

Table 2

Convection coefficients for the two cooling scenarios.

Case 1 Case 2

h (W m—2) TM (°C) h (W m—2) TM (°C)

Left 100 18 400 3

Right 100 18 30 18

Top 100 18 30 18

Bottom 0 18 0 18

Fig. 4. Comparison of SG vs. FEM results for Case 1: forced air convection. (a) Evolution of temperatures, T1 (core) and T3 (outer surface), with pulsed load profile shown in subplot; (b) temperature distribution along the centre-line (z = H/2) of the cell at denoted times; (c) Contour plot of the temperature distribution at 1060 s.

dynamics (i.e. just after a load is applied or removed), which is more difficult to capture accurately for a low-order thermal model, although the results are still in good agreement in both cases. The third sub-plot in each case shows the solution at a later time - when the problem is dominated by slow thermal dynamics - and the SG and FEM solutions are in nearly perfect agreement at these times. The full 2-D contour plots are shown in part (c) of each figure.

Again, the solution is plotted at an instant with transient dynamics to demonstrate the ability of the low-order model to accurately simulate the temperature field under these conditions. The SG model is in good agreement with the FEM solution, although the advantage of increasing the number of model states is apparent in Fig. 5(c), as the result using Ns = 9 states is in better agreement with that of the FEM solution than the result with Ns = 4 states.

r (m) r (m) r (m)

Fig. 5. Comparison of SG vs. FEM results for Case 2: end-plate liquid cooling. (a) Evolution of temperatures, T2 and T4, with pulsed load profile shown in subplot; (b) temperature distribution along the centre-line (z = H/2) of the cell at denoted times; (c) Contour plot of the temperature distribution at 240 s.

Lastly, we note that the model with 4 states is of a similar order to an equivalent circuit thermal model and so could be applied with similar computational efficiency. Moreover, a similar order model implemented using a spectral-collocation method would not give results as accurate as those presented in this section, since in the SG method the boundary conditions are implicitly satisfied by the basis functions, whereas additional equations are required to enforce the boundary conditions in the collocation case.

5.2. Frequency domain

In this section we compare the frequency response of the low-order SG models against a baseline solution obtained by using an SG model with a large number of states (Ns = 225). Specifically, we examine the impact of changes in heat generation on the core temperature, 71, given by the transfer function H(s) = T1(s)/q(s). This is calculated using eq. (65), using only the rows and columns of the B and C matrices corresponding to the heat generation input and T1 output.

The parameters chosen are the same as those in Case 1 of the previous section. However, a variety of different Biot numbers (Bi = h(rout - rin)/fcb) are obtained by varying the magnitude of the convection coefficients on each side. Convection coefficients of h = {10, 50,100, 500} Wm-2 are chosen, resulting in Biot numbers of Bi = {0.42, 2.10, 4.20 and 21.02}, respectively. The size of the model is varied by increasing the number of basis functions in both the radial and axial directions. Models with 1, 2, 3 and 5 basis functions in each direction are chosen, resulting in Ns = 1, 4, 9 and 25 states respectively.

Fig. 6(a)—(d) show the magnitude of the frequency response in

the range f = 1 x 10-4 to 1 x 10° Hz for each of the four Biot numbers, along with the error of the low order models relative to the high fidelity solution. These plots show that as the model order (i.e. the number of states) is increased, the magnitude of the error is reduced. Moreover, for all cases, there is a critical frequency at which the error becomes non-negligible, and this critical frequency increases as the model order is increased. We also note that the error increases as (i) the Biot number increases, and (ii) the perturbation frequency increases. These trends are as expected. Thus, if a particular application involves larger Biot numbers or higher frequencies (for instance, due to larger cells, more aggressive drive cycles or higher performance cooling), the model order could be increased accordingly to achieve a required accuracy.

6. Conclusions

Computationally efficient thermal models are necessary for control-oriented applications. The model presented in this paper is of a similar order to a thermal equivalent circuit (TEC) model - and so could be applied with similar computational efficiency - but has much greater spatial resolution. It could therefore provide greater accuracy than a TEC if used as part of a state-estimation scheme. Moreover, although we have presented here a model for a cylindrical cell, it can easily be modified to apply to 2-D simulation of prismatic cells. However, it is difficult to apply the SG method to more involved problems, due to the complexity of it's implementation. Hence, for configurations involving several cells or nonuniform cooling, TECs may remain favourable due to their relative simplicity.

In the companion contribution [reference to Part II], the model is

Fig. 6. Frequency response, T (s)/q(s), of the SG method with different model orders for a range of Biot number conditions as indicated. A large order (Ns = 225) SG method is used as the baseline solution.

incorporated into a state estimation scheme and the predicted temperatures at four locations (one internal and three on the cell surface) are validated experimentally against thermocouple measurements.

Acknowledgements

This work was funded by a NUI Travelling Scholarship, a UK EPSRC Doctoral Training Award, the Foley-Bejar scholarship from Balliol College, University Of Oxford, and the RCUK Energy Pro-grammes's STABLE-NET project (ref. EP/L014343/1).

Nomenclature

A, B, C, D, E system matrices

ai polynomial coefficients

cP specific heat capacity [J kg~f K_1 ]

Ck kth deg. Chebyshev poly. of llst kind

h convection coefficient [W m~2]

H cell height

H(s) system transfer function

I identity matrix

j imaginary number

kr radial thermal conductivity [W m~f K~f

kz axial thermal conductivity [W m~f K~f ]

q volumetric heat generation [W]

Ns number of model states

r radial coordinate [m]

b scaled radial coordinate [m]

rin inner radius [m]

rout outer radius [m]

s Laplace variable

T temperature [K]

T auxiliary temperature function [K]

Te boundary-lifting function [K]

t time [s]

x system state

y system output

z axial coordinate [m]

z scaled axial coordinate [m]

Abbreviations

BVP Boundary Value Problem

FDM Finite Difference Method

FEM Finite Element Method

PA Polynomial approximation

PDE Partial Differential Equation

SG Spectral-Galerkin

TEC Thermal equivalent circuit

n test function

p density [kg m~3]

f azimuthal coordinate

f basis function

J vector of basis functions

u frequency

Subscripts

b bottom edge

l left edge

r right edge

t top edge

tc ambient condition

Appendix A. Supplementary data

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

References

[1] C. Forgez, D. Vinh Do, G. Friedrich, M. Morcrette, C. Delacourt, Thermal modeling of a cylindrical LiFePO4/graphite lithium-ion battery, J. Power Sources 195 (9) (2010) 2961-2968, http://dx.doi.org/10.1016/jjpows-our.2009.10.105. http://linkinghub.elsevier.com/retrieve/pii/ S037877530901982X.

[2] N.S. Spinner, C.T. Love, S.L. Rose-Pehrsson, S.G. Tuttle, Expanding the operational limits of the single-point impedance diagnostic for internal temperature monitoring of lithium-ion batteries, Electrochim. Acta (2015), http:// dx.doi.org/10.1016/j.electacta.2015.06.003. http://linkinghub.elsevier.com/ retrieve/pii/S0013468615013602.

[3] Q. Wang, P. Ping, X. Zhao, G. Chu, J. Sun, C. Chen, Thermal runaway caused fire and explosion of lithium ion battery, J. Power Sources 208 (2012) 210-224, http://dx.doi.org/10.1016/jjpowsour.2012.02.038. http://linkinghub.elsevier. com/retrieve/pii/S0378775312003989.

[4] X. Li, F. He, L. Ma, Thermal management of cylindrical batteries investigated using wind tunnel testing and computational fluid dynamics simulation, J. Power Sources 238 (2013) 395-402, http://dx.doi.org/10.1016/jjpows-our.2013.04.073. http://linkinghub.elsevier.com/retrieve/pii/ S037877531300671X.

[5] X. Lin, H.E. Perez, J.B. Siegel, A.G. Stefanopoulou, Y. Ding, M.P. Castanier, Parameterization and observability analysis of scalable battery clusters for onboard thermal management, in: Les Rencontres Scientifiques dIFP Energies nouvelles, 2011 no. December.

[6] X. Lin, H.E. Perez, S. Mohan, J.B. Siegel, A.G. Stefanopoulou, Y. Ding, M.P. Castanier, A lumped-parameter electro-thermal model for cylindrical batteries, J. Power Sources 257 (2014) 1-11, http://dx.doi.org/10.1016/ j.jpowsour.2014.01.097. http://linkinghub.elsevier.com/retrieve/pii/ S0378775314001244.

[7] R. Mahamud, Advanced Battery Thermal Management for Electrical-drive Vehicles Using Reciprocating Cooling Flow and Spatial-resolution, Lumped-capacitance ThermalModel (Ph.D. thesis), University of Nevada, 2011.

[8] N. Damay, C. Forgez, M.-P. Bichat, G. Friedrich, A. Ospina, Thermal modeling and experimental validation of a large prismatic Li-ion battery, in: IECON 2013-39th Annual Conference of the IEEE Industrial Electronics Society, 2013, pp. 4694-4699, http://dx.doi.org/10.1109/IEC0N.2013.6699893. http:// ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6699893.

[9] H. Dai, L. Zhu, J. Zhu, X. Wei, Z. Sun, Adaptive Kalman filtering based internal temperature estimation with an equivalent electrical network thermal model for hard-cased batteries, J. Power Sources 293 (2015) 351-365, http:// dx.doi.org/10.1016/j.jpowsour.2015.05.087. http://linkinghub.elsevier.com/ retrieve/pii/S0378775315009891.

[10] F. He, L. Ma, Thermal management of batteries employing active temperature control and reciprocating cooling flow, Int. J. Heat Mass Transf. 83 (2015) 164-172, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.11.079. http:// linkinghub.elsevier.com/retrieve/pii/S0017931014010783.

[11] P. Mellor, D. Roberts, D. Turner, Lumped parameter thermal model for electrical machines of TEFC design, in: IEE Proceedings B (Electric Power Applications), vol. 138, IET, 1991, pp. 205-218.

[12] R. Wrobel, P.H. Mellor, A general cuboidal element for three-dimensional thermal modelling, IEEE Trans. Magn. 46 (8) (2010) 3197-3200, http:// dx.doi.org/10.1109/TMAG.2010.2043928. http://ieeexplore.ieee.org/lpdocs/ epic03/wrapper.htm?arnumber=5512838.

[13] N. Simpson, R. Wrobel, P.H. Mellor, A general arc-segment element for three-dimensional thermal modeling, Magn., IEEE Trans. 50 (2) (2014) 265-268.

[14] F. Qi, A. Stippich, M. Guettler, M. Neubert, R.W. De Doncker, Methodical considerations for setting up space-resolved lumped-parameter thermal models for electrical machines, in: Electrical Machines and Systems (ICEMS), 2014 17th International Conference on, IEEE, 2014, pp. 651-657.

[15] N. Simpson, R. Wrobel, P.H. Mellor, An accurate mesh-based equivalent circuit approach to thermal modeling, Magn. IEEE Trans. 50 (2) (2014) 269-272.

[16] G.H. Kim, A. Pesaran, R. Spotnitz, A three-dimensional thermal abuse model for lithium-ion cells, J. Power Sources 170 (2) (2007) 476-489, http:// dx.doi.org/10.1016/j.jpowsour.2007.04.018.

[17] A. a. Pesaran, Battery thermal models for hybrid vehicle simulations, J. Power Sources 110 (2) (2002) 377-382, http://dx.doi.org/10.1016/S0378-7753(02) 00200-8. http://linkinghub.elsevier.com/retrieve/pii/S0378775302002008.

[18] W.B. Gu, C.Y. Wang, Thermal-electrochemical modeling of battery systems, J. Electrochem. Soc. 147 (8) (2000) 2910, http://dx.doi.org/10.1149/1.1393625.

[19] V. Srinivasan, C.Y. Wang, Analysis of electrochemical and thermal behavior of Li-ion cells, J. Electrochem. Soc. 150 (1) (2003) A98, http://dx.doi.org/10.1149/ 1.1526512.

[20] K. Shah, S. Drake, D. Wetz, J. Ostanek, S. Miller, J. Heinzel, a. Jain, an experimentally validated transient thermal model for cylindrical Li-ion cells, J. Power Sources 271 (2014) 262-268, http://dx.doi.org/10.1016/jjpows-our.2014.07.118. http://linkinghub.elsevier.com/retrieve/pii/ S0378775314011720.

[21] K. Shah, A. Jain, Modeling of steady-state and transient thermal performance of a Li-ion cell with an axial fluidic channel for cooling, Int. J. Energy Res. 39 (4) (2015) 573—584.

[22] M. Muratori, N. Ma, M. Canova, Y. Guezennec, A model order reduction method for the temperature estimation in a cylindrical Li-ion battery cell, in: ASME 2010 Dynamic Systems and Control Conference, vol. 1, 2010, pp. 633—640, http:||dx.doi.org|10.1115|DSCC2010-4200.

[23] M. Muratori, M. Canova, Y. Guezennec, G. Rizzoni, E. Politecnico, V. Lambruschini, A reduced-order model for the thermal dynamics of Li-ion battery cells, in: 6th IFAC Symposium Advances in Automotive Control, Munich, Germany, 2010, pp. 10—15.

[24] Y. Kim, J.B. Siegel, A.G. Stefanopoulou, A computationally efficient thermal model of cylindrical battery cells for the estimation of radially distributed temperatures, in: American Control Conference (ACC), 2013, Washington, DC, 2013, pp. 698—703. http:||ieeexplore.ieee.org|xpls|abs_all.jsp? arnumber=6579917.

[25] Y. Kim, S. Mohan, S. Member, J.B. Siegel, A.G. Stefanopoulou, Y. Ding, The estimation of temperature distribution in cylindrical battery cells under unknown cooling conditions, IEEE Trans. Control Syst. Technol. (2014) 1—10.

[26] R.R. Richardson, P.T. Ireland, D.A. Howey, Battery internal temperature estimation by combined impedance and surface temperature measurement, J. Power Sources 265 (2014) 254—261, http:||dx.doi.org|10.1016|j.jpows-our.2014.04.129. http:||linkinghub.elsevier.com|retrieve|pii| S0378775314006302.

[27] R.R. Richardson, D.A. Howey, Sensorless battery internal temperature estimation using a Kalman filter with impedance measurement, IEEE Trans. Sustain. Energy 6 (4) (2015). http:||ieeexplore.ieee.org|xpls|abs_all.jsp? arnumber=7097077&tag=1.

[28] D.W. Hahn, M.N. Ozisik, Heat Conduction, John Wiley & Sons, 2012.

[29] V. Roscher, M. Schneider, P. Durdaut, N. Sassano, S. Pereguda, E. Mense, K.-R. Riemschneider, Synchronisation using wireless trigger-broadcast for impedance spectroscopy of battery cells, in: Sensors Applications Symposium (SAS), 2015 IEEE, IEEE, 2015, pp. 1—6.

[30] J. Shen, T. Tang, L.-L. Wang, Spectral Methods: Algorithms, Analysis and

Applications, Springer, 2011.

[31] L.N. Trefethen, Spectral Methods in MATLAB, vol. 10, Siam, 2000.

[32] V.R. Subramanian, V.D. Diwakar, D. Tapriyal, Efficient macro-micro scale coupled modeling of batteries, J. Electrochem. Soc. 152 (10) (2005) A2002, http://dx.doi.org/10.1149/1.2032427. http://jes.ecsdl.org/cgi/doi/10.1149/!. 2032427.

[33] P.W.C. Northrop, V. Ramadesigan, S. De, V.R. Subramanian, Coordinate transformation, orthogonal collocation, model reformulation and simulation of electrochemical-thermal behavior of lithium-ion battery stacks, J. Electrochem. Soc. 158 (12) (2011) A1461, http://dx.doi.org/10.1149/ 2.058112jes. http://jes.ecsdl.org/cgi/doi/10.1149/2.058112jes.

[34] L. Cai, R.E. White, Lithium ion cell modeling using orthogonal collocation on finite elements, J. Power Sources 217 (2012) 248-255, http://dx.doi.org/ 10.1016/j.jpowsour.2012.06.043. http://linkinghub.elsevier.com/retrieve/pii/ S0378775312010439.

[35] A. Bizeray, S. Zhao, S. Duncan, D. Howey, Lithium-ion battery thermal-electrochemical model-based state estimation using orthogonal collocation and a modified extended kalman filter, J. Power Sources 296 (2015) 400-412.

[36] P.W.C. Northrop, M. Pathak, D. Rife, S. De, S. Santhanagopalan, V.R. Subramanian, Efficient simulation and model reformulation of two-dimensional electrochemical thermal behavior of lithium-ion batteries, J. Electrochem. Soc. 162 (6) (2015) A940-A951, http://dx.doi.org/10.1149/ 2.0341506jes. http://jes.ecsdl.org/cgi/doi/10.1149/2.0341506jes.

[37] E.H. Doha, A.H. Bhrawy, An efficient direct solver for multidimensional elliptic Robin boundary value problems using a Legendre spectral-Galerkin method, Comput. Math. Appl. 64 (4) (2012) 558-571, http://dx.doi.org/10.1016/ j.camwa.2011.12.050. http://dx.doi.org/10.1016/j.camwa.2011.12.050.

[38] M. Fleckenstein, S. Fischer, O. Bohlen, B. Bäker, Thermal impedance spectroscopy-a method for the thermal characterization of high power battery cells, J. Power Sources 223 (2013) 259-267.

[39] A. Jarrett, I.Y. Kim, Design optimization of electric vehicle battery cooling plates for thermal performance, J. Power Sources 196 (23) (2011) 10359-10368, http://dx.doi.org/10.1016/jopowsour.2011.06.090. http:// linkinghub.elsevier.com/retrieve/pii/S0378775311013279.