A multi scale multi domain model for large format lithium-ion batteries

A multi scale multi domain (MSMD) model for large format lithium-ion battery (LIB) cells is presented. In our approach the homogenization is performed on two scales (i) from the particulate electrodes to homogenized electrode materials using an extended Newman model and (ii) from individual cell layer materials to a homogenized battery material with anisotropic electrical and thermal transport properties. Both intertwined homogenizations are necessary for considering electrochemical-thermal details related to microstructural and material features of electrode and electrolyte layers at affordable comput ational cost s. Simulation results validate the MSMD model compared to the homogenized Newman model for isothermal cases. The strength of the MSMD model is demonstrated for non-isothermal conditions, namely for a 120 Ah cell discharged with four different cooling concepts: (i) without cooling (ii) with a base plate cooling (iii) with a tab cooling and (iv) with a side cooling. As one result, temperature gradients cause a local peak discharge up to 2.8 C for a global 2 C discharge rate. © 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The increasing demand for electric mobility results in the growing relevance of large-format battery cells for electric vehicles. In this case, electrode potentials and temperatures become heterogeneous at charging and discharging, as shown by Guo et al. [1] . Modeling these heterogeneities for a lithium-ion battery cell (LIB) from the macroscale, i.e. the surface temperature of the cell case, to the microscale, i.e. the state of charge (SOC) of a single active material particle, is a challenge.
The widely accepted electrochemical LIB model is the homogenized Newman model [2] , which has been derived first by volume averaging and then formally by two-scale expansion. Herein, the distinctive features of a porous electrode structure are not spatially resolved. The potential of the solid active material phase and the liquid electrolyte phase are assumed to be continuous functions of time and space coordinates. Thus, the simulation of ionic transport is simplified to a one-dimensional transport equation, while the characteristics of the porous electrode structure are described by effective transport parameters. The lithium diffusion in the active material phase is considered by an additional model calculating solid state diffusion in spherical particles. In this way, the * Corresponding author.
E-mail address: adrian.schmidt@kit.edu (A. Schmidt). homogenized Newman model is a pseudo two-dimensional (P2D) model, where one dimension is the electrode, and the other direction is the effective spherical particle. The P2D model is state-ofthe-art for calculating the electrochemistry in LIB batteries with a low computational effort [3][4][5][6][7][8][9] .
The Newman model considers transport phenomena in electrodes and electrode pairs only at the micrometer scale. Processes on larger length scales, such as heat or electronic transport in large format battery cells, is disregarded. A way out are multi scale multi domain (MSMD) models, as introduced by the battery modeling group at the National Renewable Energy Lab (NREL) about ten years ago [10] . MSMD models calculate the electrochemistry on the micrometer scale and map e.g. the temperature distribution in the entire battery cell on the mesoscale at the same time, by comprising separate solution domains for different length scales. Each domain uses its own independent system of the variables solved in that domain. Commonly, the solution domains are (i) the cell level for physical processes in the 10 cm scale, (ii) the electrode level for electrochemical processes in the 100 μm scale and (iii) the (active material) particle level for the solid-state diffusion in the 1 μm scale.
Two challenges remain for modeling large format cells: (i) a complex model geometry representing hundred(s) of electrochemically active layers (ii) a physicochemical submodel requiring a system of nonlinear partial differential equations, i.e., as the P2D model. Electrochemical heterogeneities, always present in large format cells, would require an indefensible computational effort, if solved by a matrix of P2D submodels. A way out is using either simplified model approaches or electrochemical submodels partly coupled to the MSMD framework [ 1 , 11-18 ]. Always, the chosen electrochemical submodel has to guarantee high-quality simulation results, especially at high C-rates. Evolving heterogeneities inside a large format battery cell demand for a locally evaluated and fully coupled electrochemical submodel.
There are commercial model frameworks as Batemo Cells [19] or Simcenter BDS [20] , aiming an accurate modeling of whole battery cells and may use multi scale model couplings. But since users only have restricted opportunities to insight, modify or extend the model equations, the scientific utility of these tools is limited.
Our MSMD model shall meet all of these above described demands as perfect as possible. The electrochemical submodel is using an development of the Newman model by Ender [21] , which handles the active particles as spherical ones, but extends the P2D model to a distribution of particle sizes. Thereby, the electrochemical coupling between different fractions of the particle size distribution is considered. The computational effort of our MSMD model is reduced by homogenizing the multilayer structure in large format cells. The detailed model description including field variables und exchange variables is given in the following chapter 2. The parameterization of our MSMD model to a large format cell (120 Ah) with graphite anode and NCA/LCO blend cathode is presented in chapter 3. The model's capability is demonstrated in the results section, where a comparison with the Newman model is made for isothermal und non-isothermal cases as well as a simulation study comparing four different cooling conditions: (i) without cooling (ii) with a base plate cooling (iii) with a tab cooling and (iv) with a side cooling.

