Modeling and Simulation of the Thermal Runaway Behavior of Cylindrical Li-Ion Cells—Computing of Critical Parameters

Abstract

The thermal behavior of Li-ion cells is an important safety issue and has to be known under varying thermal conditions. The main objective of this work is to gain a better understanding of the temperature increase within the cell considering different heat sources under specified working conditions. With respect to the governing physical parameters, the major aim is to find out under which thermal conditions a so called Thermal Runaway occurs. Therefore, a mathematical electrochemical-thermal model based on the Newman model has been extended with a simple combustion model from reaction kinetics including various types of heat sources assumed to be based on an Arrhenius law. This model was realized in COMSOL Multiphysics modeling software. First simulations were performed for a cylindrical 18650 cell with a $LiCo{O}_2$-cathode to calculate the temperature increase under two simple electric load profiles and to compute critical system parameters. It has been found that the critical cell temperature $T_{\text{crit}}$, above which a thermal runaway may occur is approximately $400\mathrm{K}$, which is near the starting temperature of the decomposition of the Solid-Electrolyte-Interface in the anode at $393.15\mathrm{K}$. Furthermore, it has been found that a thermal runaway can be described in three main stages.

1. Introduction

Lithium-ion batteries (LIB) have found a wide range of applications in many consumer products in the last 25 years. LIBs are used in cell phones, tablet and laptop computers as well as in electric (EV) or hybrid electric (HEV) driven vehicles. The size and the capacity of the batteries differs with the application. A battery used in an EV or HEV is much bigger than the one used for a cell phone. However the scale up of the batteries leads to a complexity of the physical phenomenona that do not play any significant role in single or small battery systems.

The occurrence of these physical phenomena are spread over a wide range of time and length scales, from the atomic level up to the heat transfer of the whole battery. The relation of these different phenomena to each other is of great importance, especially regarding the safety of the battery systems.

The usage of model-based investigations promises a theoretical understanding of battery physics even beyond the point that is possible using experiments. Due to their general structure, the behavior of batteries is strongly affected by the interaction of physics on varying length and time scales. In the last two decades, many approaches have been reported in literature starting with the work of Newman [1,2,3] based on the theory of porous electrodes [4]. This theory is based on a electrochemical description of diffusion dynamics and charge transfer kinetics of a battery and can give a forecast of the electric response of a cell in an intercalation electrode system. This model is adequate in small battery systems. In large formatted battery cells, however, the uniformity of the electrical potential along the current collectors in the cell composites is no longer given. Additionally due to the inhomogeneity of the distribution of the temperature field with respect to the cell geometry, thermal dynamics must also be taken into account.

As consequence an extended coupled electrochemical-thermal model must be considered. First thermal approaches for Lithium batteries are suggested by Newman and Pals [5,6,7] using an inhomogeneous heat equation , where the coupling of the thermal and electrochemical model is realized with the heat generation terms. This approach can be applied to single battery cells as well as to whole batteries [8,9,10,11,12]. In parallel the multi-scale and multi-domain character of the electrochemical model is developed [13,14] by using the mathematical theory of homogenization resp. the spatial averaging theorem [15]. Newer approaches are reported in the works of [16,17,18].

First, coupling between the multi-scale, multi-domain electrochemical-thermal model is introduced by [19,20,21] and is under development until now. A general short overview can be found in [22] and some current results in [23,24,25]. The purpose of this modeling approach is to gain a better understanding of the behavior of large LIB systems, because the transport of electrical current and heat must be considered on the length scales of the electrode particles up to the length scale of the cell composite, which is in general an anisotropic medium. One central point in the modeling of LIBs is the aspects of cell safety which corresponds to the thermal stability of a LIB. Several exothermic reactions can occur inside a cell during operation as its inner temperature increases. If the heat creation is larger than the dissipated heat in the surrounding space, this may cause heat to accumulate inside the cell and chemical reactions will be accelerated which yields in a further temperature increase until a thermal runaway is reached. In these terms, a thermal runaway describes a rapid temperature increase in a very short time interval of the typical order of $100\mathrm{K}/\min$ and above [26].

This phenomenon can not be described with the mathematical models mentioned above. A first attempt to overcome this disadvantage is given in [27]. The heat equation is coupled with ordinary differential equations (ODEs), describing the temporal evolution of the concentration of the exothermic reactions based on an Arrhenius-type law. Spotnitz et al. [28,29] give a first PDE modeling of a thermal runaway including reaction kinetics. In [30,31,32], the electrochemical-thermal model is extended with reaction kinetics based on an ODE formulation. Some newer results can be found in [33,34,35]. In [33] an oven test and the influence of convection is investigated, while [34] is focused on a Lithium-Titanate Battery with a model similar as the one described in [33]. The PhD Thesis of Tanaka [35] gives a detailed insight in the modeling and simulation of a thermal runaway for different chemistry.

In Section 2 of this article a model to describe the thermal behavior of a LIB in general with all sub-models is introduced and reviewed. In Section 2.1, the thermal model is formulated and the important heat sources are identified. In the following subsections, the corresponding mathematical models behind the several heat sources are given. In Section 2.2, the abusive exothermic kinetic reactions leading to a thermal runaway are given in terms of mathematical combustion theory [36,37] followed by a simplification in Section 2.3. In Section 2.4, the electrochemical heat source is considered and the multi-scale, multi-domain electrochemical model of LIB’s [22,23], which is implemented in the Battery and Fuel Cell Module in COMSOL Multiphysics, is introduced. In the Section 3 some results simulated in COMSOL Multiphysics are shown. In Section 3.1, the model with and without exothermic contributions is compared, followed by the computation of the time, where the thermal runaway occurs in Section 3.2. Additionally a first coarse classification of the thermal runaway is given. In Section 3.3 the critical overall cell temperature, where above the thermal runaway occurs, is determined. Furthermore intervals for the convective heat transfer coefficient and environmental temperature are computed in this case. In the final section, we summarize our results and an outlook of possible future work is given.

2. Modeling the Thermal Behavior of a Lithium-Ion Battery

Since the interior of a LIB is separated from the environmental air due to the battery can, the cell can be considered as a closed system. Under working conditions, i.e., the cell is exposed to an electrical load, heat is generated inside the cell due to several electrochemical and chemical processes. If more heat can be dissipated to the environment than is generated inside the cell, the LIB works in a stable state with respect to the temperature. If more heat is generated inside the cell than can be dissipated to the environment the cell is in an unstable state with respect to the temperature, leading to a permanent temperature increase, that finally can end in a thermal runaway. The critical aspect in the case of thermal runaway is its relation to energy conservation. With the help of energy conservation, one is able to describe thermal characteristics like heat generation and heat dissipation. For our consideration we state the following assumptions with respect to the geometry of the LIB (cylindrical, pouch or prismatic):

(1)

In a certain (infinitesimal) time interval the cell is under constant charge, discharge or in relaxation.

(2)

The cell geometry resp. volume is constant, i.e., elastic processes and gas generation inside the cell are not considered.

