A Computationally Efﬁcient Multi-Scale Model for Lithium-Ion Cells

In this work, a computationally efﬁcient multi-scale and multi-dimensional model is set up to describe the electrochemical, electrical and thermal behavior for a generic pouch cell format. As solving the model in multiple spatial dimensions would require an extensive amount of computational resources, we apply effective spatial discretization techniques, namely the orthogonal collocation and Lobatto IIIA method. In order to reduce the number of electrochemical submodels, a coupling method based on node point interpolation is introduced. The proposed model shows an improvement in solution time by a factor of up to 60 while maintaining its accuracy compared to the ﬁnite element method solution. To investigate the spatial accuracy, simulation quantities such as potential distribution and temperature distribution for constant current discharge proﬁles are examined. With the aid of experimental data gained from Swagelok T-Cells, the model parameters are tuned in for discharge current rates of up to 10C and projected to a 40 Ah cell design. Due to the greatly reduced computational time, the proposed reformulated model can be used for complex physics-based simulations that are typically too computationally expensive with standard modeling approaches such as online estimation and parameter optimization. [DOI:

Lithium is viewed as one of the key materials for today and future applications in the battery sector. 1,2 Here lithium-ion cells (LIBs) and post lithium-ion technologies are seen as a promising technology to meet the growing demand in consumer electronics, electric vehicles (EVs) and energy storage systems. [3][4][5] In the EV sector, large-format cells have distinct advantages compared to smallformat cells due to their reduced share of inactive components and design simplifications of the pack and battery management system (BMS). 6 However, a large-format cell can suffer from negative effects arising from local changes of physical and chemical material properties which are induced by non-uniformities in the potential and temperature distribution across the electrodes. In addition, local gradients in temperature can cause non-uniform cell aging resulting in safety concerns for larger-cell formats. 7 These effects are often amplified by cooling conditions (e.g. surface cooling or tab cooling), 8 additional contact resistances 9 or geometrical properties such as the size of the current collectors and the placement of the cell tabs. [10][11][12][13] Various modeling approaches that take multi-scale and multidimensional (MSMD) electrochemical, electrical and/or thermal effects into account exist in literature. 9,[13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] MSMD models are usually structured in a modular model framework consisting of multiple submodels solved in separate computational domains at particle, electrode, and cell level. An additional advantage of this approach is that each submodel can be independently selected and adopted to the specific use case. At electrode level, most MSMD models are based on adoptions of the pseudo two-dimensional (p2D) porous electrode model developed by the Newman group [33][34][35][36][37] or simplifications such as the single particle model (SPM). 38,39 The p2D model resolves lithium-ion concentrations and potentials in the active material and electrolyte phase. 37 In contrast, SPM models neglect variations of lithium-ion concentration across the electrode as the active material is approximated by one representative particle resulting in a significant reduction of computation time. 38 Another popular class of MSMD models is based on an empirical approximation of the electrochemical cell behavior originally proposed by Tiedemann, Newman 40,41 and Gu. 42 Assuming a linear relationship between curz E-mail: stephan.kosch@tum.de rent and potential (linear polarization) under defined operating conditions, these type of models can be integrated in a framework with parallel electrodes 43 ensuring high accuracy and fast computational time for large-format cells. 44,45 However, like equivalent circuit models (ECM), 46 the model is only valid for a specific use case and the link to fundamental physical and chemical processes within the cell is lost.
Whereas MSMD modeling has proven to be a suitable tool to investigate non-uniform cell behavior, much effort has been made to reduce their computational cost while maintaining their ability to describe mass and charge transfer reactions. First MSMD models applied a resistive network approach to describe the potential and current distribution within the current collectors coupled to a 3D thermal model for planar 14 and cylindrical 20 cell designs. Between opposing nodes of the negative and positive current collector with the same coordinates, a p2D model calculates the current-potential relationship at the current collector-electrode interface to approximate the electrochemical behavior arising from mass transport, charge balance and reaction kinetics. The computation time can be significantly reduced by replacing the p2D model with a non-linear resistor fitted to the current-potential profile derived from the electrochemical p2D model. 14 Extensions of this approach with resistor and capacitor elements exist for spirally wound cells, including mapping techniques of the current collector domain between a two-dimensional (2D) and three-dimensional (3D) cell geometry representation. 22 Several works adopted the network approach and applied volume average techniques which can be used to reduce the number of submodels coupled to 2D/3D current collector and thermal models. 13,[15][16][17]21,26,30 In addition, coupling methods where each node point represents an electrochemical model were introduced for planar cell designs. 19,28,29,31 Guo et al. 19 compared the simulated temperature distribution of a node point coupling method and a reduced order volume average approach where one electrochemical p2D model is averaged over the computational domain. The simplified model coupling approach was extended to predict the potential and temperature distribution of a small battery module including the influence of the electrical busbars. 18 Further work included locally resolved 2D 9,23,25,27 or 3D 32 electrochemical models to describe mass and charge conservation in additional geometric directions. While a higher dimensional approach can be used to investigate local non-uniformities such as electrode edge effects of spirally wound cells, 25 the increase in number of degrees of freedom need a considerable amount of computational time. 25,27 MSMD models based on collocation method that include a 2D representation of the current collector geometry are rarely found in literature. The method is adapted for a 2D electrochemical model including the height of the current collector as an additional dimension (p3D). 27 In a more recent work, a 2D thermal model for a cylindrical cell design based on Galerkin collocation method has been published. 47 Both works show superior computational efficiency compared to standard numerical methods.
In this paper, we present a computational efficient MSMD modeling framework for large-format lithium-ion cells using model reformulation based on pseudo-spectral collocation and Lobatto IIIA method. We apply model reformulation for a coupled p2D electrochemical and 2D electrothermal model for a chosen planar cell geometry. At particle level we apply Lobatto IIIA method which has shown promising performance in p2D models. 48 In order to achieve high accuracy with minimal computational costs, we extend our previous work 26 by introducing a node point coupling method based on interpolation instead of volume average technique which can be used to reduce the number p2D models and enhance model preciseness. Hence, this work aims to foster the development of spatially resolved MSMD models that can be applied in fields such as online estimation or cell-design optimization where high accuracy and fast computation is required.

Model Reformulation
Although the porous electrode p2D model [33][34][35][36][37] can be seen as the reference for electrochemical models, a considerable amount of computational time is needed to handle the model equations which are necessary to describe mass and charge balance within the cell. In a p2D model, spatial discretization is performed in two dimensions: the (pseudo-) dimension "r" in the solid particles of the electrode, and the dimension "z" which is subdivided into a negative electrode domain, a separator domain and a positive electrode domain (Figure 1a). In both electrodes, the particle dimension "r" is coupled to the electrode dimension "z". The coupled equation system can be spatially discretized into a system of differential algebraic equations (DAEs) via finite difference method (FDM) or finite element method (FEM). This process involves the spatial discretization of each computational domain into a number of node points. Thus, the spatial discretization of typical p2D models lead to a high number of coupled DAEs with around 600 equations depending on the complexity and required accuracy of the problem. 49 Considering even a simplified 2D single layer cell representation with the additional dimensions "x" and "y" (Figure 1b), a coupled MSMD model consisting of several p2D models and spatially resolved potential and temperature distribution increases the computational burden significantly.
As shown by other researchers 27,50-57 , spectral collocation methods are a suitable solution to improve simulation efficiency in the field of battery modeling. In general, spectral methods can be classified by their basis functions which represent the expansion of the solution variable. The coefficients of the basis functions can be viewed as the spectral space representation of the approximated continuous function in a given interval. When a solution is sought over a finite number of points (collocation points) in the computational domain, the method is often called pseudo-spectral method or collocation method. 58 A major advantage of orthogonal polynomials (e.g. Chebyshev or Legendre polynomials) is that the collocation points cluster at the boundaries of the computational domain. Therefore the interpolation error near the boundaries, known as Runge's phenomenon, can be held to a minimum by properly choosing the orthogonal polynomial.
Orthogonal collocation was successfully applied for a p2D electrochemical and thermal model using the zeros of shifted Jacobi polynomials for cosine type basis functions 50 and later extended for a 2D electrochemical model (p3D) using Chebyshev polynomials. 27 Fur-ther works utilized the method in orthogonal collocation on finite elements (OCFE), 51 BMS applications [52][53][54][55] and optimal charging. 56,57 A summary of reformulation and order reduction techniques can be found in Jokar et al., 59 whereas a historical overview of pseudospectral methods and their application in chemistry and physics is given in Shizgal. 58 Orthogonal collocation.-Considering the Refs. 58, 60, 61, the standard collocation method based on polynomial interpolation can be described in the following way. To find an accurate solution of the original governing equation, the solution variable u(x) is approximated in the form of a truncated series expansion of (orthogonal) polynomial basis functions P k (x) [1] where N is the polynomial degree. Here C k is a set of expansion coefficients which can be determined by inserting u(x) into the PDE and forcing the residual of the resulting equation system to zero. A general class of basis functions P k (x) that minimizes the error of the approximation are Jacobi polynomials of which Chebyshev polynomials T k can be viewed as a special case. As Chebyshev polynomials are defined in the domain [−1,1], it is advantageous for battery models to shift the domain to [0,1] using the coordinate transformation The first few Chebyshev polynomials T k can be found in Trefethen. 61 Equations 1 and 2 lead to the generalized form of the solution variable u(x) and its p-th derivative: Applying the orthogonal collocation method, the solution is sought in the finite domain [0,1] using the zeros x i of shifted Chebyshev polynomials of the first kind which are known as Chebyshev-Gauss (CG) points. Another method of calculation is to solve T k (x) = 0 directly. In Equation 4 the subscript i indicates nodal space values. The CG collocation points x i form a set of node points (interior points) in the discrete domain 0 < x 1 < . . . < x N −1 < 1. In this work, the exterior points x 0 = 0 and x N = 1 are used to apply the boundary conditions. Equation 3 and Equation 4 can be utilized to transform the governing equation into a DAE system, which is solved by determining the expansion coefficients C k for the specified boundary conditions. However, a more convenient vector representation of the problem in form of the solution variable u(x) can be obtained by evaluating Equation 3 in nodal space values d pû dx p = D pû [5] whereû is the solution vector with x i (i = 0, 1, . . . , N ) rows and D p are the differentiation matrices of size (N + 1) × (N + 1). In consequence, the nodal space values x i now include the boundary points. The construction of the D p matrices is briefly explained in Appendix A. By reformulation of Equation 5 using the Kronecker tensor product ⊗ the system matrices for 2D grids in "x" and "y" direction can be derived from: In correspondence to the 1D case, N x and N y denotes the order of the polynomial which forms the corresponding horizontal and vertical block of grid points in a 2D domain. Consequently, the vectorû now consists of (N x + 1) × (N y + 1) rows and the solution in one point of the grid is dependent on the solution of the intersecting horizontal   Figure 1. a) Schematic of p2D Porous electrode geometry with particle dimension "r" b) Schematic of the single layer cell assembly consisting of current collectors, active materials and separator. A corresponding computational grid is presented below each submodel for clarification. and vertical block. The matrices I x and I y are the identity matrices of size (N x + 1) × (N x + 1) and (N y + 1) × (N y + 1), respectively. As the construction of the system matrices has to be done only for the 1D case to derive the 2D case, the implementation of the method is straightforward in software programs such as MATLAB. 62 We found that the introduced standard collocation method becomes unstable when a solid diffusion model with a concentration dependent diffusion coefficient is used. To approximate the highly transient concentration profiles, we apply the Lobatto IIIA method for the spatial discretization of the solid diffusion equation in each particle.
Lobatto IIIA.-As explained in Refs. 63, 64, Lobatto methods are widely used for the numerical integration of ordinary differential equations (ODEs). The Lobatto IIIA approach is a class of Runge-Kutta methods based on implicit trapezoidal rule in its second order form (s = 2). We shortly describe the method by considering the following first order differential equation subject to the initial condition u(r 0 ) = u 0 . Applying trapezoidal rule, a solution u is approximated in the discrete interval [r i , r i+1 ] at the node points i ≥ 0 by Here h is the distance between the two nodes points i and i + 1 defining each subinterval. If Equation 8 is applied for one node point r 1 we obtain 2 equations in the discrete domain r 0 < r 1 < r 2 with the corresponding solution u(r 0 ), u(r 1 ) and u(r 2 ). Considering one boundary condition, the resulting equation system is fully determined and can be solved by standard numerical techniques. Applying the Lobatto IIIA method for the spherical diffusion equation involves the reformulation of the original PDE into first order equations, which is explained in Appendix B.