Model description
The multilayer structure in large format cells, i.e., electrodes, electrolyte, separator and current collectors, is homogenized by adapting the approach of Newman and Tiedemann [2] . The schematic model structure is sketched in Fig. 1 . On the lefthand side, the homogenization of porous electrode structures to a pseudo-two-dimensional (P2D) model according to Newman is shown. Here ϕ correponds to the electric potential, c to the lithium concentration and i to the current density of solid phase (s), liquid phase (l) and the charge transfer (ct) on their interface. On the right hand side the stacked layer structure of a LIB is analogously homogenized resulting in a "homogenized cell material" with anisotropic transport parameters in the MSMD model (represented by the colored hatching). Here T represents the temperature, i ec the current density calculated by the electrochemical model and q the heat source, while the superscripts refer to the anode side (-) or the cathode side ( + ) Equivalent to the homogenization of porous electrode structures in the Newman model, the layered structure of the cell is treated as a homogenized superposition of the properties of all layers in the MSMD model. Thus, mathematically, each point of the cell consists simultaneously of positive and negative current collector, as well as porous cathode, separator and anode layers and the idea of the homogenized porous electrode of the Newman model is transferred to a homogenized cell material. Furthermore, following the Newman model, the transport properties of the materials contained in the homogenized phase are represented by effective transport parameters. However, a distinguishing feature is the directional dependence of these transport parameters. Since both, electrical and thermal transport parallel to the layers, differs strongly from the transport perpendicular to the layers, anisotropic effective transport parameters are introduced. In correspondence to the coupled diffusion model in the Newman model, a coupled electrochemical model considers the transport and exchange processes on the microscale.
Our model is subdivided into three length scales: the cell level, the electrode level and the particle level, as shown in Fig. 2 . At cell level, the field variables temperature T and current collector potentials ϕ + and ϕ − are calculated as well as the external boundary conditions (charge-/discharge currents and ambient heat transfer) are applied. These field variables on the cell level are averaged section wise and serve as boundary conditions for the electrode level. The electrode level is implemented through two extended homogenized models according to Ender [21] for cathode and anode respectively, representing an electrochemical submodel that can represent microstructural and material features of electrode and electrolyte layers including the impact of particle size distributions. It calculates the field variables of the potential in solid-and liquid phase ϕ s and ϕ l as well as the lithium concentration in the electrolyte c l . The resulting current density i ec and the generated heat q ec are calculated and returned to the cell level. At discrete node points the solution of the field variables is further transferred to the particle level, where the diffusion in the particles and the charge transfer current density i ct through their surfaces is solved. The implementation of the three submodels is presented in the following.

Electrical model
The electrical model of the cell level calculates the potential of the current collectors. Thereby, it is distinguished between the potential of the negative current collector ϕ − and the positive current collector ϕ + . For each of the two current collectors the charge transport is described by the ohmic law.
Where σ eff , + and σ eff , − are the direction-dependent effective conductivities of the associated current collector Eqs. (10) and ( (11) ) and i + respectively i − are the corresponding current densities. The current from negative to positive current collector i ec is calculated by the electrochemical model and serves as current source in the positive current collector and as a current sink in the negative current collector.

Thermal model
The temperature distribution within the cell level is calculated using the transient heat conduction equation.
Here, ρ eff is the effective volume-averaged density, c p , eff is the effective volume-averaged heat capacity, T is the temperature, and k eff is the effective thermal conductivity of the cell stack. The volume specific heat source q is composed of the ohmic heat in the current collectors and the heat source from the electrochemical processes q ec . The latter are calculated by the electrochemical model. q = i + · ∇ϕ + + i − · ∇ϕ − + q ec (6) Ohmic losses at the cell tabs cause a significant heat source Q tab [22] , which is approximated in the model by Eq. (7) .
Where I is the applied current at the tabs with the dimensions height h x width w x thickness d, σ tab represents the electrical conductivity of the tab material and σ ct is an equivalent conductivity which considers the electrical contact between the tab and the leading cable. The latter two values are taken from [22] .  Fig. 3 shows the layered structure and the corresponding coordinate system at cell level, wherein CC stands for current collector. For sake of simplicity a stacked layer structure is considered in this paper. The transport properties in x direction, perpendicular to the layers, differ strongly to the properties in y and z direction, planar to the layer structure. In the x direction the transport properties of the individual layers can be approximated as a series connection of resistors. Therefore, the effective electronic respectively thermal conductivity in this direction is described by Eq. (8) .

