A Computationally Efficient Coupled Electrochemical-Thermal Model for Large Format Cylindrical Lithium Ion Batteries

Wepresentaone-dimensional,radial,coupleddegradation-electrochemical-thermal(DET)modelofalargeformatcylindricallithiumioncell.Themodelconsistsofreducedorderequationsthatdescribetheelectrochemicalphenomena,includingthatassociated withdegradation,coupledwithanapproximatemodelofthermalbehavior.Thereducedorderelectrochemicalmodel,whichisapproximatedfromthepseudo-two-dimensional(P2D)electrochemicalmodelusingaPadéapproximationmethod,computesthe variationofelectrochemicalvariablesandheatgenerationterms.Simultaneously,acoupledthermalmodelcomputesthetemperaturedistributionintheradialdirectionofthecell.TheresultsfromDETmodelcomparefavorablytothoseobtainedfromsolvingthe 1Dradialcoupleddegradation-electrochemical-thermalpartialdifferentialequationsinCOMSOLMultiphysics,howevertheDETmodelreturnstheseresultsinsigniﬁcantlyreducedcomputationaltimes.Importantly,themodelcapabilityinprovidinginsightful informationofcelldegradationandtemperatureinacomputationallyefﬁcientmannerpavesthewayforthehealth-conscious,real-timeoptimalcontroloflargeformatcylindricalcells.

