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

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).


Introduction
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 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 1 E-mail: {robert.richardson, shi.zhao, david.howey} @eng.ox.ac.uk.
temperatures could trigger thermal runaway resulting in fires, venting and electrolyte leakage. While such incidents are rare [3], consequences include costly recalls and potential endangerment 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 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,5,6,7,8,9,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,12,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,17,18,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,25,26,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 methoda 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.

T , h forced
Spectral methods have been used in previous studies for low-order battery modelling [32,33,34,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 cylindri-cal cell based on a Chebyshev spectral-Galerkin (SG) method. The model accounts for transient heat generation, anisotropic heat conduction, 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 treatement 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 2 2 ≡ 4 states. The accuracy can be improved arbitrarily by increasing the number of basis functions in either dimension.

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

Governing equation
The governing equation is given by: where t is time,x ∈ [−1, 1] is the position coordinate and k is the thermal conductivity. The convection boundary conditions are given by where 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.

Finite sum approximation
The starting point of the SG method is to approximate the solution T of eq. (1) by a finite sum where a n are unknown solution coefficients, and φ n are the trial (or basis) functions which must satisfy the Robin boundary conditions of eq. (3).

Basis functions
In principle any Jacobi polynomial may be used for the basis 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 C n the nth degree Chebyshev polynomial of the first kind and then let {φ n (x)} N n=0 be a basis function such that where ζ n and η n are defined according to the formula in Lemma 4.3 of [30] such that they adhere to the boundary conditions, and where From this point on, φ n are considered known functions.

Residual equation
Substituting eq. (7) for T into eq. (1) leads to the residual : The principle of the SG method is to force an integral of the residual to zero by requiring where ν is a test function. Substituting eq. (7) for T and φ n for ν into eq. (13), we have State space equation With n = 0, 1, . . . , N , the above N + 1 equations representing the spatial discretisation can be written in compact form as where and Letting the outputs of the system be the temperature at the left and right boundaries, where 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 homogenisation step of the non-homogeneous boundary conditions non-trivial. However, the solution processes in both cases are equivalent.

Overview
We now describe the full 2-D SG model. A step-bystep 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.

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 (ϕ) 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 Figure 2). The resulting model is governed by the following 2-D boundary value problem [28]: 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 ρ and c p are the density and specific heat capacity respectively, and k r and k z are the anisotropic thermal conductivities in the r and z directions. The boundary conditions are given by: where {T σ,∞ ; σ = t, b, r and l} are the free-stream temperatures of the heat transfer fluid at the top, bottom, right and left surfaces 3 , and {h σ ; σ = t, b, r and l} are the corresponding convection coefficients.

Change of scale
In order to exploit the properties of the polynomial basis functions it is necessary to scale from the original (physi- We define the spectral coordinates bŷ Rearranging (29a) and (29b), the original coordinates are given by where α and β are scaling factors defined by α = 2/(r out − r in ) and β = 2/H. Substituting (30a) and (30b) into the original eq. (27), we obtain the governing equation in the spectral domain: with the boundary conditions: where,

Homogenization of the boundary conditions
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 where T e (r,ẑ) is an arbitrary function satisfying the original boundary conditions (see later), andT (r,ẑ, t) is an auxiliary function satisfying the modified problem subject to the homogeneous boundary conditions where The boundary lifting is similar to that in [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).

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, ν, and integrate over the entire domain, where (f, ν) denotes the integral of f weighted by ν, throughout the domainr,ẑ, 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 solutionT (r,ẑ, t) with a finite number of functions where x kj are unknown solution coefficients, and φr k and φẑ j 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 C k the kth degree Chebyshev polynomial of the first kind and then let {φr k (r, a ± , b ± )} N k=0 be a basis in the radial direction such that where ζr k (a ± , b ± ) and ηr k (a ± , b ± ) are defined as before such that they adhere to the boundary conditions in the radial direction. Similarly, we let {φẑ j (ẑ, c ± , d ± )} N j=0 be a basis in the axial direction such that where ζẑ j (c ± , d ± ) and ηẑ j (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.

Boundary lifting function
We now wish to determine the boundary lifting function, T e (r,ẑ). Recall that T e (r,ẑ) must satisfy the original (nonhomogeneous) boundary conditions (32a), such that when subtracted from the actual temperature, T , (eq. 34) the resulting auxiliary function,T , 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: The conditions at the vertical and horizontal sides are defined as follows. For the right side, we have where ≈ denotes weakly satisfying conditions. Substituting (43) in (44) gives Substituting forẑ = 1 and noting that the term is equal to zero (since the basis functions by definition satisfy the homogeneous boundary conditions), eq. (45) simplifies to The weak boundary condition can now be replaced by the appropriate integral equation, Through a similar process for the left, top and bottom sides, we find that (51) where f, g ẑ denotes the integral, Equations (47) and (49) and equations (50) and (51) form two linear systems which can be each solved for the corresponding d σ vectors: where {d σ = (d σ 0 , d σ 1 , ..., d σ 2 ) T ; σ = I, II, III and IV} are vectors of unknown expansion coefficient, {s σ = (s σ 0 , s σ 1 , ..., s σ N ) T ; σ = rl, and tb} are vectors of known source terms, and The P matrix for the right and left edges is defined as Thus, an explicit expression for the unknown expansion coefficients, d σ is obtained, and so the boundary lifting function, T e , can be considered known at this stage.
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, for j = 0, 1, ..., N , which in the scaled coordinates becomes An additional element must also be included in the boundary lifting function, which in the scaled domain becomes 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 where s = jω 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.

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 3,760 elements was used. The 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 [29], and the remaining thermal parameters were chosen based on typical properties of lithium iron phosphate cells ( [22,38,25]).

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.  100  18  400  3  Right  100  18  30  18  Top  100  18  30  18  Bottom 0  18  0  18 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 Figures 4 and 5. The locations of the model outputs are shown in Figure 3.
The results show that the SG method with N s = 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 N s = 9 states (3 basis functions in each direction). In Case 1, the temperatures, T 1 and T 3 (i.e. the first and third model outputs from eq. (65)) are plotted (see Fig-r9 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 Figure 5(c), as the result using N s = 9 states is in better agreement with that of the FEM solution than the result with N s = 4 states. 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.

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 (N s = 225). Specifically, we examine the impact of changes in heat generation on T 1 , given by the transfer function H(s) = T 1 (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 show the magnitude of the frequency response in the range f = 1 × 10 −4 to 1 × 10 0 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.

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.
Morevover, 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 non-uniform cooling, TECs may remain favourable due to their relative simplicity.
In the companion contribution [40], the model is 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.