Effective transport parameters of electrical and thermal model
Where d is the thickness and k the individual conductivity of the anode layer (-), the cathode layer ( + ), the respective current collectors (subscript with cc) and the separator. In y and z direction the conductivity of the individual layers can be approximated as a parallel connection of resistors according to Eq. (9) .
The separator is assumed to be electronically insulating, conducting no electronic current transport through the layers (in the x-direction) which simplifies Eq. (8) for the case of electrical conductivity σ to: Within the layers, only the corresponding current collector and the respective active material contribute to the electric current transport, which simplifies Eq. (9) to: The effective heat capacity of the cell is calculated by multiplying the effective density ρ eff and the effective specific heat capacity c p , eff . The density is the mean value of the densities of the individual layers following Eq. (12) , wherein the densities of the porous layers are the mean value of the liquid phase with the volume fraction ε l and solid phase with the volume fraction ε s ( Eq. (13) ).
The effective specific heat capacity c p , eff considers the specific heat capacity c p , i , the density ρ i and the thickness d i of each layer following Eq. (14) . Eq. (15) applies to the specific heat capacity of the porous layers.

Boundary conditions
As indicated in Fig. 3 , the boundary conditions of the model are the applied current on the positive tab I tab , + and the applied potential at the negative tab, which is usually set to ϕ tab , − = 0 .
In addition, the initial and ambient temperature T ext of the cell is assigned and the heat transfer coefficient h therm according to Eq. (16) can be defined individually on each outer surface.

Electrode and particle level
The extended homogenized (P2D) model according to Ender [21] is the electrochemical submodel in our MSMD framework. In contrast to the Newman model, the particle size distribution of the spherical particles is considered herein. Thereby, the electrochemical coupling between different sized active material particles, often a cause of inhomogeneities in technical electrodes, is included. The relevance for discharge was shown by Ender [21] for anode structures and by Schmidt et al. [23] for cathode structures, but is omitted in this paper. The extended homogenized (P2D) model is depicted in Fig. 4 . The electrode pair level describes the electronic / ionic current flow and the electrolyte diffusion between the current collectors. The particle level considers the charge transfer reaction at the electrolyte / active material interface and the solidstate diffusion in the active material.
The electrode pair level is further subdivided into cathode d + , separator d sep , and anode d − domain. Therein, i s corresponds to the electronic current density in the solid phases of the electrode layers. It should be noted that there are two i s values for the positive and the negative electrode domains. As the separator is electronically insulating, i s is set to zero in this region. The ionic current density i l is considered in the separator as well as in the porous electrodes as they are all soaked with electrolyte. At the particle level the lithium diffusion is calculated for different particle sizes at every node of the discretized electrode domain. This results in particle size dependent surface concentrations c s , i and consequently different charge transfer current densities i ct , i . The specific surface-weighted (A spec ) sum of the individual charge transfer current densities gives the total charge transfer current i ct , which serves as the coupling condition to the electrode pair model. The detailed model description and the underlying equations are presented below.

Transport processes in the electrolyte phase
The processes in the electrolyte phase are described by the theory of concentrated binary electrolyte, as developed by Newman and Thomas-Alyea [24] . The conservation of mass in the electrolyte phase leads to: Where c l is the lithium concentration in the electrolyte, D l , eff is the effective diffusion coefficient in the electrolyte, F is the Faraday constant, t + is the transference number and i ct the charge transfer current density. The electrolyte current density is described by the charge balance in the electrolyte: Where κ eff is the effective ionic conductivity, φ l the potential and ∂ln f ∂ln c l is the thermodynamic factor of the electrolyte. R g is the universal gas constant.