Lithium ion batteries are now commonly used for the storage of energy in large, renewable energy generation systems, such as wind and photovoltaic systems, to eliminate their inherent intermittency.This is due to the fact that the battery storage systems are capable of providing a rapid response to counteract the fluctuations and filter out the variabilities associated with renewable generation and therefore enhance stabilize grid performance and maximize system security benefits. 1arge format cylindrical lithium ion cells with high energy have been developed and implemented in those energy storage systems.For example, "CH75" cylindrical cells developed by Hitachi have capacity of 75Ah/cell. 2These cells can compose a large energy battery pack using relatively small number of cells, which can reduce the total number of components and therefore, increase the reliability of the pack. 2 In order to operate batteries optimally, with the intention of maintaining safety and extending battery life, a battery management system (BMS) is essential. 3Two critical functions of a BMS are the thermal management and degradation monitoring.These functions become more essential for large format cylindrical batteries in which nonuniformities in temperature and degradation occur during operation. 4urthermore, such non-uniformities exacerbate further degradation and temperature gradients within the battery.Information on heat generation is fundamentally important for managing thermal issues such as thermal runaway, electrical cell imbalance within the battery pack and poor performance at low temperatures. 5However, heat generation inside a cell is a complex process that requires the knowledge of the physical characteristics of the cell during its operation.Total heat generation is due to irreversible and reversible processes, Joule heating in the solid and electrolyte and heating from the electrode/current collector contact resistance. 6It is accompanied by changes in the electrochemical properties of the cell, such as entropy, solid and electrolyte concentrations due to chemical and electrochemical reactions and changes in the solid and electrolyte potentials.Battery degradation is mainly caused by the formation and growth of the solid electrolyte interphase layer (SEI) layer, which scavenges active lithium ions and electrolyte materials and increases battery resistance, leading to capacity and power fade, respectively. 7SEI growth is also coupled with the thermal behavior of the cell as higher temperatures increase the SEI growth rate which in turn causes higher SEI resistance and ohmic heat generation, which elevates further the temperature. 8Current BMSs rely significantly on equivalent circuit models due to their simplicity and low computational requirement. 9However, equivalent circuit models have limited insight on the electrochemical characteristics of a battery and therefore, given the intimate coupling between the two, they are not able to model battery thermal behavior precisely. 6lternatively, the pseudo-two-dimensional electrochemical model (P2D), first developed by Doyle, Fuller, and Newman, 10,11 does provide insight into battery electrochemical behavior coupled with thermal properties. 6The P2D model, however, consists of five partial differential equations (PDEs) and one algebraic equation, which cannot be practically implemented in embedded BMS applications, without simplification or order reduction, due to their high computational requirement. 9A simplified electrochemical model that has low computational overheads, whilst maintaining precision in a specific range of operation would therefore be ideal to facilitate accurate, real-time resolution and control of battery operation.][19] However, the SPM model assumes that the electrolyte concentration is constant and that the current in the electrolyte does not vary spatially, which results in poor voltage prediction capability. 12Cai and White 20 developed a reduced order model using proper orthogonal decomposition to compute the electrochemical variables at discretized locations along the model domains.Subramanian et al. 13 developed approximations of the microscale diffusion of lithium ion in the solid phase, which have been well-adapted in macroscale models in the literature.In another paper, Subramanian et al., 21 used finite difference approximations and polynomial representations to reduce a system of 12 coupled PDEs into a system of Differential Algebraic Equations (DAEs), which are amenable to be used in battery control.Smith et al. 22 derived analytic transfer functions, in the Laplace domain, for a number of electrochemical variables in a linearized P2D model.Lee et al. 9 proposed an improvement of this approach by developing a complete set of transcendental transfer functions for the linearized P2D model.These authors then used the discrete-time realization algorithm (DRA) 23 to convert these transfer functions into an optimal discrete-time statespace model which can be used in real-time applications.However, extended computational time is required to rerun the DRA to produce a new state-space model when the electrochemical parameters of the P2D change. 9Alternatively, Padé approximations provide a way of greatly simplifying complicated transcendental transfer functions in order to produce simpler rational polynomial transfer functions, which can be used in embedded applications. 24Padé approximants were first used by Fathy et al. 25 to simplify a single transcendental transfer function of Jacobsen and West's model, 26 namely the lithium ion concentration in the solid phase.Our previous work extends this outcome by presenting a reduced order model for a lithium ion battery in which Padé approximants are used to simplify all of the transcendental transfer functions associated with the variables of the linearized P2D model developed by Lee et al. 9 Our results showed that the reduced model delivers high accuracy for low computational overhead when compared to the P2D model.Importantly, the reduced model requires relatively simple rational functions to be computed, which can evaluate internal electrochemical variables at any spatial location within the cross section of the cell.Therefore, this model is ideal for coupling with thermal and degradation models of the system.
Coupled thermal and electrochemical models have been presented previously in the literature on battery thermal management. 27,28apron et al. 29 investigated the thermal behavior of large lithium ion phosphate battery using coupled thermal-electrochemical model and implemented it in COMSOL Multiphysics. 30Pals and Newman 31 presented a one-dimensional, coupled thermal-electrochemical model that calculates the heat generation rate and temperature of a single cell and then in a second paper 32 they use their model to determine the temperature distribution within a cell stack.Wang et al. 33 modelled the thermal performance of a battery pack with various cell configurations in the pack and with different cooling strategies.However, due to their complexity and computational cost these aforementioned models are implemented using powerful commercial software, embodying integrated finite element and computational fluid dynamical packages [34][35][36] such as COMSOL Multiphysics 30 or AN-SYS Fluent. 37Guo and White 12 reported that a 2-dimensional, coupled, electrochemical-thermal model of a cylindrical required approximately 20 hours to simulate a constant discharge profile.In the same work, a proposed 1-dimensional, radial model of a spirally wound cell required 20 minutes to simulate the same discharge profile.Thus, these methods are not suitable for embedded BMSs, which require much lower computational workload.An approximate thermalelectrochemical model is therefore required for real-time applications.
Lower computational burden coupled thermal-electrochemical models have been reported in literature. 28,38For example, the lumped thermal model, which is based on the assumption that the whole interior temperature of the cell is uniform, 38 has been used in battery thermal management. 6This model is simple and requires low computational overhead and is useful in the thermal management of battery packs where individual cells within the pack are considered as simple heat sources.However, this approach is not able to compute nonuniform temperature distributions such as those that occur in large format cylindrical cells, especially in the radial direction. 39Kim et al. proposed a second order polynomial approximation to the governing radial temperature of the cell 40 in which the temperature at the center of the cell is a maximum.However, due to the hollow core in spirally wound cells the maximum will not occur at the center of the cell and in large format cylindrical cells this approximation is likely to be unreasonable.In addition, previous models 40 use only a single value of material density, specific heat and thermal conductivity for the whole cell and incorporate a fitting approach to estimate these parameters and incorporate a Kalman filter to estimate the cell's temperature.The computational requirements of such an approach are likely to be relatively high.
There are a number of SEI layer growth models available in the literature.Ramadass et al. 41 proposed a model of SEI layer growth representing a continuous film formation over the surface of the solid particles, which increases the surface resistance of the particles and increases irreversible capacity loss.Safari et al. 42 developed an isothermal, multimodal, aging model for SEI growth at the carbonaceous anode material considering the solvent decomposition kinetics and solvent diffusion through SEI layer.Henrik and Göran 43 proposed a degradation model based on linear combination of two phenomena, namely, continuous SEI layer growth at the microscopic scale and SEI layer cracking due to the expansion of solid particles at the macroscopic scale.
Up to this time, SEI layer growth is still considered as "the most important but least understood solid electrolyte in rechargeable Li batteries" due to the complexity of chemical and electrochemical reactions involved and its physical properties. 7,44For the purpose of control design, in this work, we adapt the SEI growth model proposed by Ramadass et al. 41 as it is simple whilst still being able to represent the underlying physics of degradation, and it has been well-accepted in the literature.We note, however, that other degradation models can also be used in our proposed model framework.
Coupled thermal, degradation and electrochemical models are uncommon, however, Tanim et al. 8 do present a SPM coupled with thermal and SEI models.The result is a 1D spatial model that ignores nonuniform distributions of temperature and degradation along the radius of the cell, which, given their dimensions, does not seem reasonable for large-format cylindrical cells.Smith et al., 4 present an empirical degradation model coupled with a multi-dimensional, multi-scale (MD-MS) cell model for large format cylindrical batteries.However, PDE-based models are computationally expensive and are used in battery design rather than real-time, optimal control situations.These models are too complex to be implemented in real time for most control algorithms. 45n this paper, we propose an efficient 1-dimensional, radial, coupled degradation-electrochemical-thermal model of a spirally wound, cylindrical lithium ion battery, noted hereafter as the DET model.The DET model couples our previous reduced order, electrochemical model 46 with an approximate thermal profile for the cell, based on the steady state solution of the one-dimensional, radial, heat equation, and an SEI growth model. 41The model accounts for the central, electrolyte filled, hollow core and the spirally wound material.The proposed DET model is computationally efficient whilst maintaining a high degree of accuracy (compared with that obtained when using the P2D model) in simulating the radial temperature and degradation distributions within the cell over time.It is a predictive modelling tool that can be applied to a wide variety of battery applications.The model is novel in that it applies to large-format, spirally wound, cylindrical batteries and accounts for non-uniform degradation and thermal behavior coupled to electrochemical phenomena.Our proposed model is able to provide insight into the variation and evolution of the local temperature and degradation rate of each individual wind along the cell radius.For the first time, this insight information is computed using the underlying physics of degradation, rather than the fitting of empirical models.
This paper is organized in the following way.In Model domains and reduced order electrochemical model in each wind section, we introduce the model structure and briefly describe the P2D model and our previous reduced order, electrochemical model.In Degradation model for each wind section, we develop the approximation for the radial temperature distribution in the cell.We then, in Thermal approximation model section, show how to couple this approximation to our reduced order model of each wind in the cell.Results are presented in Coupled degradation-electrochemical-thermal model section where we compare the proposed DET model to the full, one-dimensional, radial implementation of the P2D model implemented in COMSOL Multiphysics.