(3)

Without loss of generality (W.l.o.g.) the emissivity $\mathsf{ϵ}$, heat transfer coefficient h, the density $\mathsf{\rho }$, heat capacity $c_{\mathrm{p}}$ and thermal conductivity $\mathsf{\kappa }$ are assumed to be constant in space and time.

2.1. The Energy Conservation

Neglecting heat convection inside the LIB, the general equation for the conservation of energy can be derived from Fourier’s Law as initial-boundary problem [36,37,38] with the parabolic differential equation:

$${\displaystyle \mathsf{\rho }{c}_{\text{p}}\frac{\partial T}{\partial t}(\mathbf{x},t)=\mathsf{\kappa }\Delta T(\mathbf{x},t)+Q_{\text{gen}}(\mathbf{x},t)}$$

(1)

for the temperature $T:(\mathbf{x},t)\in \mathsf{\Omega }\times {\mathbb{R}}^{+}\mapsto \mathbb{R}$ of the cell and the corresponding initial and boundary conditions:

$$\begin{array}{cc}\hfill T(\mathbf{x},0)& =T_0\left(\mathbf{x}\right),\forall \mathbf{x}\in \overline{\mathsf{\Omega }},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \mathbf{n}·\left(\mathsf{\kappa }\nabla T\right)& =-h(T-T_{\text{env}})-\mathsf{ϵ}\mathsf{\sigma }(T^4-T_{\text{env}}^4),\forall \mathbf{x}\in \partial \mathsf{\Omega }.\hfill \end{array}$$

(3)

In this framework $\mathbf{x}\in \overline{\mathsf{\Omega }}\subset {\mathbb{R}}^3$ is a spatial interior or surface point of the cell and $t\in {\mathbb{R}}^{+}:=\{t\in \mathbb{R},t\ge 0\}$ is the time coordinate. $T_{\text{env}}$ denotes the environmental temperature, $T_0$ is the initial temperature profile inside the cell at $t=0\mathrm{s}$ and it is assumed that the initial temperature at the boundary of the cell coincides with the environmental temperature at $t=0\mathrm{s}$. $\mathbf{n}\in {\mathbb{R}}^3$ is the outward pointing normal vector. Ω denotes the interior of the battery cell, ∂Ω represent the boundary of the cell and the closure $\overline{\mathsf{\Omega }}=\mathsf{\Omega }\cup \partial \mathsf{\Omega }$ is the complete cell. Furthermore h is the heat transfer coefficient, ϵ the emissivity, $\mathsf{\sigma }$ the Stefan-Boltzmann constant, i.e., $\mathsf{\sigma }=5.670373\times {10}^{-8}\mathrm{W}/{\mathrm{m}}^2\times {\mathrm{K}}^4$, $\mathsf{\rho }$ the density of the cell, the heat capacity is $c_{\text{p}}$ and $\mathsf{\kappa }$ denotes the thermal conductivity.

In the inhomogeneity $Q_{\text{gen}}$ different heat sources are included. These contributions come from the heat generated by exothermic kinetic reactions, Joule and Ohmic heat and from reversible and irreversible thermodynamic effects:

$$Q_{\text{gen}}(\mathbf{x},t)=Q_{\text{el-chem}}(\mathbf{x},t)+Q_{\text{exotherm}}(\mathbf{x},t)$$

(4)

For each contribution of $Q_{\text{gen}},$ an additional complete mathematical model must be formulated.

2.2. Identifying Heat Sources From Exothermic Reaction Kinetics Inside the Lithium-Ion Batteries

Exothermic reaction kinetics are closely related to thermal abuse mechanisms inside a LIB. Several exothermic chemical reactions can occur inside a battery as the temperature increases. These may generate heat, which accumulates inside the cell and accelerates the chemical reaction between the cell components, if the heat transfer to the surroundings is not sufficient [39], i.e., if the heat generation rate exceeds the dissipation rate. External conditions for a temperature rise can be external heating, over-charging or -discharging, high current charging, nail penetration, external short or others. In these cases, a thermal runaway can occur resulting in leak, smoke, gas venting, flames etc., which leads to a destruction of the cells. Several authors have given models to describe abuse behavior and thermal runaway. First attempts to describe the exothermic reaction kinetics are given by [27]. A simplified reaction-diffusion model is given in the works of [30,31,32]. The case of nail penetration is described in [40]. A catastrophe theoretic approach can be found [41], which is based on a simplified ODE model and is extended in [42]. To describe the abuse behavior including thermal runaway, one has to identify the main exothermic chemical reactions. Following [26,27,31,32,43] the general mechanism that leads to a thermal runaway can be described with respect to rising temperature in four steps as follows:

(1)

SEI decomposition reaction: At $T>T_1$, in an exothermic reaction, the solid-electrolyte interface (SEI) decomposes.

(2)

Negative solvent reaction: At $T>T_2$ an exothermic reaction between the intercalated Li-ions and the electrolyte starts (NE).

(3)

Positive solvent reaction: For $T>T_3$ an exothermic reaction between the positive material and the electrolyte takes place under the evolution of oxygen inside the cell (PE).

(4)

Electrolyte decomposition: In a final exothermic reaction the electrolyte decomposes at $T>T_4$ (ELE).

The several temperatures may differ with the cell chemistry. For a $LiCo{O}_2$- chemistry the corresponding temperatures are [32]: $T_1\in [363.15-393.15]\mathrm{K},T_2>393.15\mathrm{K},T_3>443.15\mathrm{K}$ and $T_4>473.15\mathrm{K}$. Other exothermic- or endothermic contributions are neglected in this work. For example exothermic reactions caused by a phase transition of $Co{O}_2$ from hexagonal to cubic crystal structure give a reaction enthalpy of typical 4–$7\mathrm{J}/\mathrm{g}$ [44] in comparison to the SEI decompostion, where the enthalpy is between 250–$350\mathrm{J}/\mathrm{g}$ [28]. Furthermore the produced oxygen during the phase change can lead to additional exothermic reactions with the electrolyte which can accelerate the thermal runaway [45].

A suitable theory to describe exothermic reaction kinetics is the mathematical theory of combustion [36,37]. Within this theory, if one is mainly interested in temperature dynamics, the usage of a so called solid fuel model is a possible approach. A Solid Fuel Model consists of the Equations (1)–(3) and an additional system of parabolic differential equations to describe the time-spatial evolution of the concentrations, the several constituents of the exothermic kinetic reactions with the corresponding initial and boundary equations. These equations may be derived from Fick’s Law of diffusion. In the following, we assume that the several exothermic reactions are independent from each other and are governed by a simple Arrhenius law. With ${\mathsf{\Omega }}_i,i\in \mathcal{I}:=\{\text{SEI},\text{NE},\text{PE},\text{ELE}\}\subset \mathsf{\Omega }\subset {\mathbb{R}}^3$ and ${\bigcup }_i{\mathsf{\Omega }}_i\subset \mathsf{\Omega }$, we denote the spatial domain where the exothermic reaction took place, $\partial {\mathsf{\Omega }}_i$ is the corresponding boundary and ${\overline{\mathsf{\Omega }}}_i={\mathsf{\Omega }}_i\cup \partial {\mathsf{\Omega }}_i$ the closure. Then the i-th exothermic reaction can be written as:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial c_i}{\partial t}(\mathbf{x},t)}& =& d_i\Delta c_i(\mathbf{x},t)+{\mathsf{\gamma }}_i{R}_i,(\mathbf{x},t)\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill c_i(\mathbf{x},0)& =& c_{i,0}\left(\mathbf{x}\right),\forall \mathbf{x}\in {\mathsf{\Omega }}_i\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill c_i(\mathbf{x},t)& =& 0,\forall \mathbf{x}\in \partial {\mathsf{\Omega }}_i\hfill \end{array}$$