Mathematical Model
In this work three different submodels solved in separate computational domains at particle, electrode, and cell level ( Figure 1) are employed. To describe the relevant transport mechanisms in the separator and electrode, several p2D models resolve lithium-ion concentrations and potentials in the solid phase (conductive electrode matrix) and liquid phase (electrolyte phase) between opposing nodes of the negative and positive current collector with the same coordinates. We consider charge balance in both current collectors, as the geometry of the current collectors and the tab positioning can result in additional polarization effects. The cell temperature can have a major effect on the transport mechanisms within the cell. Therefore, a 2D temperature model accounts for local changes of physical and chemical material properties due to non-uniform temperature distribution. As shown in Figure 1b, the current collector and temperature model form the 2D electrothermal model in the "x" and "y" dimension. Figure 1a, the computational domain in the "z" dimension consists of the negative electrode (index n), separator (index s), and positive electrode (index p). Each of these regions is rescaled to a domain of [0,1] and described by equations defining the liquid phase and the solid phase. The transformation into the dimensionless local coordinates is achieved by

Coordinate transformation.-As shown in
Therefore, the governing equations in each region of the porous electrode can be defined locally and the solution in whole computational domain can be described with the single dimensional variable "z". Although the submodels at cell level are described in a single region, the principles of the transformation stay the same For the sake of simplicity, we drop the subscripts related to the region of each submodel in following description of the model equations.
Porous electrode.-The electrochemical processes in the porous electrode are described by mass transport, charge balance and reaction kinetics by following the works related to Newman's porous electrode model. [33][34][35][36][37] The potential distribution ϕ l (z) and lithium-ion concentration c l (z) in the liquid phase is calculated in all areas of the domain, while the solid potential distribution ϕ s (z) is not computed in the separator area. Both phases are coupled via mass transport and conservation of charge at the particle surface, including reaction kinetics for intercalation electrodes. Thus, diffusion of lithium-ions in the solid phase c s (r, z) is limited to the particle-subdomain and concentration gradients in the electrode can only be compensated via the liquid phase.
Assuming a binary electrolyte in the liquid phase, the electrolyte solution consists of a lithium salt dissolved in organic carbonate solvents, where all species in the mixture have a mass and charge ratio of one in relation to each other. By following the principles of conservation of mass, the governing equation for diffusion-migration of lithium-ions in the electrode domain reads where ε l is the liquid phase volume fraction, D le f f is the concentration dependent liquid phase diffusivity and a s is the specific surface defined by the active area per unit volume of each particle. The subscript eff denotes an effective quantity which is scaled due to the porous electrode structure. Furthermore, it is assumed that the transference number of cations t + is constant in concentration and temperature. In Equation  10, the second term on the right-hand side represents the electrochemical reaction taking place at the interface of the solid/liquid phase in both electrodes. As there is no interface reaction in the separator, the pore wall flux at the particle surface j n becomes zero in this region. The liquid phase current density can be expressed as [11] where i l describes a flux of lithium-ions described by a modified Ohm's law using a quasi-electrostatic potential with respect to a lithium reference electrode. 33,36 Along with the diffusivity D le f f , the liquid phase conductivity κ le f f and the mean activity coefficient f ± are usually functions of concentration and temperature as pointed out in the Appendix C. Considering a single electrode reaction both equations are coupled through the pore wall flux j n which is zero in the separator region as pointed out. For a reversible reaction at the solid/liquid phase interface conservation of charge leads to ∂i l ∂z = − ∂i s ∂z = a s F j n [12] with the electronic current density i s i s = σ se f f ∂ s ∂z [13] and the solid phase electronic conductivity σ se f f of the conductive electrode matrix. Note that Equation 12 is only valid if one active material is assumed and one electron is exchanged in the reaction.
Consequently, the charge-transfer reaction at particle surface according to modified Butler-Volmer kinetics 65 reads RT η /den [14] where K r 0 is the reaction rate, η is the charge transfer overpotential and the factor 0.5 accounts for an equal anodic or cathodic contribution to the reaction. The denominator den is introduced to improve convergence and stability of the numerical solution when the local electrolyte concentration approaches zero. It can be expressed using an empirical coefficient c lim as In agreement with Mao et al., 65 we found that a value of c lim ≤ 1 mol m -3 can significantly improve numerical stabillity under high current rates with neglectable impact on model error. Although we prefer this modification due to its simplicity, such treatment may become unnecessary when intercalation kinetics are formulated based on thermodynamic theory. 66 The charge transfer overpotential is defined by the solid/liquid phase potential at the particle surface [16] and the open circuit potential U ocv . Due to the porous electrode structure, the transport parameters are commonly scaled by which is defined by the tortuosity τ and the porosity ε. The tortuosity can be described by the well-known relation τ = ε −α in which the Bruggeman coefficient α can be related to Bruggeman's work 67 and is often assumed to be α = 0.5. 35,37 However, these assumptions can underestimate the real electrode structure and the more generalized correction has shown better agreement with experimental data. 68,69 The prefactor f may be explained by individual contributions to the total tortuosity of the electrode composite (e.g. active material and conductive additive). 70,71 In accordance with experimental results, we estimated values of α = 0.5 and f = 2 for the NMC particles and α = 0.1 and f = 4 for the Graphite particles. These values are well in line with the reported data of comparable active materials and particle shapes. 68,69 The governing and boundary equations of the reformulated porous electrode model are summarized in Table I. Here the transition parameters of the boundary equations L ns and L sp arise from the scaling of the transport parameters between negative electrode, separator, and positive electrode. We explain the coupling with the current collector domain later in this chapter.
Particle.-In the solid particle-subdomain, the lithium-ion concentration in the active material particles is calculated as function of the mean particle radius in the additional dimension "r" (Figure 1a). Considering a time variant ion concentration in the spherical particles of the electrodes, the mass balance is given with [19] Here N s is the molar flux, D s is the solid phase diffusivity and c s is the concentration in the particles of the negative and positive electrode. It is assumed that the particles are initially at a uniform concentration c s (t = 0) = c s0 . The relevant boundary conditions at the center and surface of the particles are ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 129.187.254. 46 Downloaded on 2020-01-23 to IP Negative and positive electrode ( j = n, Boundary conditions negative electrode Boundary conditions positive electrode ∂ĉ ln Boundary conditions Transition parameterŝ where j n is the pore wall flux at the particle surface and R is the particle radius. As improved computational efficiency allows for more physical phenomena to be considered, we employ a solid diffusion model based on concentrated solution theory [72][73][74] in the dimensionless coordinatesr = r/R given by In Equation 22 the diffusivity D s0 is assumed to be concentration independent 72 and ξ = c s /c s,max denotes the mole fraction of lithiumions in the host material. The activity correction α(ξ), also referred to as thermodynamic factor, describes a non-ideal interaction between lithium-ions and the host material when the influence of expansion and contraction is neglected. 74 As is evident in Figure C2, the thermodynamic factor can be directly related to the open circuit potential U ocv . For the spatial discretization of the Equation 22 we adopted the transformation suggested by Urisanga et al. 48 based on Lobatto IIIA method. The reformulation of the described solid diffusion model is explained in the Appendix B.
Electrothermal model.-To calculate the potential and current distribution in the negative (index ccn) and positive (index ccp) current collector a two dimensional model based on charge balance is developed. As shown in Figure 1b, we consider a one layer representation which consist of half of a negative current collector, negative composite, separator, positive composite and half of a positive current collector. Therefore each repeating cell unit of the stacked pouch cell design consists of half of a current collector for both electrodes and the thickness is scaled by a value of 0.5 in the model. The potential distribution of the current collector is described by Poisson's equation Here cc is the local current collector potential and i n is the local transfer current density passing through the active material and separator ( Figure 1b). The thickness of the current collectors L cc and the electronic conductivity σ cc are assumed to be constant in the dimensions "x" and "y". The negative and positive current collector of the cell are made of copper (Cu) and aluminum (Al), respectively. It is assumed that the in-plane current is carried only by the current collector due to its high electronic conductivity compared to the electrode composite.
The boundary conditions at the negative and positive tab are n · −σ cc, p ∂ ccp ∂ y = i t, p = I app 0.5A t, p with A t, p = W t, p L cc, p · N cell [26] Here the first boundary condition sets the potential at the top of the negative tab to zero. The second boundary condition defines the current flowing through the tab by the applied cell current I app (A) divided by the cross sectional area of the positive tab A t, p . At the positive tab, the unit vector n is pointing in the inward direction for negative values of the applied cell current. At the remaining edges insulation boundary conditions are assumed. We summarized all boundary conditions for an exemplary 2D mesh in Figure 2 for clarification. The depth of discharge (DOD) of the full cell is dependent on the average lithium-ion concentration in the solid phase where ξ j,1 and ξ j,0 is the mole fraction of lithium-ions in each electrode for a fully charged and fully discharged cell, respectively (see Figure C2). The calculation and reformulation of c savr is briefly explained Appendix B. The governing equations for the heat balance can be derived by the principles of the conservation of energy [28] ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 129.187.254. 46 Downloaded on 2020-01-23 to IP where ρ cell is the density, c p,cell is the specific heat capacity at constant pressure, T is the temperature, and k cell is the thermal conductivity assumed to be equal along the "x" and "y" dimension. We account for additional heat generation due to contact resistances at the cell tabs by introducing the welding resistance R w with 44 where the volumetric scaling factors χ define the part of the cell volume where the heating occurs. 12 Note that R w represents a lumped thermal parameter which has to be adjusted to provide the best fit between modeling results and experimental data. The thermal quantities for the thermal model can be calculated by accounting for a series or parallel connection of thermal resistances 75 ρ cell c p,cell = j ρ j c p, j L j L cell [31] k cell = j k j L j L cell and k z = L cell j L j /k j [32] Here the subscript j indicates the individual component of the cell stack, L cell is the cell thickness and k z is the thermal conductivity in the "z" dimension. We estimated values of k z ≈ 1 W m -1 K -1 and L cell ≈ 6.7 mm for the cell used in this work. The calculation method of the cell thickness and all additional thermal properties can be found in our previous work. 12 As the thermal model is based on a 2D representation of the cell, we assume that the influence of cell internal temperature gradients are small compared to the surface temperature distribution. This assumption can be quantified by the Biot number Bi, which is an indicator of internal and external heat transfer resistance 76 To generate sufficient thermal gradients on the cell surface, we assume h = 25 W m -2 K -1 for heat flux due to convection which is in line with values used by other researchers to reproduce experimental data. 44,77 With Equation 33 the Biot number can be found around 0.07 implying that thermal uniformity condition holds true. However, for a cell with a bigger cell stack a multi-layer stack model combined with a temperature model that includes the "z" dimension may be considered. 50 The generated heat q gen in Equation 28 is given as a sum of irreversible and reversible heating, 78 ohmic heating in the current collectors and heat generated by contact resistance 79 q gen = χ s q pcm + χ n q ccn + χ p q ccp + q con [34] with q con = N cell L cell i n 2 R con [38] where q pcm is the total heat generated by the porous electrode model averaged through the layer thickness L layer . As the model is parameterized based on Swagelok T-cell data, we account for additional overpotential in a full cell design assuming a contact resistance R con located between the current collector and porous electrode matrix. 80 Note that the heat generated by contact resistance q con is transformed into a volumetric heat source via the number of cell layers N cell and the cell thickness L cell . 12 It is dependent on the local transfer current density, and therefore q con is varying in the "x" and "y" dimension. An easy way to solve the integral in Equation 35 in the collocation model is by symbolic integration of the continuous representation of the solution variable given with Equation A2 . The heat exchange with the cell's surroundings is implemented as purely convective heat loss at the surface of the cell [39] with the surrounding air of temperature T re f . Here the factor 2 accounts for the heat loss at both sides of the cell. The governing equations of the reformulated electrothermal model are listed in Table II.
Model coupling.-To effectively reduce the number of p2D models coupled with the current collector domain we introduce a node point coupling method based on interpolation instead of volume average technique. 21,26 Note that the following description assumes dimensionless coordinates as defined earlier in this chapter.
Consider 4 points forming a square in the element domain e with the edge coordinates ( Figure 3. In each edge point l = 1, . . . , 4 of the element e, the transfer current density can be associated with one p2D model.
) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 129.187.254. 46 Downloaded on 2020-01-23 to IP

Table II. Reformulated governing equations for the electrothermal model with (N x − 1) × (N y − 1) interior points.
Negative and positive current collector According to Equation 13, the coupling between the porous electrode and current collector domain is given by which implies that the solid phase current density at both current collectors are equal. To apply node point interpolation we can define the shape functions a e,l (X, Y ) with i n,e (X, Y ) = 4 l=1 a e,l (X, Y ) i n,e,l (X, Y ) [41] to derive the transfer current density in the element i n,e . The shape functions are defined to be 1 in one coordinate and 0 in all other coordinates of the element e.g. i n,e (X 1 , Y 1 ) = i n,e1 (X 1 , Y 1 ). Consequently, the transfer current density in the current collector domain can be derived from the composed solution of each element i n (X, Y ) = i n,e (X, Y ) (X, Y ) ∈ e , e = 1 . . . n e [42] where n e is the number of elements. Due to the improved spatial resolution of the transfer current density an increase in model accuracy compared to volume average techniques can be observed. The method is not bounded to rectangular elements and a variety of shape functions can be applied. To uniquely define the current collector potential we impose [43] in the positive current collector domain p . Following the model implementation described, the applied current I app assumes negative or positive values during a discharge or charge process, respectively. With Equation 43 and the imposed boundary conditions at the cell tabs charge balance in the current collector is ensured. Finally, the cell voltage can be estimated from the average voltage at the positive tab and the additional overpotential generated from the contact resistance R con by U cell =¯ ccp tab+ + R con I app W · H · N cell [44] To benchmark the reformulated model we investigate accuracy and solution time compared to the finite element solution in the software COMSOL 5.2a. The cell design is adapted from our previous work 12 with an overall nominal cell capacity of 40 Ah.
The simulation parameters values are listed in Table III and additional quantities are summarized in Appendix C.

Experimental
The electrodes used to derive parameters such as particle sizes, electrode compositions and porosities are graphite anodes and nickel-cobalt-manganese (NCM-111) cathodes. The NCM-111 (Li 1+x (Ni 1/3 Mn 1/3 Co 1/3 ) 1−x O 2 , HED-NCM-111) and electrolyte LP-57 (EC:EMC 3:7, 1M LiPF6) were received from BASF AG. The graphite was obtained from SGL Carbon GmbH, conductive carbon C-NERGY SUPER C65 from Imerys GmbH and polyvinylidene fluoride binder (PVdF, Kynar HSV 900) from Arkema, France. The N-methyl-pyrrolidone (NMP, anhydrous 99.5%) was purchased from Sigma-Aldrich. Electrode compositions of 95-0-5 (wt% Graphite-C65-PVDF) for the anodes and 96-2-2 (wt% NCM-C65-PVDF) for the cathodes were used. In case of the anode, 20 g graphite were mixed with 1.06 g PVDF (Thinky ARV-310). The solid content was at 55 wt% with mixing times of 15 min (20 g batch). The NCM-111 (30 g) was mixed in a 1-step process for 15 min with 0.625 g PVDF and 0.625 g C65 at a solid content of 65 wt%. The resulting inks were blade-coated using an automatic coater (RK Print, Germany) with gaps between 200 μm onto copper foil (10 μm, MTI) for the anode or aluminum foil (18 μm, MTI) for the cathode. The coatings were dried at ambient atmosphere either on a hot plate at 50 • C for at least 12 h or in a convection drying oven (also 50 • C, > 12 h). Subsequently 11 mm Ø (anode) and 10 mm Ø (cathode) electrodes were punched out. Compression was performed with a KBr press, for cathodes up to 2.5 t cm −2 (2 t for Ø 10 mm) and for anodes 0.8 t cm −2 (0.75 t for Ø 11 mm). The formation procedure for all produced cells consisted of three charge and discharge cycles at a C-rate of C/10. In the parameterization setup the cathode had an areal capacity of 2.53 mAh cm −2 (for 150 mAh g −1 ) and the anode of 3.15 mAh cm −2 . Full cells were assembled using Swagelok type T-cells with Lireference and glass fiber separator sheets (VWR, 250 μm compressed around 200 μm, 11 mm Ø, filled with 120 μl electrolyte LP-57) between working and counter electrode and one glass fiber separator (6 mm Ø, wetted with 20 μl electrolyte) between reference electrode and the stack.
The open circuit potential of NCM-111 was obtained in a half cell assembly vs. a lithium metal electrode with a loading of around 1.68 mAh cm −2 and an estimated porosity of 55%. The experimental data is averaged over two charge and discharge cycles at a current rate of 0.02C. For the graphite electrode, open circuit data is readily available in literature. 81 Experimental data was obtained by a battery test system from MACCOR at various constant current discharge rates (0.02C, 0.1C, 1C, 3C, 5C, 10C) at 25 • C ambient temperature. The device is capable of currents and voltages up to 5 A and 6 V with an overall accuracy of ±0.02% for the experiments carried out.

Results
Parameterization.-To test the validity of the estimated model parameters we compare the experimental data for the Swagelok T-cell described in the experimental part of this paper with the FEM solution of the p2D modeling approach. As shown in Figure 4, the model is able to follow the measured voltage profile at discharge current rates of up to 10C.
The steep voltage drop between 0.1C and 1C discharge cannot be described by the model without reducing the accuracy at higher discharge rates significantly. We reproduced the error for Swagelok Tcells with different electrode loadings. Although not the only possible explanation, this behavior may be related to the geometry of the Tcell design (11 mm Ø anode vs. 10 mm Ø cathode), which cannot be accurately described by a p2D modeling approach. However, as the focus of this work is the investigation of accuracy and computational time of the reformulated model, we found that the overall cell behavior can be adequately described. With the parameters tuned in for both electrodes, further simulations are performed to investigate the spatial accuracy of the reformulated model and FEM model.
Mesh size.-To evaluate the influence of the mesh size, simulations are performed on pouch cells with the tabs of the negative and positive current collector aligned symmetrically at the top of the cell. In the FEM model, the mesh size was refined until sufficient accuracy (error < 0.1%) was achieved. Compared to the reformulated model, the node point distribution is denser at the tabs. In the remaining area of the current collectors the location of the node points is identical for both models. The location of the node points in the reformulated model is dependent on the choice of the collocation points as explained earlier in this work. We found that the interpolation error at the tab discontinuities can be minimized by selecting the polynomial order such that the nodes in the "x" dimension are aligned with the tab locations. The two compared mesh sizes with identical tab location are presented in Figure 5.
The polynomial order used for the discretization of the geometry can be calculated by N x = n x − 1 with n x representing the number node points in the "x" direction (estimated in the same manner for the "y" direction). As a variety of basis functions can be used for proper node placement in the proposed MSMD modeling approach, the selection of the polynomial order can be done with low effort in standard software programs. For the Chebyshev polynomials used in this work we isolated the two mesh candidates which align with the tab locations in a matter of minutes. For each mesh the governing equation are imposed at the interior points and the associated boundary conditions are forced to be satisfied at the exterior points ( Figure 2).
We found that a number of node points in the "z" dimension for the negative electrode, separator and positive electrode of n n = 8, n s = 5 and n p = 7 is sufficient to approximate the relevant transport mechanisms in the cell. The resulting DAE system consists of 14 ordinary differential equations and 6 algebraic equations of which some can be eliminated. Note that the chosen discretization is in line with the findings of other research groups. 50,55,57 The number of node points in the negative and positive particle is set to n rn = 30 and n r p = 20 due to the arising non-linear equations of solid phase diffusion. Thus, a high number of nodes is necessary when a concentration dependent diffusion coefficient is employed. 48 To minimize the influence of the electrochemical model on the potential and temperature distribution, the number of p2D models is set to one for the comparison of the mesh size. In consequence, the transfer current density is identical for the reformulated and the FEM model with a value of 126.26 A m −2 for a 5C discharge. To evaluate the spatial accuracy in temperature and potential, we evaluate the peak model error at the end of discharge when the 2.7 V cutoff voltage is reached. Furthermore, a discharge rate of 5C is defined as the highest possible current rate of the projected 40 Ah cell. As the geometrical and material parameters are comparable to the cell design used in our previous work, 12 this discharge current can be seen as an upper bound for a manufactured cell of this format. Additional geometrical quantities describing the tab location can be found in Table CI. We emphasize that the governing equations of the FEM and reformulated model are identical but differ in their spatial discretization for the analysis carried out.
Mesh size.-Potential distribution.-The potential difference at the end of discharge in the negative and positive current collector is illustrated in Figure 6.
The potential imbalance is mainly caused due to the geometrical and electrical properties of the current collectors. In consequence, a faster change of potential can be observed in the tab region. As evident, the reformulated model is able to accurately predict the potential distribution in both current collectors. A discharge rate of 5C results in a potential drop of approximately 11 mV in the negative current collector ( Figure 6a) and 13 mV in positive current collector (Figure 6b). The largest deviation in potential compared to the FEM solution is located near the cell tabs as a result of interpolation error at the tab discontinuities ( Figure 7). For the coarser mesh the maximum potential error assumes values of roughly 0.7 mV and 0.8 mV for the negative and positive current collector, respectively. In case of the finer mesh an error of around 0.2 mV and 0.8 mV can be observed. The overall error in potential distribution throughout the whole domain is less than 5% for both mesh sizes, which is remarkable given that the cell tabs are approximated with a relatively low number of node points.  Figure 8, the temperature distribution is very similar for all investigated models.
Due to the higher resolution of the heat generation at the tabs, the reformulated model with the finer mesh resembles the FEM solution more accurately in this region. It is notable, that the introduced temperature imbalance of around 10 • C can be resolved by both mesh sizes.
The error in temperature is presented in Figure 9 for both mesh sizes. With a local maximum error of 1 • C, the approximation of temperature distribution with the coarser mesh is still adequate. The small temperature offset of 0.2 • C for both mesh sizes can be attributed to the combined error of ohmic heating in the current collectors and heat generated by the porous electrode model, which accumulates during the discharge of the cell. If the imposed thermal boundaries at the tabs are neglected, the temperature distribution of the reformulated and the FEM model is nearly identical. We identified the maximum error in heat generated by the porous electrode model with around 3% compared to the FEM model, implying that the discretization of the reformulated p2D model is sufficient to predict the overall cell behavior.
Mesh size.-Potential and temperature.-As the coarser mesh size is able to approximate the cell behavior with sufficient accuracy, we compare the tab voltage and surface temperature values for both models at different discharge rates ( Figure 10). As shown in Figure 10a, the reformulated model is able to follow the average voltage at the positive tab throughout the discharge of the cell.
The voltage error of both models at 3C and 5C discharge rate does not exceed 15 mV until the defined 2.7 V cutoff voltage is reached. At the end of the discharge, the error tends to be higher due to the steep change of the electrodes' open circuit potential. The small negative slope of the voltage error can be attributed to the Lobatto IIIA second order method used to reformulate the solid phase diffusion. As the approximation between the nodes is based on piecewise trapezoidal rule, the method underestimates or overestimates the concentration profile in the particles due to its mostly increasing or decreasing shape. Therefore, a higher order approach may be used when long continuous charge or discharge phases are simulated. However, we found that the improvement of accuracy by applying a fourth order method can be neglected for the applied discharge rates. The evolution of the maximum T max , minimum T min and mean T mean surface temperature and the error for the highest discharge rate of 5C are presented in Figure 10b. With a maximum error of 1 • C for all discharge rates, the simulated surface temperature values of both models are well in line. The sudden increase in maximum temperature error at the beginning of the 5C discharge is mostly due to interpolation error at the tab discontinuities as show in Figure 9. After the first phase of discharge the error levels off at 0.8 • C and remains almost constant. of p2D models (1x p2D, 9x p2D, 25x p2D) to the full distribution FEM model. The full distribution FEM model consists of 150 p2D models which are placed at every node point of the mesh shown in Figure 5a. With a maximum deviation in voltage error of around ±5 mV we observed no significant difference in tab voltage compared to the data shown in Figure 10a for all investigated placement cases. Consequently, we focused our analysis on the temperature behavior as presented in Figure 11. With a maximum error of −0.6 • C (T min ) and 1.2 • C (T max ) for a discharge rate of 5C, the reformulated 1x p2D model is still able to predict the temperature behavior compared to the full distribution FEM model (Figure 11a). For the 1x p2D, the rapid change in temperature error at the end of the discharge can be explained by the local change of the transfer current in the full distribution FEM model. The DOD imbalance causes a shift in local potential, resulting in a higher transfer current in bottom region of the current collector. In consequence, the increased heat generation reduces the imbalance in surface temperature, which can be observed for all distributed models in Figure 11. When a higher number of p2D models is employed the increase in spatial accuracy results in a more uniform error distribution throughout the discharge. The error in maximum temperature reduces from 1.2 • C (Figure 11b) to 0.9 • C (Figure 11c) whereas the error in minimum temperature stays at roughly 0.3 • C. The same trend can be observed when the comparison is based on FEM models. The results indicate that even a low number of p2D models may be sufficient to describe the cell behavior. As the employment of more p2D models increases the number the DAEs significantly, it is evident that a balance between computational burden and required accuracy has to be found. Note, that due to the applied node point coupling method, the p2D models can be freely distributed which can further reduce the number the equations.
Transfer current density and state of discharge.-The imbalance of the transfer current and DOD for the 5C discharge rate is shown in the Figure 12. We observed that the peak values of imbalance for both quantities arise in the endpoints of the center line of the positive tab for the investigated cell, which is where the near maximum values in potential and temperature difference can be observed (see Figure 6 and Figure 8).
During the first phase of discharge, both models show nearly equal results with a value of around 140 A m −2 for the transfer current density at the positive tab and approximately 120 A m −2 at the opposite site of the tab for an average transfer current density of i n,avr =126.26 A m −2 . In the reformulated model, the fluctuations in the beginning of discharge are related to the solid phase diffusion due to approximation error introduced by the concentration dependent diffusivity at high discharge rates. In the second stage of discharge, the transfer current is higher in regions which have been less discharged due to the shift in the local potential at roughly 30 Ah discharge capacity. The maximum DOD difference throughout the discharge can be observed in the positive electrode and is presented in the bottom of Figure 12. It indicates an error between both models of roughly 0.3 percentage points in this phase of discharge in the positive electrode. This error is primarily caused by the interpolation error of the reformulated model at the cell tabs. Note that the DOD difference in the negative electrode is around 1.5 percentage point lower whereas the error between both models is nearly neglectable. The results show that the spatial resolution of the reformulated model is sufficient to approximate the cell behavior with good accuracy.
Computation time.-All model equations are solved by means of the software COMSOL Multiphysics 5.2a on a laptop computer (i5, 2.5 GHz, 2 cores) with 8 Gb of RAM. Both models were solved by a direct method (MUMPS) and a backward differential formula solver (BDF) for time-stepping with a relative tolerance of 1 × 10 −4 . The solution time to solve both models for a constant current discharge of 1C is illustrated in Figure 13. To give a fair measure of the solution time, we closed all unnecessary tasks and averaged the results over 5 consecutive simulation runs.
The reformulated model shows an improvement in solution time of around 60 to 20 times compared to the FEM model dependent on the number of p2D models and the step time. As shown in Figure 13a, a higher increase in mean solution time for the reformulated model of factor 60 (1x p2D), 45 (9x p2D) to 16 (25x p2D) can be observed for a step time of 20 s. For a step time of 1 s, the same trend can be identified with values of 50 (1x p2D), 39 (9x p2D) and 21 (25x p2D). The results prove the superior performance of the introduced reformulated model approach. Note, that the model equations and underlying physics are the same for both models. For further comparison we listed the degrees of freedom and mean solution time for a step time of 1 s in Table IV.
The degrees of freedom (DOFs) are reduced by a factor of 9 compared to the FEM model. Each p2D model has approximately 3000 DOFs in the FEM model and around 500 DOFs in the reformulated model. A relatively high number of nodes in each solid particle is needed for both models to approximate the steep gradients in concentration caused by the concentration dependent diffusivities. For a larger number of DOFs a reduction in computation time difference can be observed. The simulations indicate a decrease in convergence with increasing DOFs. This may be related to the overall dense system matrices of the reformulated model due to the global nature of spectral collocation method. Furthermore, global convergence may be negatively affected by the introduced discontinuities at the tabs. However, we emphasize that in our analysis no specialized solvers are used that may be more efficient when dealing with dense system matrices. For a moderate number of DOFs the results of the proposed model are promising with possible applications in online estimation, parameter optimization and fast localized simulation.

Conclusions
In this work, we present a two-dimensional multi-scale and multidimensional model based on spectral collocation method which is capable of accurate approximation of cell behavior compared to the finite element solution. The reformulated model approach improved computation time by a factor of 20 to 60 dependent on the number of physics based submodels. Two mesh candidates for a planar large-format Li-ion cell are evaluated to determine the model error in the presence of discontinuities at the cell tabs. The analysis shows that even for a low number of node points a high spatial accuracy in  46 Downloaded on 2020-01-23 to IP temperature and potential distribution can be achieved. This allows for fast and accurate simulations under a variety of thermal and electrical boundary conditions, which can be used for cell design optimization to reduce adverse effects arising from non-uniform current density and state of charge throughout the electrode. Furthermore, we show that a low number of physics based models can be sufficient to describe the cell behavior compared to a full distribution model approach. The proposed coupling method based on node point interpolation allows for free distribution of physics based models between the current collectors to further reduce the computational burden. Since the model is able to account for additional non-idealities arising from mass transport in the solid phase, a high degree of flexibility when simulating different electrode materials is ensured.
Future work will involve the application of the presented modeling framework for cell design optimization, parameter estimation and the study of non-uniform aging behavior. b a Figure C1. Entropy of the a) negative 65  The temperature dependency of the solid phase diffusivity D s0 and the reaction rate K r 0 is assumed to follow the Arrhenius equation where T re f is the reference temperature and Ea is the activation energy. The corresponding parameter values are listed in Tables III and CI. The entropic coefficient of the reversible heat is derived from literature 65,83 and shown in Figure C1.
The open circuit potential and the thermodynamic factor of both active materials are defined as a function of the mole fraction of lithium-ions as shown in Figure C2. by using non-linear least squares curve-fitting of experimental data with an adjusted Rsquared value of around 0.99. It is noteworthy, that the resulting concentration dependent diffusivity is in accordance with literature data. 84 c savr Average Li+ concentration in particle (mol m -3 ) c pcell Cell heat capacity (J kg -1 K -1 ) D l Liquid phase diffusivity (m 2 s -1 ) D s Solid phase diffusivity (m 2 s -1 ) Ea

List of Symbols
Activation energy (J mol -1 ) f Bruggeman prefactor (-) f ± Mean activity coefficient (-) h Heat flux due to convection (W m -2 K -1 ) H Cell height (m) i l Liquid phase current density (A m -2 ) i s Solid phase current density (A m -2 ) i n Transfer current density (A m -2 ) i t Tab current density (A m -2 ) I app Applied cell current (A m) j n Pore wall flux at the particle surface (mol m -2 s -1 ) k cell In-plane thermal conductivity (W m -1 K -1 )