Transport processes in the active material phase
The diffusion in the active material is based on the assumption of spherical particles. Thus, Fick's law can be formulated as a onedimensional problem by transformation into spherical coordinates as shown in Eq. (19) . This equation has to be solved separately for each of the N particles in the considered particle size distribution. Further, c s , i is the lithium concentration in particle i , and D s is the diffusion coefficient in the active material. The individual charge transfer current density i ct , i defines the flux of lithium-ions at the outer surface of each particle and determines the boundary condition for the diffusion equation according to Eq. (21) . Therein, the flux of lithium-ions must be scaled from the surface-volume ratio of a sphere ( 3 /R i ) to that of the active material (A spec , i /ε s , i ).
The electron transport in the active material is described by Ohm's law along the x-axis, where σ s , eff is the effective electronic conductivity in the solid phase: [a] experimentally determined value, temperature and SOC dependent.
[c] combined value as wetted separator

Charge transfer reaction
As coupling condition between active material and electrolyte phase the charge transfer kinetics according to Butler-Volmer is applied Eqs. (22) and ( (23) ). These are as well evaluated individually for each particle size. Here, k BV corresponds to the reaction rate constant, α to the charge transfer coefficient, η to the overpotential of charge transfer reaction, R SEI to the area specific resistance of solid electrolyte interface and φ ocv to the open circuit potential of the electrode. The charge transfer currents of the individual particles i ct , i are then summed up according to Eq. (24) to calculate the total charge transfer current density, which appears in the Eqs. (17) , (18) and (21) . Here, the active surfaces fraction of each particular particle size A spec , i / A spec is taken into account.

Effective transport parameters and temperature dependence
The properties of the porous microstructure are included in the model by effective transport parameters, considering the volume fraction of the respective phase ε and the lengthening of the transport path due to the tortuosity τ according to Eqs. (25) - (27) .
The temperature dependence of diffusion, charge transfer and SEI resistance are implemented by the Arrhenius Eq. (31) . The corresponding quantities X are calculated based on their value at 25 °C, where E act is the activation energy and k b is the Boltzmann constant.

Experimental
The model is parameterized to a prismatic large format (120 Ah) cell with graphite as anode material and an NCA/LCO blend as cathode material. All model parameters are given in Fig. 5 , Table 1 and Table 2 . Whereas microstructural parameters as well as most electrochemical and thermal parameters are measured by our group [ 21 , 27-29 , 34 ], further necessary parameters are from literature [ 25 , 26 , 30-33 ]. In the following, the measurement or selection of the model parameters will be discussed in more detail.

Cell level
In this section the experimental determination of the thermal transport properties according to Table 1 is presented. While Table 1 shows the thermal transport properties of the battery components for a state of charge (SOC) of 50 % and a temperature of  25 °C, these values have been measured and implemented in dependence of the temperature and the state of charge. The thermal conductivity of the materials is obtained by multiplying the density ρ, the specific heat capacity c p as well as the thermal diffusivity a according to Eq. (29) [35] . The devices and temperature ranges for the experimentally determined density and heat capacity parameters are listed in Table 3 . A more detailed description of the methodology for the determination of the specific heat capacity of porous electrode coatings is given by Loges et al. [36] .
For the challenging task of determining the thermal diffusivity and the effective thermal conductivity of porous electrodes a new experimental method was developed using a 467-HyperFlash from Netzsch in a temperature range of -20 °C -60 °C. Hereby, the electrode sample is thermally excited on the front side by a light pulse of a xenon flash lamp. The induced heat pulse penetrates the sample and an infrared sensor detects the temperature increase of the backside [ 37 , 38 ]. The average thermal conductivity is determined by evaluating the resulting time-dependent signal using a suitable model.
The commonly used adiabatic model of Parker et al. [37] assumes that the energy of the laser pulse is completely absorbed on the front of the sample without penetrating the material. This does not hold true for porous electrodes since the xenon flash partially penetrates surface and leads to a premature temperature increase on the sample backside. Therefore, the McMasters penetration model [38] is applied to determine the thermal diffusivity of the electrode stack, which considers the thickness of the absorption layer, predicts the premature temperature increase and is therefore suitable for porous electrode stacks. Liquid electrolytes evaporate even below room temperature. For the considered temperature range, helium ( ∼0,15 W (m K) −1 [39] ) has a similar thermal conductivity as the electrolyte (e.g. LP30: ∼0,18 W (m K) −1 [40] ) and is therefore found to be a suitable substitute filling fluid. Due to the negligible density of helium (0,16 kg m −3 [41] ) compared to the bulk material of the coating ( ρ s ρ He ), the effective thermal conductivity of the electrode sample ( Eq. 30 ) results by inserting the Eqs. (14) and (15) in Eq. (29) .
The total thermal resistance of the electrode sample can be approximated by a series connection of thermal resistances of the current collector and the electrode coating. Therefore, with knowledge of the thermal conductivity of the current collector material, Eq. (31) yields the effective thermal conductivity of the electrode coating.