(7)

Here $c_{i,0}\left(\mathbf{x}\right)$ is the initial concentration, $d_i$ is the diffusion coefficient, ${\mathsf{\gamma }}_i$ the stoichiometric coefficient. Then the exothermic heat sources can be computed assuming an Arrhenius-law:

$$\begin{array}{ccc}\hfill Q_i(\mathbf{x},t)& =& q_i{R}_i(\mathbf{x},t)\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill R_i(\mathbf{x},t)& =& {\displaystyle A_i{c}_i^{m_i}(\mathbf{x},t)exp\left(-\frac{E_{a,i}}{RT(\mathbf{x},t)}\right)}\hfill \end{array}$$

(9)

where $q_i$ is the heat release in $\mathrm{J}·{\mathrm{g}}^{-1}$, $R_i$ the reaction rate in $1/\mathrm{s}$, $A_i$ the frequency factor in $1/\mathrm{s}$, $m_i$ the dimensionless reaction order, $E_{a,i}$ the activation energy in $\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$ and R the universal gas constant, i.e., $R=8.314462\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}·\mathrm{K}$.

This gives the total exothermic heat source as:

$${\displaystyle Q_{\text{exotherm}}(\mathbf{x},t)=\sum _{i\in \mathcal{I}}{Q}_i(\mathbf{x},t)}$$

(10)

2.3. Simplification of Combustion Model: The Constant Fuel Model

For our simulation the solid fuel model will be simplified to a constant fuel model [38,52] by neglecting all time-spatial dynamics of the concentrations $c_i$, if one assumes that:

$$\begin{array}{c}\hfill c_i(\mathbf{x},t)=c_{i,0},\forall t>0,x\in {\mathsf{\Omega }}_i\end{array}$$

(11)

$$\begin{array}{c}\hfill m_i=1\end{array}$$

(12)

This reduces the model to the Equations (1)–(3) and the exothermic heat sources defined in Equations (8) and (9) are only dependent on temperature T.

2.4. The Electrochemical Heat Source of a LIB and the Electrochemical Model

If one applies the thermal model described in last two sections to a LIB, this will coincide with a so called oven-trial for a Li-ion cell [31,32]. The cell is placed in a hot environment and the environmental temperature is increased with a specific profile. The cell is then heated only from the outside, if no current load is applied. If the environmental temperature is held constant and the cell is exposed to a current load, the cell will be heated from the inside. In addition to the exothermic heat contributions, additional electrochemical heat sources like Joule/Ohmic heating and reversible and irreversible effects will also occur. In summary these heat sources are given as

$$Q_{\text{el-chem}}=\sum _{j\in \mathcal{J}}{Q}_j,\mathcal{J}:=\{Joule,Irrev,Rev,Ohm\}$$

(13)

with

$$\begin{array}{cc}\hfill Q_{\text{Joule}}& =U·I\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill Q_{\text{Irrev}}& =I\left(U_{\text{eq}}-U\right)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill Q_{\text{Rev}}& {\displaystyle =IT\frac{d{U}_{\text{eq}}}{dT}}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill Q_{\text{Ohm}}& =I^2{R}_{\text{in}}\hfill \end{array}$$

(17)

where I is the applied current profile in A, U the cell voltage and $U_{\text{eq}}$ the equilibrium voltage of the cell in V. $R_{\text{in}}$ is the internal resistance of the LIB in Ω. The critical variable in the electrochemical heat source is the equilibrium voltage and derivation of the equilibrium voltage with respect to the temperature. In the Ohmic heat source, the internal resistance must be computed from the first partial derivative of the electric potential in the liquid and solid phase and the concentration of the Li-ions in the cell. These are additional variables to compute. These quantities must be computed with a separate mathematical model based on the porous electrode model developed by Newman [4]. This type of model has been under constant development in the last twenty years and describes a multi-scale multi-domain model (MSMD) to resolve the spatial-temporal electrochemical dynamics of a LIB with the Equations (1)–(3). The MSMD approach takes the physical and geometrical structure on different length scales and different geometrical domains into account. The general scheme is given in Figure 1. The three domains of interest and their inclusion in the LIB are identified (a) the particle domain; (b) the electrode domain and (c) the cell domain. To every domain an associated coordinate system taking into account the corresponding length scale and geometry is used.

At the electrode-electrolyte interface charge transfer kinetics have to be solved. The transport of the Li-ions is modeled with a diffusion mechanism and the migration and diffusion of the Li-ions through the liquid electrolyte can be evaluated. The charge balances in the solid and liquid phase are also resolved in the corresponding domains, i.e., cathode, anode and separator. The model is able to predict the electrical behavior of a LIB cell for given electrical load. In this extended form with the Equations (1)–(3) the temperature-distribution is solved over the whole cell geometry. The electrochemical dynamics of a LIB can mathematically be described as coupled system of elliptic and parabolic differential equations for the concentration of the Li-ions and the electric potential in the electrochemical active domain ${\mathsf{\Omega }}_{\text{el-chem}}\subset \mathsf{\Omega }$. From the concentractions of the Li-ions and electrical potentials in the subdomains the corresponding heat sources in the Equations (14)–(17) have to be computed with a suitable homogenization using the spatial averaging theorem [15] from current densities and heat fluxes. This model is implemented in the Battery and Fuel Cell Module of COMSOL Multiphysics. A short review can be found in [22] and a more detailed review in [46].

3. Simulations

First simulation results using COMSOL Multiphysics Version 4.4 are given in this chapter. The model was implemented in the Battery and Fuel Cell Module of COMSOL Multiphysics. This module includes some predefined models of a LIB as described in the last section. The user defines his simulation problem including the concrete models, their coupling, the cell geometry, the simulation parameters and the spatial-temporal discretisation.

The purpose of this first simulation is to get a better insight in the thermal behavior of a LIB under normal and abuse conditions. The standard battery model in this module is not able to resolve thermal runaway. In the following the standard model is labeled as Model A. Therefore, an extension of this model with additional contributions is required. This extension is described in the last section together with the battery model and implemented in COMSOL Multiphysics in our already existing model in the most simple case first. This model is labeled as Model B in the following.