Model Domains and Reduced Order Electrochemical Model in Each Wind
A schematic of a cylindrical lithium ion cell and its cross-section is given in Fig. 1a.The cell has an outer radius R a and inner radius R 0 .The cell structure includes a centrally located hollow core, which is filled with electrolyte.The wound (or "jelly roll") structure of the cell consists of repeated electrode, separator and current collector layers 39 as shown in Fig. 1b.A single "wind" is also defined in Fig. 1b and consists of a sandwich of 8 layers beginning with the negative current collector and ending with the negative electrode.
The proposed model framework is multi-scale and includes a macroscopic domain, characterized by the radial variable, r (0 ≤ r ≤ R a ), in which the temperature distribution is solved and a microscopic  Temperature, T (r, t ), is assumed to vary along the radial coordinate, r, of the cell in such a way that the temperature within each wind (at the microscale) is uniform in space, however, the temperature across winds (at the macroscale) can vary.As the temperature in a wind is assumed to be uniform in space, and given that the applied current density is identically symmetric over each half domain of the wind, the electrochemical process occurring inside the three main domains in the first half of the wind (namely, the negative electrode, the separator and the positive electrode as shown in Fig. 1c) will be identical to the corresponding ones in three "mirrored" domains in the other half of the wind (namely, the positive electrode, the separator and the negative electrode).Therefore, we only need to solve for the electrochemical variables in one half of the wind as these values are then "mirrored" in the remaining half of the wind.
For each wind, we solve an electrochemical model and degradation model in sub-domain z (the microscopic domain) associated with the corresponding T (r, t ).Specifically, the input of each wind is the temperature T (r, t ), which affects the electrochemical processes including the transport of ions, reaction kinetics and SEI growth.The output of the electrochemical model are a series of heat generation rates, which in turn form the input of the thermal model in the domain r (macroscopic domain).
The electrochemical characteristics within each wind is solved using a Padé approximation based, reduced order, P2D model 46 for five state variables; including the electric potential ϕ s (z, t ) in the solid electrode, the electric potential ϕ e (z, t ) in the electrolyte, the lithium concentration of the active material c s (z, r s , t ) of the positive and negative electrodes, the lithium concentration c e (z, t ) in the electrolyte, and the molar fluxes j(z, t ) of the charge that flows between the active material in each electrode and electrolyte.In this model Padé approximations are used to simplify the complicated transcendental transfer functions that result from the Laplace transform solution of the linearized P2D model. 9,22We note that the Padé approximations are used to simplify all the equations of the linearized P2D model in Ref. 46 including ϕ s (z, t ), ϕ e (z, t ), c e (z, t ), and j(z, t ).For the lithium concentration in the solid phase c s (z, r s , t ), we only consider the transfer function of the lithium ion in the surface of the spherical particle c s (z, R s , t ).These transfer functions will compute the variation of the electrochemical variables in the domain z and the time t shown in Fig. 1c.Further details of the Padé model equations are given in Appendix E and in our previous work in Ref. 46.Corrections for nonlinear behavior, as given by Lee, Chemistruck and Plett, 9 are incorporated in the model.This approach reduces the complex nature of the above transfer functions to ones that contain only rational functions of simple polynomials in the Laplace domain. 24In this way we can develop a system of transfer functions that represent the variables of the P2D model for each wind.These transfer functions have the advantage in that they can be easily implemented and are computationally efficient.Here, they are used to calculate the voltage of each wind and the volume average heat generation rate in each of the domains shown in Fig. 1c.We discuss the coupling of these transfer functions with the temperature and degradation models in Coupled degradation-electrochemical-thermal model section, but first we introduce these models in the following sections.