Electrode and particle level
In the following the parameterization of the electrochemical model according to Table 2 is explained. The microstructure parameters (volume fractions, porosity, tortuosity, active surface area and particle size distributions) were obtained from focused ion beam (FIB) tomography [28] in case of the cathode and from X-ray tomography [42] in case of the anode.
The electronic conductivity of the active materials σ s was determined using the method described in [29] . The exchange coefficient k BV is determined by combining electrochemical impedance spectroscopy (hereinafter referred to as EIS) with an equivalent circuit model fit, as described by Costard et al. [43] . The impedance spectra of the individual electrodes are measured in an in-house developed cell housing and fitted to a physical motivated transmission line model. The fit enables the determination of the charge transfer resistance and, in the case of the anode, the SEI resistance of the electrode. With information about the active surface area of the electrode, the exchange current density and finally the exchange coefficient can be calculated. This method was applied in [27] for the modeled material system at different tem peratures and thus not only the exchange current density and the SEI resistance but also their activation energy was determined. The equilibrium potential in Fig. 5 b and c was measured in experimental cells of the type ECC-PAT-Core from EL-CELL in halfcell configuration against metallic lithium. Slow charge/discharge experiments with C/40 were performed four times in two individual experimental cells for each electrode. Thus, the reproducibility as well as the requirement for a quasi-static state of the electrodes was verified before averaging the curves to yield the equilibrium potential.
The parameters for a solution of 1 M LiPF 6 in EC:EMC (3:7 w:w) are taken from [44] . Therein, functions for the description of the ionic conductivity κ, the salt diffusivity D l , the thermodynamic factor ( 1 + ∂ ln f /∂ ln c l ) and the transference number t + are introduced as a function of temperature and salt concentration. These equations from [44] were implemented into the model without displaying them here once again.
In Fig. 6 literature values of lithium diffusion coefficients for NCA, LCO and graphite are summarized. The literature values for NCA and LCO range from 10 -7 cm ² s −1 [ 53 , 57 , 59 ] to 10 -12 cm ² s −1 [61][62][63] . Therefore, an average literature value of 10 -10 cm ² s −1 [ 55 , 58 ] is chosen for the simulations in this work. The literature values for the diffusion coefficient in graphite also vary between 10 -5 cm ² s −1 [46] and 10 -11 cm ² s −1 [47] . An average literature value of 10 -9 cm ² s −1 is assumed according to Nishizawa [ 48 , 52 ]. It should be noted that these are by far the most uncertain parameters applied in the model having a significant influence on the simulation results since low diffusion coefficients lowering the discharge capacity as shown in [21] . The comparative methodology of cooling concepts shown here is still applicable since all modeled cells use the same parameters. A more detailed discussion about the validity is given at the end of the results section.

Model implementation
Our MSMD model was implemented in COMSOL Multiphysics 5.3 using the Batteries & Fuel Cells module, the Heat Transfer module and the Transport of Diluted Species interface. The model setup is controlled by MATLAB via the COMSOL Livelink for MATLAB interface. Thereby the partitioning of the battery cell into sections, for each of which an electrochemical submodel is implemented (compare Fig. 2 ), is fully arbitrary. In this work the cell is divided into 4 × 4 × 4 sections. Thus, 64 electrochemical submodels are solved in parallel when calculating the model. At cell level we use an extruded free quad mesh with 635 vertices, while the electrode level is divided in 46 mesh elements, and each particle contains 10 mesh elements. The resulting 66,384 degrees of freedom are solved fully coupled using MUMPS (multifrontal massively paral-lel sparse direct solver) with a relative tolerance of 10 −3 and BDF (backward differentiation formula) free time stepping. The simulation of a complete discharge with this setup takes approximately 17 min and 5 GB RAM on a standard laptop with an Intel i7-8550U CPU (4 × 1.8 GHz). A subdivision into 3 × 3 × 3 electrochemical sections takes approximately 8 min and 3 GB RAM, while a model with 5 × 5 × 5 electrochemical submodels takes approximately 41 min and 9 GB RAM. For the shown cases, there was no significant improvement in accuracy between the 4 × 4 × 4 model and the 5 × 5 × 5 model. The resolution chosen in this paper therefore represents a good trade-off between computing time and electrochemical resolution. Generally, the model provides short computing times and can be executed on standard PC's and laptops. Furthermore, the complex interactions between the temperature distribution and the local electrochemistry can be modeled, which will be exemplified in two presented studies: A comparison of our MSMD model with the Newman type P2D model for model validation and a simulation study of four different cooling concepts for benchmarking.