We consider a cylindrical 18650 cell with $LiCo{O}_2$ chemistry, i.e., in the COMSOL Multiphysics material database a $L{i}_x{C}_6$ anode, a $L{i}_{x}Co{O}_2$ cathode and as electrolyte $1:1\text{EC}:\text{DEC}$ with a $LiP{F}_6$ salt is chosen. The internal spiral wounds are not resolved in this work, which means, that the geometric considerations will be restricted to the $r-z$-plane. Only the additional heat sources are used to extend the temperature equation with additional heat sources (Constant fuel assumption) from the reaction-diffusion system of the exothermic reactions inside the LIB. These heat contributions are directly implemented with the electrochemical heat sources in the MSMD approach.

It is assumed, that a single LIB is placed in an Accelerating Rate Calorimeter (ARC) and the environment of the LIB in the ARC is filled with air. Experimental measurements inside the ARC gave typical values for h in the range of $h\approx [0,30]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$. During the simulations the convective heat transfer coefficient h is varied in the interval of $h\in [0,24]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$. This setup represents adiabatic conditions for $h=0\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$. The non adiabatic conditions with positive h describes the free flow of gases and steams for $h\in [5,35]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$ in general, stationary gases and air in the environment are described by h in the interval $h\in [2,10]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$ [47]. Therefore only natural convection is considered in this work.

In a cylindrical geometry, the equations for the particle and electrode domain model remains unchanged. Only on cell domain the corresponding equations for the potentials at the current collector and the temperature field have to be reformulated in cylindrical coordinates in the $r-z-$plane, i.e., $\mathbf{X}=(r,z)$:

$$\begin{array}{cc}\hfill -\frac{1}{r}{\mathsf{ϵ}}_{\pm }{\mathsf{\sigma }}_{\pm }^{\text{eff}}\frac{\partial }{\partial r}\left(r\frac{\partial {\mathsf{\varphi }}_{\pm }}{\partial r}\right)-{\mathsf{ϵ}}_{\pm }{\mathsf{\sigma }}_{\pm }^{\text{eff}}\frac{{\partial }^2{\mathsf{\varphi }}_{\pm }}{\partial z^2}& =J_{\pm }\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \left({\mathsf{\rho }}_{\text{batt}}{c}_{\text{p},\text{batt}}T\right)}{\partial t}-\left(\frac{1}{r}\frac{\partial }{\partial r}\left({\mathsf{\kappa }}_{\text{T},\text{rad},\text{batt}}r\frac{\partial T}{\partial r}\right)+{\mathsf{\kappa }}_{\text{T},\text{z},\text{batt}}\frac{{\partial }^{2}T}{\partial z^2}\right)}& =Q_{\text{gen}}\hfill \end{array}$$

(19)

The initial and boundary conditions of the temperature field T are given in the Equations (2) and (3). Since the electrodes and the separator possess a porous structure, where the pores are filled with the electrolyte, the thermal parameter in the heat Equation (19) have to be considered as an effective parameter. Therefore in each of the three porous domains the total value of the thermal parameter is computed with a mixture formula containing the parameter value of the corresponding solid phase and the liquid electrolyte. Then the effective thermal conductivity (radial-, z- direction) ${\mathsf{\kappa }}_{\text{T},\{\text{rad},\text{z}\},\text{batt}}$, density ${\mathsf{\rho }}_{\text{batt}}$ and heat capacity $c_{\text{p},\text{batt}}$ are given as follows [30]:

$$\begin{array}{ccc}\hfill {\mathsf{\kappa }}_{\text{T},\text{z},\text{batt}}& =& {\displaystyle \frac{{\displaystyle \sum _{i\in \mathcal{I}}{\mathsf{\kappa }}_{T,i}{L}_i}}{L_{\text{batt}}}}\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\mathsf{\kappa }}_{\text{T},\text{rad},\text{batt}}& =& {\displaystyle \frac{L_{\text{batt}}}{{\displaystyle \sum _{i\in \mathcal{I}}\frac{L_i}{{\mathsf{\kappa }}_{T,i}}}}}\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\mathsf{\rho }}_{\text{batt}}& =& {\displaystyle \sum _{i\in \mathcal{I}}\frac{{\mathsf{\rho }}_i{L}_i}{L_{\text{batt}}}}\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill c_{\text{p},\text{batt}}& =& {\displaystyle \sum _{i\in \mathcal{I}}\frac{c_{p,i}{L}_i}{L_{\text{batt}}}}\hfill \end{array}$$

(23)

where $\mathcal{I}=\{pos.,neg.,c{c}_{\text{pos}},c{c}_{\text{neg}},sep\}$ denotes the index set to the corresponding parameter in cathode, anode, separator and current collector.

The last four equations can be regarded as special cases from corresponding formulas of the spiral geometry restricted to the $r-z-$plane [48,49]. Therefore, this thermal model can be considered as a serial thermal network in $r-$direction and a parallel thermal network in $z-$direction. With the parameters from

Table A1. Simulation parameters.