Degradation Model for Each Wind
A SEI layer on the surface of the particles that form the negative electrode is one of the major causes of capacity loss and impedance rise for lithium ion battery. 7In this work we account for this by including model equations, based on those developed previously by Ramadass et al., 41 which describe the growth of the SEI layer during charging.
The irreversible side reaction, which forms a film at the solidelectrolyte interface of the negative electrode, is given by, S + 2Li + + 2e − → P, [1]   where, S is the solvent of the electrolyte and P is the "product" that forms the SEI layer.The ohmic resistance in the SEI film is related to the thickness of the film, δ film (z, t ), namely, where, and M p is the molecular weight of the SEI film, ρ p is the density of the film, κ p is the conductivity of the film.
Here, J film (z, t ) denotes the molar flux of the irreversible side Reaction 1 and is expressed as, where i 0,film is the exchange current density for the side Reaction 1, R is the universal gas constant and F is the Faraday's constant.
The overpotential, η film (z, t ), of the irreversible side Reaction 1 is calculated by, In this work, we solve this degradation model at the microscopic scale and couple it with the Padé model and the thermal model.Details of this coupling will be discussed in Coupled degradation-electrochemicalthermal model section.

Thermal Approximation Model
In this section we derive an approximation for the radial temperature distribution within the cell depicted in Fig. 1.We begin in Governing equation of the thermal model and heat generation terms sub-section, by recalling the governing equation of heat conduction, suitable for a 1-dimenasional radial domain, r, and we also briefly describe the heat generation terms in each microscopic domain, z.These generation terms are calculated using the electrochemical variables mentioned in Model domains and reduced order electrochemical model in each wind section.In Governing equation of the thermal model and heat generation terms sub-section we introduce an average thermal conductivity for the wound material within the cell.This conductivity is used to approximate the governing thermal equation as shown in Approximation of the thermal model sub-section.