Isothermal
The correct implementation of our model is validated versus the Newman (P2D) model by using the model parameters (cf. Table 2 ) for both. Since the Newman model does not provide thermal transport paths and thus is isothermal, the temperature of the "MSMD model isothermal" is set to the initial temperature T ext = cons tant = 25 °C . The calculated discharge curves are compared in Fig. 7 a showing excellent agreement between the MSMD model isothermal and the Newman model. Among them, the deviation for 1 C is as small as 4.4 mV from SOC 100% to SOC 0%, respectively 8.8 mV for 2 C. It originates from the ohmic drop in the current collectors, which is unconsidered in the Newman model.

Non-isothermal
The influence of the cell's self-heating is demonstrated by a non-isothermal MSMD model calculation under adiabatic boundary conditions.
The non-isothermal MSMD model shows significant deviations from the isothermal approach with an increasing voltage difference from -4.4 mV to 90.7 mV between SOC 100% and SOC 0%. (cf. Fig. 7 a and b). Naturally, the voltage differences further increase from -1 mV to over 200 mV at a discharge with 2 C. The underlying self-heating process elevates the averaged cell temperature from 25 °C to about 57 °C at 1 C, respectively 67 °C for 2 C (cf. Fig. 8 ). This goes hand-in-hand with decreasing internal resistances, which are considered in the non-isothermal MSMD model.
Diffusion, transport and transfer of lithium and lithium-ions are improved in the electrolyte, the active material and at the interfaces between them. This is connected to a lower overvoltage and thus a higher discharge capacity. The distinct influence of the thermally activated processes on the battery behavior and thus the relevance of non-isothermal modeling will become even clearer in the next section.