Initial values
Li-concentration neg. solid phase	$c_{\text{s},\text{neg},0}=7917\frac{\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}$	Li-concentration pos. solid phase	$c_{\text{s},\text{pos},0}=16002\frac{\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}$
Initial temperature	$T_0\in [273.15,373,15]\mathrm{K}$	Li-concentration liquid phase	$c_{\text{e},0}=2000\frac{\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}$
Li-diffusivity
Neg. solid phase	$D_{\text{s},\text{neg}}=3.9\times {10}^{-14}\frac{{\mathrm{m}}^2}{\mathrm{s}}$	Pos. solid phase	$D_{\text{s},\text{pos}}=1\times {10}^{-13}\frac{{\mathrm{m}}^2}{\mathrm{s}}$
Volume frcation
Neg. solid phase	${\mathsf{ϵ}}_{\text{s},\text{neg}}=0.471$	Pos. solid phase	${\mathsf{ϵ}}_{\text{s},\text{pos}}=0.297$
Neg. liquid phase	${\mathsf{ϵ}}_{\text{e},\text{neg}}=0.357$	Pos. liquid phase	${\mathsf{ϵ}}_{\text{e},\text{pos}}=0.444$
Maximal concentration
Neg. solid phase	$c_{\text{s},\text{max},\text{neg}}=26390\frac{\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}$	Pos. solid phase	$c_{\text{s},\text{max},\text{pos}}=22860\frac{\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}$
Thermal conductivity
Neg. electrode	${\mathsf{\kappa }}_{\text{T},\text{neg}}=1.04\frac{\mathrm{W}}{\mathrm{m}·\mathrm{K}}$	Pos. electrode	${\mathsf{\kappa }}_{\text{T},\text{pos}}=1.58\frac{\mathrm{W}}{\mathrm{m}·\mathrm{K}}$
Neg. current collector	${\mathsf{\kappa }}_{\text{T},\text{cc},\text{neg}}=298.15\frac{\mathrm{W}}{\mathrm{m}·\mathrm{K}}$	Pos. current collector	${\mathsf{\kappa }}_{\text{T},\text{cc},\text{pos}}=170\frac{\mathrm{W}}{\mathrm{m}·\mathrm{K}}$
Density
Neg. electrode	${\mathsf{\rho }}_{\text{neg}}=1347.33\frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}$	Pos. electrode	${\mathsf{\rho }}_{\text{pos}}=2328.5\frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}$
Neg. current collector	${\mathsf{\rho }}_{\text{neg},\text{cc}}=8933\frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}$	Pos. current collector	${\mathsf{\rho }}_{\text{pos},\text{cc}}=2770\frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}$
Particle radius
Neg. particle	$r_{\text{p},\text{neg}}=12.5\times {10}^{-6}\mathrm{m}$	Pos. particle	$r_{\text{p},\text{pos}}=8\times {10}^{-6}\mathrm{m}$
Reaction rate coefficient
Neg. electrode	$k_{\text{neg}}=2\times {10}^{-11}\frac{\mathrm{m}}{\mathrm{s}}$	Pos. electrode	$k_{\text{pos}}=2\times {10}^{-11}\frac{\mathrm{m}}{\mathrm{s}}$
Neg. electrode	$c_{\text{p},\text{neg}}=1437.4\frac{\mathrm{J}}{\mathrm{K}\mathrm{g}·\mathrm{K}}$	Pos. electrode	$c_{\text{p},\text{pos}}=1269.21\frac{\mathrm{J}}{\mathrm{K}\mathrm{g}·\mathrm{K}}$
Neg. current collector	$c_{\text{p},\text{cc},\text{neg}}=385\frac{\mathrm{J}}{\mathrm{K}\mathrm{g}·\mathrm{K}}$	Pos. current collector	$c_{\text{p},\text{cc},\text{pos}}=875\frac{\mathrm{J}}{\mathrm{K}\mathrm{g}·\mathrm{K}}$
Separator
Density	${\mathsf{\rho }}_{\text{sep}}=1008.98\frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}$	Thermal conductivity	${\mathsf{\kappa }}_{\text{T},\text{sep}}=0.344\frac{\mathrm{W}}{\mathrm{m}·\mathrm{K}}$
Heat capacity	$c_{\text{p},\text{sep}}=1978.16\frac{\mathrm{J}}{\mathrm{K}\mathrm{g}·\mathrm{K}}$
Battery geometry
Radius	$r_{\text{batt}}=9\mathrm{m}\mathrm{m}$	Height	$h_{\text{batt}}=65\mathrm{m}\mathrm{m}$
Mandrel radius	$r_{\text{mand}}=2\mathrm{m}\mathrm{m}$
Cell thickness	$L_{\text{batt}}=157\mathsf{\mu }\mathrm{m}$	Thickness battery canister	$d_{\text{can}}=0.25\mathrm{m}\mathrm{m}$
Thickness neg. current collector	$L_{\text{neg},\text{cc}}=7\mathsf{\mu }\mathrm{m}$	Thickness pos. current collector	$L_{\text{pos},\text{cc}}=10\mathsf{\mu }\mathrm{m}$
Length pos./neg. electrode	$L_{\text{pos}/\text{neg}}=55\mathsf{\mu }\mathrm{m}$	Length separator	$L_{\text{sep}}=30\mathsf{\mu }\mathrm{m}$

in the appendix, the corresponding effective parameters are given in

Table 1. Effective thermal parameter used in simulation in comparison with effective thermal parameter given in literature.

Parameter	Zhang [50]	Drake [51]	Simulation
${\mathsf{\kappa }}_{\text{T},\text{z},\text{batt}}$	Initial: $1.1\mathrm{W}/\mathrm{m}·\mathrm{K}$	$0.2\pm 0.01\mathrm{W}/\mathrm{m}·\mathrm{K}$	$0.898\mathrm{W}/\mathrm{m}·\mathrm{K}$
	Optimized: $0.48\mathrm{W}/\mathrm{m}·\mathrm{K}$
${\mathsf{\kappa }}_{\text{T},\text{rad},\text{batt}}$	Initial: $18.7\mathrm{W}/\mathrm{m}·\mathrm{K}$	$30.4\pm 1.5\mathrm{W}/\mathrm{m}·\mathrm{K}$	$27.544\mathrm{W}/\mathrm{m}·\mathrm{K}$
	Optimized: $21\mathrm{W}/\mathrm{m}·\mathrm{K}$
$c_{\text{p},\text{batt}}$	Initial: $913.60\mathrm{J}/\mathrm{k}\mathrm{g}·\mathrm{K}$	$1720\pm 86\mathrm{J}/\mathrm{k}\mathrm{g}·\mathrm{K}$	$1389.70\mathrm{J}/\mathrm{k}\mathrm{g}·\mathrm{K}$
	Optimized: $1243\mathrm{J}/\mathrm{k}\mathrm{g}·\mathrm{K}$
${\mathsf{\rho }}_{\text{batt}}$	-	-	$2172.99\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$

. In this table, the effective parameters are compared with parameters given in literature [50,51].

The main physical parameters of the simulation for both models can be found in Table A1 and in

Table A2. Exothermic simulation parameters.

$x\in \{\mathit{sei},\mathit{pe},\mathit{ne},e\},$	Reaction Heat	Frequency	Activation Energy	Volume Content
$y\in \{c,p,e\}$	${\displaystyle H_x/\left[\frac{\mathrm{J}}{\mathrm{kg}}\right]}$	factor ${\displaystyle A_x/\left[\frac{1}{\mathrm{s}}\right]}$	${\displaystyle E_A/\left[\frac{\mathrm{J}}{\mathrm{mol}}\right]}$	${\displaystyle W_y/\left[\frac{\mathrm{kg}}{{\mathrm{m}}^3}\right]}$
SEI reaction	$2.57\times {10}^5$	$1.667\times {10}^{15}$	$1.3508\times {10}^5$	$1.39\times {10}^3$
Neg. solvent reaction	$1.714\times {10}^6$	$2.5\times {10}^{13}$	$1.3508\times {10}^5$	$1.39\times {10}^3$
Pos. solvent reaction	$3.14\times {10}^5$	$6.667\times {10}^{13}$	$1.396\times {10}^5$	$1.3\times {10}^3$
Electrolyte decomp.	$1.55\times {10}^5$	$5.14\times {10}^{25}$	$2.74\times {10}^5$	$5\times {10}^2$

in the appendix. The additional simulation parameters of the exothermic reaction in Model B for the abuse contributions are taken from [30,31].