Governing equation of the thermal model and heat generation
terms.-It is known that temperature variations occurring along the radius of cylindrical cells are caused by heat conduction in the radial and azimuthal directions. 38Chen et al. 47 reported that the heat flow in the azimuthal direction of cylindrical lithium ion cells that have a high number (20 or more) of winds can be negligible. 47They showed that when the number of winds is high, heat flow is predominantly radial and therefore a 1-dimensional, radial model can be used in preference to a 2-dimensional model to accurately represent the temperature distribution of such cells. 47In this work, we adopt that finding and use a 1-dimensional, radial, thermal model to calculate the temperature distribution within the cell.
The governing equation for the radial temperature distribution within the cell is given by, 48 ρc p ∂T (r, t ∂T (r, t ) ∂r + q(r, t ), [6]   where, ρ [kg/m 3 ] is the material density, c p [J/kgK] is the specific heat, λ[W/mK] is the thermal conductivity.Boundary conditions for (6) are shown in Fig. 2 as BC1, BC2 and BC3.The average heat generation rate per unit volume q(r, t ) [W/m 3 ] in each domain of each wind consists of reversible and irreversible heating, Joule heating in the solid and solution and resistive heating at the electrode/current collector interface.These heating terms are calculated using the electrochemical variables that have been solved using Padé model at the microscopic scale as described in the Model domains and reduced order electrochemical model in each wind section.
Irreversible heating, q i , is due to electrochemical chemical reactions and is given by, 6 q i = a s F j(z, t ) η(z, t ).[7]   where, a s is the specific interfacial area.Reversible heating, q r , results from a change in entropy, namely, 6 q r = a s F j(z, t )T ∂OCP ∂T .[8]   where OCP is the Open Circuit Potential of the of the electrode.( , ) ( , ) Joule heating in the solid, q s J , is a result of ohmic resistances and is given by, 6 q s J = where, σ eff = σε s is the effective solid conductivity, σ is the solid conductivity, ε s is the volume fraction of the solid phase, L i is the thickness of the electrode, and i ∈ {neg, pos}.Joule heating in the electrolyte, q e J , is ohmically related to the variation of electrolyte concentration and potential, namely, 6 ∂z ∂φ e (z, t ) ∂z .
[10] where, κ eff = κε Brug e is the effective electrolyte conductivity, κ is the electrolyte ionic conductivity, ε e is the electrolyte volume fraction, Brug is the Bruggeman coefficient, κ D eff is defined as + is the transference number, f ± is the mean molar activity coefficient, and i ∈ {neg, sep, pos}.
The heat generation rate per unit volume, qcollector , from current collector is, 6 qcollector (r i , t where, R collector is current collector resistance, A c is the cell surface area, and V c is the cell volume. Given the above equations the average heat generation term in each electrode is calculated as, qneg (r i , t ) = 1 0 q neg i + q neg r + q s,neg J + q e,neg J dz, [12]   qpos (r i , t ) = 1 0 q pos i + q pos r + q s,pos J + q e,pos J dz. [13]   whilst the average heat generation term in the separator is calculated only from the Joule heating in the electrolyte, namely, dz. [14]   The average heat generation rate per unit volume of a wind, q, in ( 6) is calculated by summing the products of the average (volumetric) heat generation in each domain ( qcollector for current collector, qpos for positive electrode, qneg for negative electrode or qsep for separator) with the volume of that particular domain and dividing this sum by the total volume of each wind.
Average thermal conductivity for the wound material.-Inorder to simplify (6) in the next sub-section, we first introduce an average thermal conductivity λ wound [W/mK] for the wound material (R 0 ≤ r ≤ R a ) as shown in Fig. 2. The average conductivity, λ wound , is expressed in (15) and is based on the concept of discrete thermal resistances (in each domain of Fig. 1c) in series. 49iven the temperature difference between the points at r i + r layer j and r i as shown in Fig. 2, the heat transfer rate can be expressed as: 49 q transfer i = T (r i ) − T r i + r layer j 2πH cell λ layer j ln r i + r layer j r i , [15]   where, H cell is the cell height, λ layer j is the thermal conductivity of each of the domains, r i is the radius of each wind and r layer j is the thickness of each of the domains as shown in Fig. 2.
Expression (15) has a similar form to Ohm's law in which the potential difference (T (r i ) − T (r i + r layer j )) across a conductor (2πH cell λ layer j /ln ) between two points (r i + r layer j and r i ) is directly proportional to the current (q transfer i ) through it.Therefore, the heat transfer problem in a composite cylinder, which consists of multiple domains, each with a different thermal conductivityλ layer j , can be represented as a system of thermal resistances connected in series as shown in Fig. 2. The thermal resistance of a single domain in the composite cylinder is given as: 49 The total resistance of the wound material, R wound , can then be expressed as: 49 Treating the whole wound material region as a single layer as shown in Fig. 2 allows us to define the average thermal conductivity λ wound as: We can therefore obtain: Approximation of the thermal model.-Inorder to solve the governing PDE in (6), some methods such as finite Fourier transforms 30 and spatial discretization methods 30,37 can be used; however, they require high computational workloads. 12For the feasibility of using DET model in control design, we seek for an approximation profile for (6), which is based on its steady state solution.
At steady state (∂T (r, t )/∂t ≈ 0) and assuming a uniform constant source term q, (6) becomes, which has the general solution, 49 Noting (21) we propose that an approximation of the temperature within the cell is given by, [22] where the initial temperature along the cell radius is equal to the ambient temperature, T amb , namely, T (r, 0) = T amb .[23]   We note that the approximation given in ( 23) satisfies a symmetry boundary condition at r = 0 (BC1), namely, The boundary condition at r = R a (BC2) describes the convective transfer of heat at the outer surface of the cylindrical battery, namely, where h [W/(m 2 K)] is the convection coefficient.Thermal continuity is assumed at the boundary r = R 0 (BC3 in Fig. 2), thus, and Substitution of the approximate form of T (r, t ) given in (22) into the boundary conditions (25) to (27) gives,

Ra
= hA(t ) [28]   Solving for B(t ), C(t ) and D(t ) in terms of A(t ), yields, Now, integrating both sides of ( 6 qdomain rdr.[30]   Substituting ( 22) and ( 29) into ( 30) and integrating we obtain, Equation 31 is the first order ordinary differential equation (ODE) for A(t ).A more detailed derivation of ( 31) and the form of the notations ξ , μ , and δ are given in Appendix A. Solving (31) for A(t ) then allows for B(t ), C(t ) and D(t ) to be determined from (29), which in turn determines the radial temperature distribution within the battery according to (22).Our approach considers the thermal properties such as material density, the specific heat and the thermal conductivity of different layer in the wounded material; however, it requires low computational workload since there is only one ODE function of A(t ) needed to be solved.
The parameters within the Padé model are temperature dependent.These parameters include the solid phase diffusion coefficient D s , the reaction rate κ, the electrolyte diffusion coefficient D e and the conductivity σ.These temperature dependencies can be calculated using an Arrhenius law, 6 namely, where, ψ represents one of the parameters D s , κ, D e , and σ, at temperature, T, and ψ ref is the reference value of that parameter at the reference temperature T ref . 37he temperature dependence in the open circuit potential OCP of the electrode is approximated by a first-order Taylor series expansion, namely, 48

OCP
where, OCP ref is the OCP of the electrode at the reference temperature T ref .

Coupled Degradation-Electrochemical-Thermal Model
In this section, we will show how to couple the thermal approximation to our Padé model and degradation model of each wind in the cell.To illustrate how the DET model is developed, we use the cell dimension of the large-format cell reported by Lee et al. in Ref. 39.The cylindrical cell has the inner radius R 0 = 4mm and the outer radius R a = 22.5mm.The parameters are obtained from Ref. 9 and are given in Appendix C. The radial geometry of the cylindrical battery consisting of 20 winds with 8 layers in each wind as illustrated in Fig. 1.As mentioned in Model domains and reduced order electrochemical model in each wind section, due to the assumption that temperature within one wind is uniform in space, electrochemical process in each wind is represented by one Padé model.Since the negative electrodes and positive electrodes of the 20 winds connect to the same negative current collector and positive current collector, respectively, we model this by connecting 20 Padé models in parallel as shown in Fig. 3.The output voltage and the input current of each Padé model must satisfy (34) according to Kirchhoff' circuit laws, namely, . [34] The coupling of the degradation, Padé and thermal models is shown schematically in Fig. 4. We see that the degradation model at each point z of each wind is coupled with the Padé model.The heat generation rates qcollector , qpos , qneg and qsep calculated from the electrochemical variables of the Padé model are used to simulate the temperature variation T (r, t ) at the macroscopic scale r.In turn the temperature T (r, t ) in each wind coupled with the Padé model and the degradation model.In our large format cell, consisting of 20 winds, there will be 41 coupled "sub-systems" of equations, 20 of which are Padé sub-systems, 20 are degradation sub-systems and one of which is a thermal sub-system as shown in Fig. 4.

Simulation and Results
In this section, demonstrate the performance of the DET model by comparing its outcomes to those of a full P2D model that has been modified to include degradation and thermal behavior.In the full P2D model, instead of approximating the PDE governing equations as presented in the Model domains and reduced order electrochemical model in each wind, and Thermal approximation model sections, we directly solve the coupled governing PDE equations using COMSOL Multiphysics. 30Specifically, we solve the P2D model coupled with the degradation model mentioned in Degradation model for each wind section to simulate the electrochemical variables on the microscopic scale.Simultaneously, we solve (6) which governs the temperature on the macroscopic scale.We use the same parameters that are used for the DET model simulation implemented in MATLAB/Simulink. 50These parameters are adapted from Ref. 9 and are listed in Appendix C.
Figures 5 and 6 show simulation results from the DET model, for the temperature distribution along the cell radius (here expressed in terms of the wind number) for different values of the convection coefficient h, plotted over time.Fig. 5 is for an input current of 1C for a single charge-rest-discharge-rest cycle and Figure 6 is for an input current associated with a UDDS (Urbane Dynamometer Driving Schedule) profile. 51Both of these current profiles are given explicitly in Refs.14,46.From each figure we observe that larger convection coefficient values h, result in more non-uniform temperature profiles, spatially.In Fig. 5 we can see that the temperature difference between the outer wind and the inner wind can be as high as 4°C.Alternatively, low values of h, result in higher cell temperatures over time, however, the temperature distribution along the cell radius is nearly uniform in these cases.
Figures 7 and 8 show the temperatures at r = R 0 and r = R a and the cell voltage over time resulting from the charge-rest-discharge-rest cycle and the UDDS cycle, respectively, for the DET and P2D models.Curves for different convection coefficient values are shown.We see

Wind 1 Wind 2
Wind n th I app,1 Padé model for wind 1

Padé model for wind 2
Padé model for wind n th ... ...  that in all of these cases the DET model results compares favorably to the full P2D model results.
The computation times for each of these applied current profiles is shown in Table I.All simulations were carried out on the same desktop computer with an Intel Core i7-6700 CPU running at 3.4GHz and 16GB RAM.We observe that the DET model requires approximately three minutes to simulate the 1C charge-rest-discharge-rest cycle profile in comparison to the reference model, which takes approximately eight minutes.The time savings of the DET model are more evident in the case of the UDDS profiles, with the DET model solving in approxi- mately 90 seconds, and the reference model taking approximately two hours.Fig. 9 shows the evolution of the degradation at the inner and outer winds for a periodic 1C charge-rest-discharge-rest current profile and h = 100 W /m 2 /K. Figure 9c shows the SEI resistance at r = R 0 and r = R a , which we can see increase over time.We also note that the SEI resistance of the inner wind is higher than that at the outer wind.This is due to the non-uniform temperature distribution across the cell.The higher temperatures at the inner radii facilitate the degradation reaction leading to thicker SEI layers and thus, higher resistances.Figure 9d shows the maximum lithium ion concentrations at r = R 0 and r = R a , which we can see decrease over time as lithium ions are removed by the degradation reaction.Furthermore, consistent with our observations for Fig. 9c, the maximum concentration in the inner wind decreases faster than that for the outer wind.Over extended charge cycles we observe that both the SEI film resistance and the maximum lithium concentration of the inner wind and outer wind diverge significantly.
Our DET model has the capability of predicting spatial and temporal temperature and degradation distributions for a spiral wound cell.Based on this information we can design optimal control algorithms that minimize, over time, distributions that degrade the cell's state-ofhealth.Examples of such algorithms have been given in the literature for situations where cell temperature is spatially uniform.Offline trajectory optimization is reported in Ref. 52     development of such algorithms for large format cells is beyond the scope of this paper but will be addressed in our future works.

Conclusions
This paper presents a new and computationally efficient onedimensional, radial, coupled degradation-electrochemical-thermal model of a large format, spirally wound cylindrical cell.We employ a Padé approximant model to compute the variation of electrochemical variables and heat generation terms in each wind and couple it with a degradation model.We derive an approximate model for the radial temperature distribution of the cell that is, in turn, coupled with the Padé and degradation models.The model is able to return accurate predictions when compared to those of a full P2D model, but in a fraction of the computation time.

Acknowlegment
This work was supported by the Australian Research Council Discovery grant DP160101325.The authors also acknowledge informative conversations with Professor Colin Please from the University of Oxford during the writing of this paper.
Here qdomain is the heat generation rate of each layer in each wind and is calculated in ( 11)-( 14     (A5) The overpotential of the intercalation reactions(i ∈ {neg, pos}): η i (z, t ) = ϕ i s (z, t ) − ϕ i e (z, t ) − U i OCP ( c i s,e (z, t )) − F R i film j i (z, t ).(A6) The potential difference between the positive and the negative current collectors yields the cell voltage: V cell (t ) = ϕ

Figure 1 .
Figure 1.(a) Illustration of a cylindrical lithium ion battery with spirally wound design and its cross-sectional view, (b) the component layers in each wind and (c) the domains that constitute half of each wind and the corresponding Padé model transfer functions.46

Figure 2 .
Figure 2. Illustration of cross section of the cylindrical lithium ion battery with boundary conditions.

Figure 4 .
Figure 4.A schematic representation of the DET model.

Figure 5 .
Figure 5. Simulation result of temperature distribution at different convection coefficients from DET model during a 1C charge-rest-discharge-rest cycle.

Figure 6 .
Figure 6.Simulation result of temperature distribution at different convection coefficients from the DET model during UDDS cycles.

Figure 7 .
Figure 7.Comparison of temperature variation between DET model and the full P2D model in COMSOL at r = R a (outer wind) and r = R 0 (inner wind) at different convection coefficients during a 1C charge-rest-discharge-rest cycle.

Figure 8 .
Figure 8.Comparison of temperature variation between DET model and the full P2D model in COMSOL at r = R a (outer wind) and r = R 0 (inner wind) at different convection coefficients during UDDS cycles.

Figure 9 .
Figure 9. Evolution SEI resistance R film and maximum concentration of lithium ion c s,max in solid phase in negative electrode of inner and outer winds (h = 100 W /m 2 /K).

Table I . Comparison of simulation time.
and online Nonlinear Model Predictive Control reported in Refs.53, 54 for these simpler cases.The