Simulation study of different cooling concepts
Four cooling scenarios are compared for a large format cell (120 Ah) using the non-isothermal MSMD model: (i) without cool-ing (ii) with a base plate cooling (iii) with a tab cooling and (iv) with a side cooling (cf. Fig. 8 a-d). Case (i) as the simplest scenario is only applicable for small charge and discharge rates. Case (ii) is a state-of-the art cooling concept which applies for numerous automotive applications because of its still simple and cost-effective realization. The cases (iii) and (iv) are rather complex in realization but are herein assessed to be more efficient and cause less thermal inhomogeneities.
For comparison, a complete discharge is simulated, starting at a homogeneous temperature of T = 25 • C with a discharge rate of 2 C. In the no cooling scenario, adiabatic boundary conditions are applied, assuming the cell is surrounded by other cells with the same temperature evolution. The other thermal boundary conditions are taken from Worwood et al. [64] : For base plate cooling and tab cooling, the dissipation heat is carried out by a liquid cooling system, assuming a heat pipe system with a water glycol mixture. The heat coefficient at the emphasized surfaces is 875 W m -2 K -1 . Side cooling is realized by an air-cooling system with an assumed heat transfer coefficient of 60 W m -2 K -1 . Fig. 8 e-h presents the local cell temperatures after a 2 C discharge with different cooling concepts (i) to (iv). Fig. 8 i-j shows the related average temperature progression as well as the temperature inhomogeneity during a 2 C discharge. The inhomogeneity of SOC, C-rate and temperature at a time t i is calculated according to Eq. (32) .
The no cooling case (i) results in the highest average temperature of up to 67.4 °C, with a rather small inhomogeneity of about 5 K. The maximum temperature is located at the positive tab (the right tab in the image).
The base cooling case (ii) results in the lowest average temperature being the only concept staying below 50 °C. At the same time, however, base cooling causes the largest inhomogeneity up to 25.6 K. The maximum temperature is located at the positive tab with up to 58.8 °C, while the lowest temperature is at the bottom with only 34.3 °C at the same time.
The tab cooling case (iii) results in an average temperature of up to 55.4 °C, caused by rather small cooling surfaces. The temperature inhomogeneities remain until SOC 50 % below 5 K. Thereafter, a pronounced increase up to 16.4 K is predicted (cf. Fig. 8 j). It becomes obvious, that the temperature increases with the distance from the tabs, the warmest area is at the bottom of the cell.
The side cooling case (iv) has a similar average temperature course as the base cooling with a maximum of 51.1 °C. The temperature inhomogeneity raises up to 15.5 K, which is lowest for all actively cooled concepts. The temperature gradient proceeds from the center to the outer cell surfaces, with the warmest point located at the positive tab.
The non-isothermal MSMD model impressively differentiates between the four cooling concepts. It is demonstrated, that calculating the cell's average temperature or measuring the cell's temperature at a single point is not sufficient, as the spatial distribution of the possible temperatures varies significantly.
As a next step, consequences of the raising temperatures on the electrochemical cell behavior are presented. Fig. 9 . shows the course of (a) cell voltage (b) SOC inhomogeneity and (c) C-rate inhomogeneity during a 2 C discharge at four characteristic time points t1-t4.
The course of the cell voltage (cf. Fig. 9 a) is only minor influenced by the chosen cooling concept.
The SOC inhomogeneities are shown in Fig. 9 b. Base cooling has the highest maximum with up to 8.5 %, followed by 4.2 % for side cooling, 2.9 % without cooling and 1.6 % for tab cooling.
The C-rate inhomogeneities are shown in Fig. 9 c. The courses of no cooling, base cooling and side cooling agree qualitatively, with three maxima located at SOC 50 % (t1), SOC 24 % (t2) and SOC 4 % (t4), and a minimum at SOC 13 % (t3). In contrast, the tab cooling differs by the absence of the peak at t1, and a slightly shifted maximum after t2.
Initially, the cell is in equilibrium state and the SOC in Fig. 10 b is uniformly at 100 %. At the same time t0, the C-rate in Fig. 10 a differs by 4 %, arising from the ohmic voltage losses in the current collectors: the voltage between positive and negative current collectors drops gradually with increasing distance from the tabs. Since the temperature at t0 is the same, the C-rate distribution holds true for all cooling concepts (cf. Fig. 9 c).  The faster discharge in the tab area until time t1 ( Fig. 10 c) is associated with a lower SOC ( Fig. 10 d). This is in turn connected with a locally lower equilibrium voltage, counteracting the elevated discharge rate and thus also the prevailing SOC inhomogeneity. Crucial for the counteracting process of the SOC inhomogeneities is therefore the derivation of the equilibrium voltage with the SOC, which is rather small at t1 (revealed by the flat voltage curve in Fig. 9 a). As a result, the C-rate inhomogeneities reach their first maximum and the SOC inhomogeneities increase further (cf. Fig. 9 b and c).
At time t2, the discharge curve (cf. Fig. 9 a) is much steeper, revealing larger equilibrium voltage differences induced by the fluctuating local SOC in Fig. 10 f. Thus, the counteracting process against the SOC inhomogeneities predominate, the SOC inhomogeneities decrease (cf. Fig. 9 b), resulting in a reversal of the C-rate inhomogeneities: Fig. 10 e now shows an excessed discharge rate at the base region of the cell.
At time t3 = 1564s, the discharge curve flattens briefly. The SOC inhomogeneities of the cell momentarily persist, resulting in a temporary homogeneous discharge (cf. Fig. 9 ).
From time t4, the cell voltage curve is very steep, revealing sudden large equilibrium voltage differences induced by the SOC inhomogeneities in Fig. 10 h. The resulting fast decay of the SOC inhomogeneities leads to a strong excess of the discharge rate up to 2.8 C in Fig. 10 e, respectively to the sharp peak at time t4 in Fig. 9 c.
The sharp peak in Fig. 9 c at time t4 is evident in all cooling concepts. Fig. 11 depicts the local distribution of the elevated discharge rates for the cooling concepts. The point of highest discharge rate correlates significantly with the point of lowest temperature of the cell in Fig. 8 e-h: For no cooling and base cooling at the bottom of the cell, for tab cooling at the tabs and for side cooling at the outer surfaces.
Furthermore, a quantitative correlation exists between the temperature inhomogeneity in Fig. 8 j and the electrochemical inhomogeneities in Fig. 9 . The base cooling has the highest inhomogeneities throughout the discharge, followed by the side cooling. The smallest thermal inhomogeneities of tab cooling in the first half of the discharge and the subsequent exceeding of the no cooling curve are likewise evident in the electrochemical inhomogeneities: In Fig. 9 b the tab cooling exceeds the no cooling curve at SOC 24 %; in Fig. 9 c the peak of tab cooling at time t2 is the smallest, while at time t4 it exceeds the no cooling curve. Moreover, isothermal simulations reveal SOC inhomogeneities below 0.5 % and C-rate inhomogeneities below 0.11 h −1 , which further underlines the strong temperature-dependent nature of the inhomogeneities.
While comparable previous model studies [ 1 , 12 , 13 , 18 , 22 ] focus on thermal inhomogeneities, the last part of our study impressively demonstrates how inhomogeneous cell conditions and local discharge rates are interlinked. Doubtless, this kind of timewise overloading is associated with local cell aging. This is experimentally confirmed by Werner et al. in [65] and [66] , wherein the cycle aging of consumer cells (3.2 Ah) is studied at different cooling concepts. In particular, temperature gradients strongly accelerate aging of capacity and polarization resistance.
The plausibility of our simulation results is underlined by further experimental results. After a cell discharge with 1 C at an ambient temperature of 25 °C, a temperature increase of 11.15 K and 21.45 K is reported in [67] and [68] at the surface of consumer cells (1.65 Ah pouch cell and 1.8 Ah cylidrical cell) without cooling. Our MSMD model predicts a temperature increase of 13,6 K for base Fig. 11. Local discharge rate distribution of the four cooling concepts during 2C discharge at t4 = 1728 s. cooling and 32.8 K without cooling, using the LIB cell parameters given in Tables 1 and 2 . Since the temperature increase depends very much on the selected cooling concept, also the heat generation is calculated. Here, Vaidyanathan et al. [69] reports a heat generation of 11.2 mW cm −3 at a 0.5 C discharge of a LCO/LNO cell, which is in good agreement with this work (9.35 mW cm −3 @ 0.5C discharge and a NCA/LCO cathode). Nevertheless, the experimental validation with a large format prismatic LIB cell of the same material composition is planned in future.