The simulations were performed with the two current load profiles with a constant current density of $i=\pm 80\mathrm{A}/{\mathrm{m}}^2$ as seen in Figure 2. In the first profile labeled Profile I shown in Figure 2a there are no relaxation periods between the charging and discharging pulses with duration of $250\mathrm{s}$ each. In the second profile (Profile II) relaxation periods of $250\mathrm{s}$ between charging and discharging pulses are used as shown in Figure 2b. A complete current cycle lasts $500\mathrm{s}$ and $1000\mathrm{s}$. The maximal simulation time is $t_{\text{sim}}=5\times {10}^4\mathrm{s}$ and $t_{\text{sim}}={10}^5\mathrm{s}$. For the time integration in COMSOL Multiphysics, a BDF integration scheme is chosen with a minimum order of 1 and a maximum order of 5 using a variable step size with a maximum time step of $1\mathrm{s}$ and an absolute tolerance of $0.001$. Since the spatial discretisation in COMSOL Multiphysics is based on the FEM method, we have used an adaptive spatial discretisation in the three models of the particle domain, the electrode domain and the cell domain respectively. The model of the particle domain is solved automatically in the battery module of COMSOL Multiphysics. Therefore only a spatial discretisation for the electrode domain and the cell domain is required. In the electrode domain, the maximum element size in the discretisation is chosen as $1\mu m.$ In total the discretisation contains 168 elements. Quadratic basis functions were chosen for the one-dimensional finite element discretisation in the electrode domain. In the cell domain the spatial discretisation is performed in the $r-z$-plane using 2266 triangular elements with the element size in the interval $[3.9\times {10}^{-4},8.45\times {10}^{-2}]\mathrm{m}$ and quadratic basis functions. The FEM discretisation is shown in Figure 3.

The main purpose of this section is to compute the overall mean cell temperature, below which the LIB shows normal behavior and above which the cell shows a thermal runaway with respect to some system parameters like the environmental temperature or the convective heat transfer coefficient. Furthermore, critical values for such parameters will be computed from the simulations and a first coarse classification of the thermal runaway will be given.

3.1. Model Comparison

In a first step Model A and Model B are compared under normal and abusive conditions with Profile (I). The following three cases are considered:

(1)

$h=0\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$: In this case the LIB is unable to dissipate any generated heat to the environment. Therefore, a thermal runaway will always occur after a sufficiently large time.

(2)

h positive and a thermal runaway occurs after a sufficiently large time: In this case, more heat is generated inside the LIB, as can be dissipated to the environment.

(3)

h positive and no thermal runaway occurs after a sufficiently large time: In this case the generated heat can be dissipated completely to the environment.

The first case coincides with adiabatic boundary conditions while the other two cases coincide with isoperibolic boundary conditions. Additionally for this comparison an environmental temperature of $T_{\text{env}}=293.15\mathrm{K}$ is chosen as a representative temperature in the interval of $[273.15,343.15]\mathrm{K}$. This temperature coincides with the minimum temperature where stable environmental conditions can be adjusted in the ARC. As corresponding convective heat transfer coefficients the values $h=0,5,10\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ are chosen. The spatial averaged mean cell temperature $\overline{T}$ of all interior points is considered W.l.o.g. the overall mean cell temperature $\overline{T}$ can be replaced by the mean temperature ${\overline{T}}_{\text{surf}}$ of the cell surface or the temperature $T(\mathbf{x},t)$ at an arbitrary interior point of the cell. , i.e.,

$${\displaystyle \overline{T}\left(t\right)=\frac{1}{\mathsf{\Omega }}{\int }_{\mathsf{\Omega }}T(\mathbf{x},t)\mathrm{d}\mathsf{\Omega }}$$

(24)

Since the simulation results are stored in a non-equidistant time grid, the mean temperature $\overline{T}$ is interpolated to an equidistant time grid with $\Delta t=0.1\mathrm{s}$ using splines with Matlab^® 2015b. In Figure 4a–c the Time- Temperature plots for Model A (blue) and Model B (red) are shown with respect to $h=0,5,10\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$. In Figure 4d–f the corresponding Temperature-Heating Rate plots are given on a semi-logarithmic scale. From Figure 4 one can see that in (a) and (b) both models shows similar behavior until $t\approx 1900\mathrm{s}$ resp. $t\approx 3000\mathrm{s}$. The corresponding cell temperature is approximately $400\mathrm{K}$ (Lower horizontal line in (a), (b), first vertical line in (d), (e)). This temperature is near to the start temperature of SEI decomposition, which is given in [26] to $393.15\mathrm{K}$. After this time Model B exhibits a greater increase in temperature as Model A. Finally, the simulation of Model B shows a blow-up in solution which corresponds to a thermal runaway in a LIB. For $T\approx 460\mathrm{K}$, the reaction for the electrolyte decomposition starts [26]. This is marked with the upper horizontal line in (a), (b) and the second vertical line in (d), (e) in Figure 4. The corresponding heating rate is $\approx 1.624\mathrm{K}/\mathrm{s}=97.44\mathrm{K}/\mathrm{m}\mathrm{i}\mathrm{n}$ which correspond very well with the start of the thermal runaway given in [26] with $100\mathrm{K}/\mathrm{m}\mathrm{i}\mathrm{n}$. In (c), both models have no significant difference in their time evolution of the temperature and the both solutions remain bounded. Therefore, in the graphs (a), (b) and (d), (e) of Figure 4 one can split the time evolution of the temperature in the three zones (I),(II) and (III).

First consider the plots (a) and (b) of Figure 4. In Zone (I) the cell shows a normal behavior. In this zone, the heat generation is governed by the electrochemical, Joule and Ohmic heat contributions. Exothermic contributions are negligible. In Zone (II), the cell temperature is high enough, so that the exothermic heat contributions are no longer negligible. This zone can be considered as an intermediate zone or the onset of the thermal runaway. In Zone (III), the exothermic heat contributions are dominant and thermal runaway takes place. Electrochemical heat contributions can be neglected in Zone (III).

Next consider in Figure 4 in the plots (d), (e), (f) the corresponding Temperature-Heating Rate trajectory. Since the current profiles Profile I and Profile II are of a rectangular shape, the first derivative $\text{d}\overline{T}/\text{d}t$ will be a piecewise continuous function with discontinuities at the time, where the current profiles jumps. Therefore, we consider $\text{d}\overline{T}/\text{d}t$ in the sense of distributions. At the time of the discontinuities one has to take the limit from right and left to get a derivation from left and right.

The graphs in (d) and (e) can be separated into three zones again (Vertical black lines). In the first zone (I) Model A and Model B show similar behavior with an increase in the mean cell temperature $\overline{T}$ while the heating rate $\text{d}\overline{T}/\text{d}t$ remains almost constant or has a small decrease with respect to some outlier points. These points coincide with the time, where the current profile jumps from the charging state to the discharging state and vice versa. At the end of zone (I) in Model B, the heating rate begins to rise. This is justified due to heat contribution of the exothermic reaction, which is no longer negligible. The general shape of the curve from Model B remains unchanged. This is the transition to zone (II). In the zones (II) and (III) one observes an increase in temperature T and an increase in the heating rate $\text{d}\overline{T}/\text{d}t$, which is justified in the accelerating of the exothermic reactions. The transition from zone (II) to (III) is reached, when the heating rate exceeds the critical value of $\text{d}\overline{T}/\text{d}t\ge 1.67\mathrm{K}/\mathrm{s}$ [26]. The corresponding temperatures are computed from the plots (d) and (e) in Figure 4 as the intersection the horizontal black line with the graph of Model B (red). This result is compared with the upper black line in the plots (a) and (b) with the graph of Model B (red). The results are given in the last column in