Conclusions
Our multi scale multi domain model (MSMD) for large sized lithium-ion battery cells applies separate solution domains for (i) the cell level, (ii) the electrode level and (iii) the particle level. We introduce novel homogenization approaches on two scales: (1) from the particulate electrodes to homogenized electrode materials using an extended Newman model and (2) from different material layers in the cell to a homogenized battery material with anisotropic electrical and thermal transport properties. In fact, the low RAM requirement makes it executable on standard laptops at affordable computational times.
Discharge characteristics are simulated using the following parameters: (a) cell geometry, dimensions and capacity, (b) thermal transport properties, (c) open circuit voltage of the NCA/LCO cathode, (d) microstructural and electrochemical parameters and (e) open circuit voltage of the graphite anode.
For the isothermal case, our MSMD model is in excellent agreement with discharge simulations made with the Newman model, but is superior for the non-isothermal case, as it considers selfheating effects.
This superiority is demonstrated for a 2 C discharge of a 120 Ah LIB cell, while applying four different cooling concepts: (i) without cooling (ii) with a base plate cooling (iii) with a tab cooling and (iv) with a side cooling. The arising temperature gradients are calculated and, i.e., the coupled SOC inhomogeneities and locally differing dischar ge rates are evaluated for all cooling concepts. For example, a local peak discharge rate of 2.8 C is proven for base cooling, potentially connected with excessive aging of the LIB cell.
In conclusion, our comparative study confirms, that induced inhomogeneous discharge rates originate from the interaction of (i) voltage losses in the current collectors (ii) emerging and decaying SOC inhomogeneities (iii) the slope of the discharge curve and (iv) the temperature profile of the cell. Even more important is the outcome, that the local distribution and the magnitude of excessive discharge rates develops counterintuitively. This underlines the necessity for fully coupled MSMD models.
The MSMD model presented is not only suitable for the identification of the best cooling concept for a specific cell design, but also for discovering the optimum cell design at given external boundary conditions and for defining safe operating conditions at locally arising overloads.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.