Table 2. Classification of the time evolution of the temperature during thermal runaway.

Zone	Abraham [26]	Peng [32]	Simulation
(I)	$T<393\mathrm{K}$	$T<393\mathrm{K}$	$T<400\mathrm{K}$
(II)	$T\in (393,453)\mathrm{K}$	$T\in (393,473)\mathrm{K}$	$T\in (400,464)\mathrm{K}$ for $h=0\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$
			$T\in (400,461)\mathrm{K}$ for $h=4\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$
(III)	$T>453\mathrm{K}$	$T>473\mathrm{K}$	$T>464\mathrm{K}$ for $h=0\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$
			$T>461\mathrm{K}$ for $h=4\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$
Chemistry	$LiN{i}_{0.8}C{o}_{0.15}A{l}_{0.05}{O}_2$	$LiCo{O}_2$	$LiCo{O}_2$

. Next, we compare our simulation results for Model B with the experimental results from Abraham et al. [26]. In their work, the thermal abuse of a 18650 cell with a $LiN{i}_{0.8}C{o}_{0.15}A{l}_{0.05}{O}_2$ chemistry is considered. The time evolution of the temperature is also grouped into three zones. A comparison between the experimental results from [26], the results given in [32] and our simulation results is given in Table 2. From a physical point of view one can interpret the zones (I)–(III) as follows: If the $\overline{T}-\text{d}\overline{T}/\text{d}t$ curve remains in zone (I) and stays below certain temperature values the cell works in safe conditions. If the $\overline{T}-\text{d}\overline{T}/\text{d}t$-curve changes to zone (II) thermal runaway starts, and both the heating rate and temperature rise. In zone (III) the heating rate exceeds a critical heating rate value. This is the thermal runaway, in which the heating rate tends to infinity. In Figure 4c,f both models stay in zone (I), no thermal runaway occurs and no significant difference is to be observed in both models.

In summary one can see the difference between Model A and Model B due to the exothermic heat sources in the temperature development over the time and corresponding heating rate. Furthermore Model B gives a good agreement with experimental results from literature [26,32] and a coarse classification of the time evolution of a thermal runaway event has been made. This topic will be considered in more detail in a subsequent publication.

3.2. Computing the Thermal Runaway Time $t_{TR}$

In this and the following section only Model B is considered. The environmental temperature is $T_{\text{env}}\in [273.15,373.15]\mathrm{K}$ and steps of $\Delta T_{\text{env}}=20\mathrm{K}$ were employed. In the simulations for each $T_{\text{env}}$, the heat transfer coefficient h is in the interval $[0,10]\mathrm{W}/({\mathrm{m}}^2·\mathrm{K})$ and the emissivity is equal to $0.02$. The time required for thermal runaway to occur with respect to the environmental temperature and convective heat transfer coefficient is computed in the simulation. For this purpose, the term Thermal runaway time has to be precisely defined.

Definition 1. 

The Thermal Runaway Time ${\mathbf{t}}_{\mathbf{TR}}$ of a Lithium ion cell is defined as the time, at which the heating rate exceeds a predefined threshold i.e., $\exists t_{\text{krit}}\in {\mathbb{R}}^{+}$ such that:

$${\displaystyle \frac{\text{d}T}{\text{d}t}\left(t\right)>\frac{\text{d}T}{\text{d}t}{|}_{\text{thres}}\forall t>t_{\text{krit}}}$$

(25)

for sufficiently high cell temperatures.

For an equidistant time grid of $\Delta t=0.1\mathrm{s}$, the heating rate $\text{d}\overline{T}/\text{d}t$ is approximated by the Equation (26).

$$\frac{\text{d}\overline{T}}{\text{d}t}\approx T^{\prime }:=\frac{\Delta \overline{T}}{\Delta t}=\frac{\overline{T}\left(t_2\right)-\overline{T}\left(t_1\right)}{t_2-t_1}.$$

(26)

The threshold is chosen here as:

$$T^{\prime }{|}_{\text{thres}}:={10}^2\mathrm{K}/\mathrm{s}.$$

(27)

The thermal runaway time $t_{TR}$ is chosen as the time at which the condition $T^{\prime }\ge T^{\prime }{|}_{\text{thres}}$ is achieved.

In

Table 3. Thermal runaway time $t_{\text{TR}}/\left[\mathrm{s}\right]$ with respect to heat transfer coefficient h and the environment al temperature $T_{\text{env}}$. (a) Profile (I), (b) Profile (II).

(a) $h/T_{\text{env}}$:	0	2	4	6	8	10	(b) $h/T_{\text{env}}$:	0	2	4	6	8	10
$273.15$	2234	2672	3431	5385	-	-	$273.15$	4173	6943	-	-	-	-
$293.15$	1899	2204	2671	3544	6376	-	$293.15$	3490	4847	10684	-	-	-
$313.15$	1585	1787	2059	2478	3272	5656	$313.15$	2821	3695	5673	-	-	-
$333.15$	1283	1398	1556	1758	2057	2539	$333.15$	2246	2708	3517	5491	-	-
$353.15$	984	1052	1133	1232	1346	1490	$353.15$	1671	1850	2217	2719	3677	6664
$373.15$	1116	1199	1282	1408	1636	1818	$373.15$	710	741	777	812	852	897

and Figure 5 the time $t_{TR}$ for the current profiles (I) and (II) is given with respect to $T_{\text{env}}$ and h. The following two cases are considered:

(1)

For $T_{\text{env}}\in [273.15,373.15]\mathrm{K}$ with $\Delta T_{\text{env}}=20\mathrm{K}$ the environmental temperature is fixed, then for $h\in \{0,1,2,3,4,5,6,8,10\}\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$ the thermal runaway temperature $t_{TR}$ is determined (Figure 5a,c).

(2)

For $h\in \{0,1,2,3,4,5,6,8,10\}\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$ the convective heat transfer coefficient is fixed, the environmental temperature $T_{\text{env}}\in [273.15,373.15]\mathrm{K}$ with $\Delta T_{\text{env}}=20\mathrm{K}$ is varied and the thermal runaway time $t_{TR}$ is computed (Figure 5b,d).

In case (1) for a fixed environmental temperature, the thermal runaway time $t_{\text{TR}}$ increases with increasing convective heat transfer coefficient h. In case (2) for a fixed convective heat transfer coefficient, the thermal runaway time $t_{\text{TR}}$ decreases with increasing environmental temperature.

3.3. Critical parameter intervals

In the last section the thermal runaway time $t_{\text{TR}}$ was computed. In this section the time evolution of the mean cell temperature $\overline{T}$ was computed, again in the cases (1) and (2) from Section 3.2. The aim is to calculate the critical mean cell temperature ${\overline{T}}_{\text{cell}}$ with respect to the parameter $T_{\text{env}}$ and h in the sense, that below the critical temperature the LIB is under safe conditions and above the LIB will show a thermal runaway. Therefore, multiple simulations of the current profiles (I) and (II) are performed. In all simulations the environmental temperature is in the interval $T_{\text{env}}\in [273.15,373.15]\mathrm{K}$ and the convective heat transfer coefficient h is in the interval $h\in [0,24]\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$. The purpose of the simulations was to compute numerically approximately intervals for the critical environmental temperature $T_{\text{env},\text{crit}}$ and the critical convective heat transfer h, analog to the cases (1) and (2) from Section 3.2 for fixed environmental temperature and variable heat transfer coefficient h (case(1)) and vice versa (case (2)). The critical parameter intervals are determined with a bisection procedure, where successive the intervals of the parameter under investigations is shrinking with respect to the fact, if a thermal runaway in the simulation occurs or not.

The results are plotted in the Figure 6 and Figure 7. In the subplots (a) and (c) for $T_{\text{env}}=273.15\mathrm{K}$, the convective heat transfer coefficient is varied in the interval $h\in [0,24]\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$ for the current profiles (I) and (II). In the subplots (b) and (d), the environmental temperature is varied over the interval $T_{\text{env}}=[273.15,373.15]\mathrm{K}$ for $h=8\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$ and $h=4\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$ for the current profiles (I) and (II). In

Table 4. Critical parameter intervals: Profile (I) (a) $T_{\text{env}}$ fixed, (b) h fixed, Profile (II) (c) $T_{\text{env}}$ fixed, (d) h fixed.

Profile (I)	($a_1$)	($a_2$)	Profile (II)	($c_1$)	($c_2$)
$T_{\text{env}}=$	$273.15\mathrm{K}$	$293.15\mathrm{K}$	$T_{\text{env}}=$	$273.15\mathrm{K}$	$293.15\mathrm{K}$
$h\in$	$[7.4725,7.48]\mathrm{W}/{\mathrm{m}}^2\mathrm{K}$	$[8.945,8.95]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$h\in$	$[3.78125,3.796875]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$[4.605,4.61]\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$
$\Delta \mathrm{h}=$	$0.0075\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$0.005\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$\Delta h=$	$0.015625\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$0.005\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$
Profile (I)	($b_1$)	($b_2$)	Profile (II)	($d_1$)	($d_2$)
$h=$	$8\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$10\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$h=$	$4\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$	$5\mathrm{W}/{\mathrm{m}}^2·\mathrm{K}$
$T_{\text{env}}\in$	$[281.05,281.15]K$	$[304,304.0625]K$	$T_{\text{env}}\in$	$[278.15,278.62]\mathrm{K}$	$[300.725,300.7625]\mathrm{K}$
$\Delta T_{\text{env}}=$	$0.1\mathrm{K}$	$0.0625\mathrm{K}$	$\Delta T_{\text{env}}=$	$0.47\mathrm{K}$	$0.0375\mathrm{K}$

some approximate intervals for the critical parameter are given.

Moreover from the temperature curves of the mean cell temperature $\overline{T}$ in the Figure 6 and Figure 7 the critical mean cell temperature is ${\overline{T}}_{\text{crit}}\approx 400\mathrm{K}$. Since the development of a complete mathematical theory is far beyond the scope of this work an extension of the theories of Semenov [52] and Frank-Kamenetskii [38] and the corresponding critical parameter will be considered in future work. In summary the following general statement can be given from the observation of the simulation results:

Statement 1. 

Let ${\overline{T}}_{\text{stat}}$ be a the asymptotic temperature profile if the limit ${\displaystyle {\overline{T}}_{\text{stat}}\left(t\right)=\underset{t\to \infty }{lim}\overline{T}\left(t\right)}$ exists. Then the following holds:

• (a) If the environmental temperature $T_{\text{env}}$ is constant and fixed and the convective heat transfer coefficient h is allowed to vary, then $\exists h_{\text{crit}}$ with:

  -

  ${\displaystyle \underset{t\to \infty }{lim}\overline{T}\left(t\right)={\overline{T}}_{\text{stat}}(t,T_{\text{env}})\forall h\ge h_{\text{crit}},}$

  -

  ${\displaystyle \underset{t\to \infty }{lim}\overline{T}\left(t\right)\to \infty \forall h<h_{\text{crit}},}$

• (b) If the convective heat transfer coefficient h is constant and fixed and the environmental temperature $T_{\text{env}}$ is allowed to vary, then $\exists T_{\text{env},\text{crit}}$ with:

  -

  ${\displaystyle \underset{t\to \infty }{lim}\overline{T}\left(t\right)={\overline{T}}_{\text{stat}}(t,h)\forall T_{\text{env}}\le T_{\text{env},\text{crit}},}$

  -

  ${\displaystyle \underset{t\to \infty }{lim}\overline{T}\left(t\right)\to \infty \forall T_{\text{env}}>T_{\text{env},\text{crit}},}$

If ${\overline{T}}_{\text{stat}}(t,h)$ or ${\overline{T}}_{\text{stat}}(t,T_{\text{env}})$ is bounded then it is (almost) constant for the current profile (I) and periodic for the current profile (II).

4. Conclusions and Outlook

In this work a model for the simulation of a thermal runaway in LIBs for a given current profile was introduced. The thermal model was extended with additional heat sources coming from several exothermic reactions inside the LIB due to the constant fuel assumption. In the electrochemical part of our model we took the length scales for the essential dynamics into account. This results in a coupled system of elliptic and parabolic partial differential equations on the particle, the electrode and the cell domain respectively. The coupling of the equations of the particle domain with the the equations of the electrode domains, as well as the the equations of the electrode domain with the equations of the cell domain is solved with suitable homogenization procedures using the spatial averaging theorem which is based on the assumption of spatial homogeneity. Despite the fact that the advantages of this model in the spatial resolution are not used in this work, it is the intention of the authors to consider spatial localized effects in the dynamics of a thermal runaway in future work with this model. In this work we show the difference between the old model Model A and the model with additional contributions for exothermic reactions Model B. Furthermore, we have shown that Model B is able to show a thermal runaway in contrary to Model A. We defined and computed the thermal runaway time with respect to some system parameters as well as the critical cell temperature which can be interpreted as a marker for the cell stability with respect to thermal dynamics. Furthermore, a coarse classification of thermal runaway is given with the Model B. In summary, in these simulations, which were performed using COMSOL Multiphysics, we have shown that these additional exothermic terms in the numerical FEM model can resolve such phenomena for a specific current load profile. Due to the constant fuel assumption, this model is not able to resolve the decay period of a thermal runaway after all fuel is consumed. Also, the change in the amount of fuel is not considered in this model. For future work, this first model describing the thermal runaway of LIB’s will be improved by skipping the constant fuel assumption. A detailed consideration of the underlying battery chemistry and corresponding exothermic reaction kinetics is planned in addition to determining a criteria to predict a thermal runaway with the information given by the simulations. Furthermore, it is planned to compare the simulations of this model with experimental results in a future work.