A Multiphysics Model Simulating the Electrochemical, Thermal, and Thermal Runaway Behaviors of Lithium Polymer Battery

: Lithium-ion batteries (LIBs) have circumvented the energy storage landscape for decades. However, safety concerns about liquid–electrolyte-based LIBs have challenged their mobilization. Lithium polymer (LiPo) batteries have gained rising interest due to their high thermal stability. Despite an array of commercially available LiPo batteries, limited studies have ventured into modeling. Numerical simulations allow low-cost optimization of existing battery designs through parameter analysis and material conﬁguration, leading to safer and more energy-efﬁcient batteries. This work examined the electrochemical, thermal, and thermal runaway behavior of a lithium cobalt oxide cathode, graphite anode, and poly(vinylidene ﬂuoride-hexaﬂuoropropylene) electrolyte pouch-type LiPo battery using COMSOL Multiphysics ® , and validated results with experimental data. The simulated potential curve exhibited strong agreement with experiment data, while the temperature proﬁle during discharge displayed qualitative discrepancies rationalized by the reversible heat generation. Thermal runaway simulations via oven tests revealed that the highest heat generation is from the cathode–electrolyte reaction, while the solid electrolyte interface decomposition initiates the heat generation process. These results suggest a thorough selection of cathode and electrolyte material to heighten battery safety. Overall, the developed models can be utilized as design tools to investigate various chemistries and designs to estimate the behavior and performance of batteries.


Introduction
Energy demand is increasing as new technologies are being developed rapidly. Researchers make vigorous efforts to tap into energy storage to keep up with the rise in energy demand. Nearly all portable devices operate with batteries, and though no holy grail battery has been developed, there are boundless areas for improvement in their design and structure. Among these batteries, lithium-ion batteries (LIBs) have been the focus of studies because of their high energy capacity [1]. LIBs deliver energy by intercalating and de-intercalating lithium ions in a rocking chair mechanism ( Figure S1). Lithium polymer (LiPo) batteries that use gel polymer electrolytes (GPEs) are of increasing interest because of their non-leakage, flexibility, stability, and lower resistance between electrode and electrolyte compared to solid electrolytes, and because of their high ionic conductivity [2]. As a result of heightened interest, numerous studies on the material modification and venturing of new materials have paved the way for different LiPo battery designs and frameworks. Despite an increase in diverse LiPo batteries, there is a limited volume of studies about simultaneously modeling these batteries' behaviors. fire at worst. To mitigate thermal runaway, research pursuits aim to increase the onset temperature and the tipping point temperature while reducing the temperature peak.
Battery management systems must suppress thermal runaways to avoid fire and/or battery explosion. Furthermore, as the demand for electrochemical battery storage transition from mobile use to transportation and even health care, user safety is an important factor in battery performance. Thus, there have been rigorous efforts to assist with the experimental optimization of battery design for safety purposes through the analysis of thermal mechanisms by using mathematical models [7,8,[21][22][23]. One study by Hatchard et al. simulated the thermal runaway of a lithium cobalt oxide (LCO) battery using a Fortran numerical model [24]. Although their model validated cylindrical oven test results, prismatic-type batteries lacked experimental comparison. Another study capturing LCO batteries under thermal runaway was conducted by Zhang et al. using a multiphysics model [20]. Although the temperature profile shows close agreement, the validation was referenced to previous literature findings. Other studies consider using a thermal resistance model for thermal runaway such as the works of Chen et al. [24]. Thermal runaway has been successfully modelled. However, thermal runaway is factored by multiple phenomena such as electrochemical and thermal behaviors.
To address the challenges in the venture for novel materials for reliable and safe energy storage technologies, endeavors for modeling and optimization should be given equal and utmost importance. Thus far, little effort has been made to design a robust yet straightforward mathematical model to understand the behavior of commercially available battery models. Additionally, studies that have modeled commercially available batteries often validate experimental results with past literature, lacking validation from actual battery tests.
This study reports the first multiphysics model for a commercially available mobile LiPo battery with an LCO, LiCoO 2 , cathode, a graphite anode, and a poly(vinylidene fluoride-hexafluoropropylene) (PVdF-HFP) electrolyte, investigating the electrochemical, thermal, and thermal runaway (through external heating) behaviors and validating them with experimental results. The mathematical model can be utilized further to understand the battery's behavior under different conditions. The model also intends to create opportunities for improving battery design and optimizing battery material composition for enhanced performance and heightened safety.

Experimental Methods
A commercial LiPo battery, GEB 585460, produced by General Electronics Battery Co., Ltd. (Shenzhen, China), was used as the framework battery for the study. The schematics of the experimental research are illustrated in Figure S3. The battery specifications are summarized in Table 1. The electrochemical analysis was performed using an AUTOLAB PGSTAT302N potentiostat. Working voltage at a constant ambient temperature was determined using the potentiostat under galvanostatic discharge operation, where the working electrode was placed on the positive terminal of the cell and counter and reference electrodes on the negative terminal. Thermal analysis by analyzing the battery's time-dependent temperature during discharge was performed using two K-type stainless steel thermocouples (WRNT-10) situated at the center of the battery cell on two sides, as illustrated in Figure  S4. The temperature was logged using a Fluke 2625A HYDRA data logger. The average temperature from the thermocouples designates the average battery surface temperature. The electrochemical and thermal analysis of the battery was conducted simultaneously at a discharge rate of 1.0 C. Thermal runaway experimental analysis was carried out on an oven test coupled with an identical 2-thermocouple arrangement in the thermal analysis to measure the battery surface temperature. The oven tests for the battery proceeded at two temperature points-at (i) 403.15 K, where the PVdF separator breaks down [25], and at (ii) 423.15 K, where the onset temperature of thermal runaway of LCO-graphite takes place [19].

Numerical Methods
The electrochemical model provides the necessary insights to understand the relationship between the current distribution along with the battery phases, the different current densities, and the battery's overpotential. The fundamental equations for the electrochemical model of the LiPo battery are summarized in Table 2. Table 2. LiPo battery modeling equations for electrochemical behavior.
Ionic conductivity can be determined as a function of lithium-ion concentration in the electrolyte phase using Equation (12) [12].
The lithium diffusion coefficient in the electrolyte phase as a function of temperature and lithium concentration in the electrolyte phase is expressed in Equation (13) [4].

of 24
The equilibrium potential of the electrodes, expressed as the open circuit potential (OCP), is expressed as a function of the surface state of charge. The equations were obtained from literature based on experimental fitting with units of volts [26]. The OCP of the electrodes can be determined from the stoichiometric coefficient, which designates the ratio of the solid-phase lithium-ion concentration on the local electrode surface to the maximum available lithium-ion concentration the electrode can hold [27]. The OCP, as a function of lithiation for the cathode and anode, is illustrated in Figure 1, established from their respective theoretical equations [26]. κ = 4.1253x10 + 5.007x10 c Li − 4.7212 10 (c Li ) + 1.5094 10 -3 (c Li E ) 3 − 1.6018 10 -4 (c Li E ) 4 .
The lithium diffusion coefficient in the electrolyte phase as a function of temperature and lithium concentration in the electrolyte phase is expressed in Equation (13) [4].
The equilibrium potential of the electrodes, expressed as the open circuit potential (OCP), is expressed as a function of the surface state of charge. The equations were obtained from literature based on experimental fitting with units of volts [26]. The OCP of the electrodes can be determined from the stoichiometric coefficient, which designates the ratio of the solid-phase lithium-ion concentration on the local electrode surface to the maximum available lithium-ion concentration the electrode can hold [27]. The OCP, as a function of lithiation for the cathode and anode, is illustrated in Figure 1, established from their respective theoretical equations [26]. For the positive electrode, the equilibrium potential is expressed in Equation (14), where the surface state of charge of the LCO cathode, , is the stoichiometric coefficient, x, of LixCoO2 [26].
The cell's total current density can be computed by taking the sum of the current density from both the solid and electrolyte phase.
The thermal model allows us to understand the battery's thermal behavior under various conditions. Heat generation develops from the chemical reactions inside the battery cell. The thermal model enables the prediction of the battery's thermal distribution. The equations for satisfying the energy balance for LiPo battery thermal modeling are summarized in Table 3. For the positive electrode, the equilibrium potential is expressed in Equation (14), where the surface state of charge of the LCO cathode, x p , is the stoichiometric coefficient, x, of Li x CoO 2 [26].
For the negative electrode, the equilibrium potential is expressed in Equation (15), where the surface state of charge of the graphite anode, x n , is the stoichiometric coefficient, y, of Li y C 6 [26]. U eq,n = 0.13966+0.68920 exp (−49.20361x n )+0.41903 exp (−254.40067x n )− exp (49.97886x n − 43.37888) − 0.028221arctan(22.52300x n − 3.65328) −0.01308arctan(28.34801x n − 13.43960) (15) The cell's total current density can be computed by taking the sum of the current density from both the solid and electrolyte phase.
The thermal model allows us to understand the battery's thermal behavior under various conditions. Heat generation develops from the chemical reactions inside the battery cell. The thermal model enables the prediction of the battery's thermal distribution. The equations for satisfying the energy balance for LiPo battery thermal modeling are summarized in Table 3. Table 3. LiPo battery thermal modeling equations.

Name Equation Reference
General energy balance ρc p ∂T ∂t = λ∇ 2 T + q (16) [22,28] Battery heat generation q = q B (17) [11] Heat generation in battery model q B = q rev +q irrev (18) [20] Reversible heat q rev = ai n T ∂U eq ∂T (19) [11] Irreversible heat q irrev = q pol +q ohm (20) [20] Polarization heat q pol = ai n η (21) [11] Ohmic heat [28] As the battery was allowed to dissipate heat to the environment, the heat transfer in the battery predominantly accounts for the convective heat transfer to the surrounding air and the radiative heat transfer through the surface areas. However, heat transfer through radiation is negligible and is often assumed to be zero for ease of computation. Thus, according to Newton's cooling law, the boundary conditions for the thermal model are described as follows (Equations (24) and (25)): The entropic coefficients to account for the reversible heat (Equation (19)) use literaturederived values grounded on experimental work. The entropic coefficient equations are a function of the degree of lithiation or the surface state of charge and are expressed in volts per Kelvin [29][30][31]. The entropic coefficient as a function of lithiation is illustrated in Figure 2. For the positive electrode, the entropic coefficient is expressed in Equation (26), where the degree of lithiation of the LCO cathode, x p , is the stoichiometric coefficient, x, of Li x CoO 2 [26]. The equation was derived from fitting experimental works of the variation of half-cell potential vs. Li/Li+ to the curves of Guo et al. [26]. The heat coming from the battery's electrochemical model is vital to the heat generation for the thermal model. In contrast, the thermal model's temperature distribution affects the electrochemical model's current density. This coupling of both models allows a synergistic approach to investigate the multiple physical phenomena of For the negative electrode, the entropic coefficient is expressed in Equation (27), where the degree of lithiation of the graphite anode, x n , is the stoichiometric coefficient, y, of Li y C 6 [32]. The equation was based on the experimental data reported by Reynier et al. [33] by determining the variation of the half-cell potential vs. Li/Li+ and fitted to the works of Rheinfeld et al. [32].
The thermal model uses a 3-dimensional approach to heighten model reliability. The overall heat energy balance expressed in Equation (16) was used as the foundation of the thermal analysis up to the thermal runaway modeling. The thermal model neglects the negligible heat dissipation from radiative heat transfer [23]. The irreversible heat generation is limited to active polarization heat and ohmic heat. In this study, the heat of mixing and enthalpy heating were not considered because of their insignificant amount. The effect of the heat of mixing is negligible for LIBs because of the constrained concentration gradients present [34].
The heat coming from the battery's electrochemical model is vital to the heat generation for the thermal model. In contrast, the thermal model's temperature distribution affects the electrochemical model's current density. This coupling of both models allows a synergistic approach to investigate the multiple physical phenomena of the battery's dynamics. The physical and thermodynamic parameters used in the coupled electrochemical-thermal models are summarized in Table 4. Transference number of lithium-ion (t Li ) 0.435 [28] Cathodic charge transfer coefficient (α C ) 0.5 0.5 [38] Anodic charge transfer coefficient (α A ) 0.5 0.5 [38] Electronic conductivity in the solid phase (σ S ) 1 3.8 × 10 7 2 10 2 100 1 6 × 10 7 1 - [35] 2 - [38] Ionic conductivity (κ) (Equation (12)) [8] Energies 2023, 16  Li x CoO 2 Equilibrium potential (U eq,p ) (Equation (14)) [26] Li y C 6 Equilibrium potential (U eq,n ) (Equation (15) Thermal runaway occurs when there is a spontaneous increase in the cell's temperature due to heat generation from multiple chemical reactions. Thermal runaway is one of the battery failure modes that can inflict serious harm to the user. Recognizing and predicting when a battery cell will undergo thermal runaway is vital. Table 4 summarizes the framework equations for thermal runaway modeling for a LiPo battery. Table 5 summarizes the framework equations for thermal runaway modeling for a LiPo battery. The overall heat generation from the thermal runaway chemical reactions proposed by Kim et al.'s mathematical model was built from the Arrhenius equation's principles [40]. The total heat from Equation (17) is then substituted to q TR , designated as the total thermal runaway heat generation. The total heat generation is the sum of the heat (Equation (28)) from (i) SEI decomposition reaction (Equation (29)), (ii) reaction of the negative electrode and the electrolyte (Equation (30)), (iii) reaction of the positive electrode and the electrolyte (Equation (31)), (iv) electrolyte decomposition (Equation (32)), and (v) the reaction between the binder and the negative electrode [41]. The reaction between the binder and negative electrode rarely occurs because lithium is consumed more quickly and more completely during the reaction between the electrolyte and the negative electrode [34]. Therefore, the heat from the exothermal reaction between the binder and the negative electrode was neglected in modeling. Table 5. LiPo battery thermal runaway modeling equations.

Name Equation Reference
Heat generation in thermal runaway q TR = q SEI +q a−e +q c−e +q e The mathematical simulations were achieved using COMSOL Multiphysics 5.4 software coupling the batteries and fuel cells, and the heat transfer module. The models used the same geometry and material composition as the commercial LiPo battery for experimental analysis. The P2D electrochemical and 3D thermal model, including the thermal Energies 2023, 16, 2642 9 of 24 runaway model, was coupled with temperature and heat generation variables. The thermal runaway used the ODE domain to model and solve the exothermic chemical reactions. The mathematical model uses the backward differentiation formula as the time-stepping method with a maximum order of five and a minimum order of two, using a relative tolerance of 10 −6 and an absolute error tolerance of 10 −8 , as suggested by the work of Zhang et al. [20]. Working voltage and discharge time were analyzed for the electrochemical model, while temperature as a function of time was analyzed for thermal and thermal runaway models. The electrochemical and thermal model employed a discharge C-rate of 1.0 C. For the thermal runaway model, the heat generation rate as a function of time was also investigated.
The simulation was run using a computer workstation with an Intel(R) Core i5-10300H CPU at 2.50 GHz, 2496 Mhz, 4 Cores, and 8 Logical Processors. The schematic illustration of the models' geometries is presented in Figure S5. The convergence of the computational solution of the partial differential equations is highly dependent on the number and geometry of the grids and the mesh analysis due to the battery's high aspect ratio via the finite element method. To provide reliability and evaluate the computational speed and quality of the model, a preliminary grid sensitivity study was conducted at varying physics-controlled mesh at COMSOL Multiphysics. The RMSE, comparing experimental and simulated discharge potentials at a different mesh, was calculated using a 1 C discharge rate and 100 s time step computation, tabulated in Table 6. The grid independence study demonstrates that the RMSE was roughly 0.06 V, and the lowest RMSE was achieved using a normal mesh type. As such, a normal physics-controlled mesh was utilized in the study as fine mesh would require four times the computational time with roughly the same results as a coarse-type running half the duration of the normal-type while finer mesh results in computer memory overcapacity.

Electrochemical Behavior
The discharge curve for both the experiment and the simulation can be seen in Figure 3. The calculated RMSE for the 1.0 C rate is 0.06 while the 0.5 C rate is 0.07, demonstrating that the generated model can be used as a reliable first approximation for battery behavior prediction. Figure 3 shows that the experimental discharge curve for both the 0.5 C and 1.0 C rate tests drops to the cut-off voltage earlier than the simulation. This may be factored by an overestimate of the specific surface of the cell affected by the imposed current density [43]. Additionally, it could be attributed to the aging of the cell since it was cycled for a certain period prior to the experiments [44]. However, in general, the current model is able to represent the electrochemical behavior of the cell during discharge. Nonetheless, incorporating descriptions of bulk lithium diffusion and surface kinetics may enhance the accuracy of the simulation.  To investigate the battery behavior at higher-discharge rates, a discha a 10 C rate was simulated, which can be seen in Figure S6a. Compare simulated profile of the 0.5 and 1.0 C-rates in Figure 3, it is evident that th activation losses as the C-rate increases. Figure S6b shows the distribution lithium-ion concentration inserted at the electrodes during the simulated 1 rate at various times. A slight decrease in curvature can be seen on the nega near the separator, especially during the early stages of discharge. As the p further from the separator, the average lithium-ion concentration beco diffuses evenly. During the end of discharge, all lithium ions inserted from electrode transfer to the positive electrode. The distribution of the averag concentration in the electrolyte can be seen in Figure S6c during the 10 discharge. At the onset of the discharge process, the lithium-ion concen electrolyte increases drastically in the negative electrode. The opposite can positive electrode, where a sharp decrease is evident. The incremental chan ion concentration decreases slightly as time increases, and the difference profound as the end of discharge draws near. Furthermore, it is eviden distance draws closer to the separator, the lithium-ion concentration of electrode decreases, which is opposite to the case of the positive electrode.

Thermal Behavior
The electrochemical and thermal behaviors were simultaneou experimentally and computationally. The calculated RMSE was 1.5153 K us of 100 s for the 1.0 C discharge. Figure 4 illustrates the average battery surfac during 1.0 C discharge. The experimental and simulated battery surface Experimental discharge curves (broken lines) and simulated curves (solid lines). 1 C discharge C-rate (black lines) and 0.5 C discharge C-rate (red lines).
In both C-rates, a slight decrease in simulated working potential can be seen at t = 1000 s for 1.0 C and t = 2000 s for 0.5 C, this abrupt change in slope denotes that a two-step discharge profile is evident for the battery chemistry [45][46][47]. This phenomenon was not profound in the experimental discharge except for a very minimal curve observed at 0.5 C at t = 4000 s. The voltage plateau defines the operating voltage the battery is suited for. The voltage plateaus are observed in both C-rates, but the plateau at the 0.5 C rate is more defined. A decrease in the voltage plateau is evident as an increase in the discharge rate was imposed. This is primarily due to higher polarization at the higher discharge C-rate, which leads to more significant polarization heat generation, further propagating to higher battery temperature operation.
To investigate the battery behavior at higher-discharge rates, a discharge profile for a 10 C rate was simulated, which can be seen in Figure S6a. Compared to both the simulated profile of the 0.5 and 1.0 C-rates in Figure 3, it is evident that there are higher activation losses as the C-rate increases. Figure S6b shows the distribution of the average lithiumion concentration inserted at the electrodes during the simulated 10 C discharge rate at various times. A slight decrease in curvature can be seen on the negative electrode near the separator, especially during the early stages of discharge. As the position moves further from the separator, the average lithium-ion concentration becomes flat as it diffuses evenly. During the end of discharge, all lithium ions inserted from the negative electrode transfer to the positive electrode. The distribution of the average lithium-ion concentration in the electrolyte can be seen in Figure S6c during the 10 C simulated discharge. At the onset of the discharge process, the lithium-ion concentration in the electrolyte increases drastically in the negative electrode. The opposite can be seen in the positive electrode, where a sharp decrease is evident. The incremental change in lithium-ion concentration decreases slightly as time increases, and the difference becomes less profound as the end of discharge draws near. Furthermore, it is evident that, as the distance draws closer to the separator, the lithium-ion concentration of the negative electrode decreases, which is opposite to the case of the positive electrode.

Thermal Behavior
The electrochemical and thermal behaviors were simultaneously analyzed experimentally and computationally. The calculated RMSE was 1.5153 K using a timestep of 100 s for the 1.0 C discharge. Figure 4 illustrates the average battery surface temperature during 1.0 C discharge. The experimental and simulated battery surface temperature starts at about 303.5 K and reaches a maximum temperature of 311 K at the end of discharge, corresponding to nearly a 7.50 K increase in cell temperature. It can be seen in Figure 4 that the simulated results do not agree well with the experimental results, thus having a maximum temperature difference of around 3.0 K at about t = 600 s. It should not be demerited that the RMSE of 1.5 K can be at par with modern simulations with RMSE values ranging from 1.0 K [48] to 1.75 K [49]. The sources of discrepancy can be traced to various factors. One of these is that the electrochemical and thermal properties used in the literature from several years ago may no longer hold accurate values for modern batteries. Another factor is that most models of batteries have been focused on cylindrical batteries, such as 18650 and coin-type batteries, and less on pouch-type batteries [50].
Energies 2023, 16, x FOR PEER REVIEW 11 of 24 used in the literature from several years ago may no longer hold accurate values for modern batteries. Another factor is that most models of batteries have been focused on cylindrical batteries, such as 18650 and coin-type batteries, and less on pouch-type batteries [50]. Other sources of error may be due to the limitations of the simulated model itself. The generated model presumed uniform heat generation within the active battery material, which may cause inaccuracies. Assuming homogeneity of the system in calculating the thermodynamic properties may also factor into the errors. The thermodynamic property of the active battery material used is a volumetric average of the component materials' properties. This assumption considers the active cell material to be a unified body that speeds up the calculation by reducing the number of grids [51]. A study by Capron examined the effects of homogenous cells and discrete cells, where the cell is split into smaller regions [48]. It was found that both the discharge curve and surface temperature profile of the cell have good agreement, with experimental results validating the reliability of a homogeneous cell with a simpler simulation design. However, the temperature profile of the internal region of the cell was captured fully using the discrete model. The merit of multiple layering is that it allows the determination of the point-bypoint temperature gradient of the battery during operation, and portrays the temperature flow pattern at a certain time [13]. This analysis would require high-precision thermal cameras to determine experimental battery temperature and compare it with 3D simulation [52].
The accuracy of selecting the thermal boundary condition may also contribute to the discrepancies. Furthermore, the entropic change coefficient is a crucial parameter in predicting the heat generation inside the battery. The literature-provided entropic change coefficients were usually additionally fitted to experimental curves to heighten accuracy [26]. Unlike the graphite anode with various available entropic change coefficient equations, such as in the works of Rheinfeld et al. and Guo et al., the LCO has been given Other sources of error may be due to the limitations of the simulated model itself. The generated model presumed uniform heat generation within the active battery material, which may cause inaccuracies. Assuming homogeneity of the system in calculating the thermodynamic properties may also factor into the errors. The thermodynamic property of the active battery material used is a volumetric average of the component materials' properties. This assumption considers the active cell material to be a unified body that speeds up the calculation by reducing the number of grids [51]. A study by Capron examined the effects of homogenous cells and discrete cells, where the cell is split into smaller regions [48]. It was found that both the discharge curve and surface temperature profile of the cell have good agreement, with experimental results validating the reliability of a homogeneous cell with a simpler simulation design. However, the temperature profile of the internal region of the cell was captured fully using the discrete model. The merit of multiple layering is that it allows the determination of the point-by-point temperature gradient of the battery during operation, and portrays the temperature flow pattern at a certain time [13]. This analysis would require high-precision thermal cameras to determine experimental battery temperature and compare it with 3D simulation [52].
The accuracy of selecting the thermal boundary condition may also contribute to the discrepancies. Furthermore, the entropic change coefficient is a crucial parameter in predicting the heat generation inside the battery. The literature-provided entropic change coefficients were usually additionally fitted to experimental curves to heighten accuracy [26]. Unlike the graphite anode with various available entropic change coefficient equations, such as in the works of Rheinfeld et al. and Guo et al., the LCO has been given less attention. Currently, the fundamental entropic change coefficient used for the LCO cathode is a modification of what was used back in the 1990s. [53] This study, however, was limited to using the commonly used theoretical equations for its parameters and comparing the results to experimental data. Recent experimental studies of the LCO cathode have shown endothermic heat release from SOC > 0.8 and exothermic heat release from SOC < 0.8 [54]. This scenario cannot be reflected in Figure 2a, which shows that modern commercial batteries require rigorous experimental determination of critical electrochemical and thermal parameters to secure the model's reliability and accuracy.
The electrochemical and thermal models were coupled with temperature and heat generation. The simulation results of the discharge and temperature profile of the 1.0 C rate can be seen in an overlapping figure in Figure S7. A slight increase in working potential at around t = 800 s may be due to the sudden temperature decrease of the battery owing to the endothermic reaction.
The general heat balance equation of Equation (16) determines the battery temperature. It is factored heavily on the heat generation coming from the battery consisting of the reversible heat generation and irreversible heat generation denoted in Equations (19) and (20), respectively. The effects of these two primary heat sources on the average battery surface temperature can be seen in Figure 5, where the values were determined using numerical simulation. The simulations exhibit an average heat generation of 0.90 W m −3 during 1 C discharge. As C rate increases, current density and ohmic heating also increase, thus resulting in greater energy lost as heat. This effect was validated as the average total heat generation at 10 C discharge rose to 13.20 W m −3 , roughly fifteen times that of 1 C. Additionally, 0.5 C discharge presented an average total heat generation of 0.39 W m −3, which is around half that of 1 C discharge. This value is expected since the lower current density imposes lower ohmic heating. Determining the optimum and maximum operating C-rate of the battery to prevent overheating is, thus, of high importance. The electrochemical and thermal models were coupled with temperature and heat generation. The simulation results of the discharge and temperature profile of the 1.0 C rate can be seen in an overlapping figure in Figure S7. A slight increase in working potential at around t = 800 s may be due to the sudden temperature decrease of the battery owing to the endothermic reaction.
The general heat balance equation of Equation (16) determines the battery temperature. It is factored heavily on the heat generation coming from the battery consisting of the reversible heat generation and irreversible heat generation denoted in Equations (19) and (20), respectively. The effects of these two primary heat sources on the average battery surface temperature can be seen in Figure 5, where the values were determined using numerical simulation. The simulations exhibit an average heat generation of 0.90 W m −3 during 1 C discharge. As C rate increases, current density and ohmic heating also increase, thus resulting in greater energy lost as heat. This effect was validated as the average total heat generation at 10 C discharge rose to 13.20 W m −3 , roughly fifteen times that of 1 C. Additionally, 0.5 C discharge presented an average total heat generation of 0.39 W m −3, which is around half that of 1 C discharge. This value is expected since the lower current density imposes lower ohmic heating. Determining the optimum and maximum operating C-rate of the battery to prevent overheating is, thus, of high importance. During discharge, the surface SOC of the positive electrode increases while the surface SOC of the negative electrode decreases. During the initial discharge, both electrodes exhibit a negative entropic change coefficient in the regions where the SOC of the positive electrode ranges from 0.4 (initial) to 0.5 and of the negative electrode from 0.76 (initial) to 0.1. This negative entropic change coefficient exhibits exothermic heat release; the reversible heat generation and the irreversible heat generation both rise as t increases to 550 s and 600 s, respectively. The combined effect of the two heat sources explains the rapid temperature rise from the initial to t = 750 s. This initial temperature rise effect is consistent with those reported in the literature [55]. heat generation and the irreversible heat generation both rise as t increases to 550 s and 600 s, respectively. The combined effect of the two heat sources explains the rapid temperature rise from the initial to t = 750 s. This initial temperature rise effect is consistent with those reported in the literature [55].
After the initial temperature rise associated with the increase in both reversible and irreversible heat, the simulated temperature drops from t = 750 s to t = 1200 s. This is mainly due to the reversible heat drop to a negative value designating an endothermic reaction. Unlike the irreversible heat that is always associated with heat production, reversible heat may either absorb or release heat. The entropic change coefficient of the positive electrode exhibits a positive value around SOC = 0.5 to 0.55, indicating an endothermic reaction. The battery heat absorption can be seen when the reversible heat suddenly drops after t = 550 s, and a negative reversible heat generation value is present at t = 800 s to t = 1600 s. Furthermore, as the battery temperature increases, the rising temperature gradient between the cell and surrounding accelerates heat dissipation of the battery, which in turn prevents temperature rise. In line with experimental results, from t = 1500 s to t = 2500 s, a slow temperature increase was present. Additionally, a stagnant temperature of about 308.5 K from t = 1800 s to t = 2400 s was observed. This phenomenon may be explained by a combination of possible endothermic heat by the electrodes and increased heat dissipation through convection. The former can be determined experimentally through potentiometric and calorimetric techniques [56].
During the final stage of discharge, the local surface SOC of the positive electrode reaches 1.0 while that of the negative electrode drops to zero. Despite the entropic heat coefficient of the negative electrode having a positive value signifying endothermal reaction, the negative exothermal value of the positive electrode holds a more significant amount than the reversible heat generation. The reversible heat generation rises from t = 1000 s until it becomes positive again at t = 1600 s. In return, the temperature rises back up from t = 1200 s, due to an escalation of both irreversible heat generation and reversible heat generation. The irreversible heat quickly increases due to a sharp surge in the magnitude of the overpotential as the discharge comes to an end. Consequently, the battery temperature also rises.
Overall, the thermal model provides viable results for a quick-glance view of the battery's temperature profile during operation. The simulation can be improved by (i) incorporating a non-homogeneous active battery by developing a layer-by-layer cell to impose point-by-point temperature mapping, (ii) determining experimentally the physical properties of cell materials (i.e., thermal conductivity) to improve the accuracy and specificity of the model, (iii) including non-uniform heat generation to account for temperature gradients within the active material and, more importantly, (iv) fine-tuning the entropic change coefficient of the electrodes by joint experimental and numerical optimization to alter the overall curve of the temperature profile to tolerably match the experimental results. Figure 6 shows the experimental and simulated evolution of the average battery surface temperature during the oven tests conducted at 403.15 and 423.15 K. The initial battery temperature was approximately 301.15 K for both tests. During the 403.15 K oven test, a considerable temperature spike was not observed in the experimental and simulated data. The 403.15 K oven test has an RMSE value of 3.04 K using a timestep of 100 s. The principal error was in the t = 500 to 700 s region, where a maximum temperature difference of 9 K was observed. As seen in Figure 6, both experimental results developed a slight temperature increase at around t = 700 s. As both of the experimental results have a higher average surface temperature above their respective oven temperature, thermal runaway was triggered by the thermal abuse of external heating of the oven. This may indicate a sudden heat generation due to SEI decomposition and/or anode-electrolyte reaction, indicating the presence of the tipping point temperature. The simulated battery temperature during the 423.15 K oven test sho accelerated temperature peak starting at around t = 650 s and ultimately peaking a t = 700 s. The temperature at the peak was about 773.15 K, which is 350 K higher t oven temperature, indicating the occurrence of thermal runaway. However, the r spike effect was not reflected in the 423.15 K experimental battery oven test. Th occurrence may be due to (1) the theoretical equations primarily based on cylindri cells, (2) homogenous average battery thermodynamic properties, (3) the condition imposed in the model, and (4) thermal expansion and mechanical stre considered in the model. The theoretical equations were adapted from the w Hatchard et al., where the simulated temperature profile was compare experimental data using 18650 cylindrical cells [57]. In their study, a temperatu was not observed during an oven test performed at 423.15 K. However, their showed otherwise, which is reflected in this study. Despite the difference, the stressed that the equations might be applied to pouch-type cells without compa experimental results from such geometry. Furthermore, the differences in cell ra thicknesses could affect the occurrence of thermal runaway owing to the differenc surface area-to-volume ratio, which consequently affects the dissipation of the heat generation. The thermodynamic properties of the cell play a crucial role accuracy of both the thermal and thermal runaway models. The model used in th assumes homogeneous average thermodynamic properties in the active battery m in a block-type geometry rather than the actual layer-by-layer of a pouch-type This assumption results in faster numerical computation by sacrificing the accu internal heat transfer [58]. To improve model reliability, a layer-by-layer geometr pouch cell can be carried out to account for the series of active cell materials as con by Macdonald et al. for an LFP battery [59]. The results of their work were agreement with the experimental findings of Joachin et al. [60]. These studies h the importance of incorporating geometry-specific conditions to improve however, they also challenge the model to a longer computational time and The simulated battery temperature during the 423.15 K oven test showed an accelerated temperature peak starting at around t = 650 s and ultimately peaking at about t = 700 s. The temperature at the peak was about 773.15 K, which is 350 K higher than the oven temperature, indicating the occurrence of thermal runaway. However, the runaway spike effect was not reflected in the 423.15 K experimental battery oven test. This nonoccurrence may be due to (1) the theoretical equations primarily based on cylindrical type cells, (2) homogenous average battery thermodynamic properties, (3) the no-vent condition imposed in the model, and (4) thermal expansion and mechanical stresses not considered in the model. The theoretical equations were adapted from the works of Hatchard et al., where the simulated temperature profile was compared with experimental data using 18650 cylindrical cells [57]. In their study, a temperature spike was not observed during an oven test performed at 423.15 K. However, their model showed otherwise, which is reflected in this study. Despite the difference, the authors stressed that the equations might be applied to pouch-type cells without comparison to experimental results from such geometry. Furthermore, the differences in cell radii and thicknesses could affect the occurrence of thermal runaway owing to the differences in the surface area-to-volume ratio, which consequently affects the dissipation of the internal heat generation. The thermodynamic properties of the cell play a crucial role in the accuracy of both the thermal and thermal runaway models. The model used in this study assumes homogeneous average thermodynamic properties in the active battery material in a block-type geometry rather than the actual layer-by-layer of a pouch-type battery. This assumption results in faster numerical computation by sacrificing the accuracy of internal heat transfer [58]. To improve model reliability, a layer-by-layer geometry of the pouch cell can be carried out to account for the series of active cell materials as conducted by Macdonald et al. for an LFP battery [59]. The results of their work were in close agreement with the experimental findings of Joachin et al. [60]. These studies highlight the importance of incorporating geometry-specific conditions to improve results; however, they also challenge the model to a longer computational time and greater complexity of the simulation equations.

Thermal Runaway Behavior
The venting conditions of a 18650 cell were studied by Coman et al. The simulations of thermal runaway with and without venting electrolyte were compared to determine the effect of gas pressure. The venting parameter, which is another heat dissipation factor in the cell, changes when a temperature spike occurs [41]. The simulation results under gasventing conditions corresponded well with the experimental time of the temperature spike. However, the simulated value of the temperature spike overshoots the experimental result by around two-fold, regardless of whether venting or no-vent conditions were imposed. Furthermore, imposing the vent condition requires experimental determination of the ejecta properties specific to the battery electrolyte. Overall, the venting condition is a powerful addition to the model to improve the agreement of the time of occurrence of thermal runaways.
The volumetric changes in battery configuration during thermal runaway can also deviate from its behavior. The expansion of the battery due to thermal stresses changes the heat transfer mechanism, as illustrated by the study of Xu et al., in which thermal gradients were affected by thermal expansion [61]. However, as in many studies, Xu et al. did not consider the temperature variations on volume changes of the battery geometry during operation. Additionally, simulations with and without thermal expansion consideration were not compared. The mechanical abuse of the cell components alters the temperature profile during thermal runaway [62]. In another study, Liu et al. developed a thermal runaway model due to short-circuit from mechanical abuse overloading [7]. An auxiliary mechanical model was added to the conventional thermal runaway model and resulted in a temperature difference graph similar to the experimental results. The added mechanical stress in the model indirectly affects the thermal runaway model through the short-circuit model, adding another heat generation variable.
The high temperature difference between the simulated and experimental results during thermal runaway makes it impractical to compute the RMSE for the 423.15 K oven test. The difference in temperature profiles of the two simulated oven tests demonstrates that different exothermic reactions occurred at the battery, where the severity of one reaction may have a more significant impact depending on the temperature the battery was exposed to. The 403.15 K oven test demonstrated sluggish exothermic reactions as the heat generated by the runaway reactions is quickly released in time without a temperature spike reflected in Figure 6. However, the simulated 423.15 K oven test illustrates more intensive exothermic reactions leading to intensive heat accumulation, ultimately raising the cell temperature at an accelerated pace, as seen in Figure 6. A supplementary oven test at 443.15 K ( Figure  S8) illustrates the non-occurrence of a temperature spike for the experiment. Furthermore, results showed a higher simulated temperature peak compared to the 423.15 K oven test due to the higher oven temperature. Moreover, the occurrence of a temperature spike at simulation occurred earlier than during the 423.15 K oven test.
To verify that the heat generation due to thermal runaway was within the active cell component, a 3D temperature distribution of the battery at pre-and post-thermal runaway was developed. Figure S9 shows the temperature distribution of the battery surface during the 423.15 K oven test. The temperature profile before thermal runaway ( Figure S9a) shows that a higher temperature was observed at the current collectors than the battery material; this may be due to the difference in thermal conductivity and geometry of the current collectors and the metal cover of the battery. The temperature gradient demonstrates possible heat transfer from the current collector to active battery material. As the battery proceeds to be further heated, triggering a temperature spike due to thermal runaway, the active battery material experiences higher temperatures than the current collectors ( Figure S9b). The homogenous average thermodynamic properties of the active battery material assumed in this model may affect the temperature distribution, as developing a layer-by-layer model could provide a different temperature profile at the cross-section of the battery that can hinder the speed of occurrence of thermal runaway. Furthermore, the effect of the series of exothermic reactions increases the battery temperature above that of the oven temperature and begins heat transfer to the current collector.
The simulated temperature profile of the thermal runaway by oven test is qualitatively analogous to the works of Hatchard et al. [57]. Compared to their works, this study has different quantitative values, such as the time of temperature spike and maximum cell temperature. Differences arose due to the differences in (1) [19] and be observed in the developed simulation. However, thermal runaway may not happen even at 423.15 K, depending on battery geometry [57].
The temperature the battery is exposed to plays a critical factor in determining the occurrence and severity of the exothermic reactions. Figure 7 shows the comparison between the dimensionless chemical parameters representing the degree of these reactions during the oven tests at 403.15 and 423.15 K. The SEI decomposition reaction degree is characterized by c SEI in Equation (30). c SEI , which represents the dimensionless amount of Li-containing meta-stable species in the SEI, had an initial value of 0.15 [34]. Liu et al. proposed that the SEI decomposes exothermically at 363. 15-393.15 K [34]. On the one hand, the SEI decomposition during the 403.15 K oven test started at t = 100 s, became more defined at around t = 400 s, and was completed around t = 1200 s. On the other hand, during the 423.15 K oven test, SEI decomposition started earlier around t = 80 s, became profound at t = 200 s, and ended earlier around t = 600 s. The time of occurrence of the SEI decomposition and the corresponding rate of reaction depend on the oven temperature.
(3) employed thermodynamic parameters. Despite being quant qualitative similarities are present in which a temperature spike indic thermal runaway. Experimental literature findings, such as the polym down at 403.15 K [25], have been reflected by the presence of SEI deco runaway can occur in LiPo batteries at 423.15 K [19] and be observe simulation. However, thermal runaway may not happen even at 423. battery geometry [57].
The temperature the battery is exposed to plays a critical factor occurrence and severity of the exothermic reactions. Figure 7 show between the dimensionless chemical parameters representing th reactions during the oven tests at 403.15 and 423.15 K. The SEI deco degree is characterized by in Equation (30). , which represent amount of Li-containing meta-stable species in the SEI, had an initia Liu et al. proposed that the SEI decomposes exothermically at 363.15-3 one hand, the SEI decomposition during the 403.15 K oven test started more defined at around t = 400 s, and was completed around t = 1200 s during the 423.15 K oven test, SEI decomposition started earlier arou profound at t = 200 s, and ended earlier around t = 600 s. The time of o decomposition and the corresponding rate of reaction depend on the  The anode-electrolyte reaction occurs as the reaction proceeds to temperatures above the 393.15 K threshold of the SEI decomposition. The dimensionless c a parameter represents the amount of Li intercalated within the graphite carbon with an initial value set at 0.75.
The anode-electrolyte reaction during the 403.15 K oven test began around t = 175 s and proceeded sluggishly. Meanwhile, the 423.15 K oven test showed the anode-electrolyte reaction earlier, at around t = 135 s and declined slowly thereafter. The value of c a drastically decreased as t approached 520 s, dropped quickly at t = 650 s, and ultimately reached a zero value at t = 700 s. This rapid decrease is due to the temperature spike as observed in Figure 6, where the anode-electrolyte reaction accelerates due to the high battery temperature.
The cathode-electrolyte reaction is represented by α c−e , which is the conversion degree of the reaction. The α c−e has an initial value of 0.04. Figure 7 shows the evolution of the value of α c−e . The cathode-electrolyte reaction during the 403.15 K oven test started around t = 310 s and then increased slowly. A considerable increase was observed around t = 650 s, where the corresponding battery temperature was approximately 438.15 K. The conversion degree reached unity, indicating the completeness of the reaction during the 423.15 K oven test at t = 700 s, where the temperature peaked as seen in Figure 6. The cathode-electrolyte reaction accelerated as the temperature was increased, especially when a threshold temperature was met.
The electrolyte decomposition reaction, which occurs when the temperature of the cell is about 473.15 K, is represented by the dimensionless concentration of electrolytes, c e , with an initial value of unity. Electrolyte decomposition did not occur for the 403.15 K oven test as the cell temperature did not meet the threshold criterion. A sharp decline in c e was observed for the 423.15 K oven test at t = 700 s, where the temperature spike occurred. The elevated battery temperature triggered the electrolyte decomposition reaction and exhausted the reaction quickly for about 50 s. The 403.15 K oven test demonstrated two immediate exothermic reactions, namely the SEI decomposition and anode-electrolyte reactions. On one hand, the cathode-electrolyte reaction for the 403.15 K oven test was extremely sluggish and is negligible. On the other hand, the 423.15 K oven test exhibited all four thermal runaway exothermic reactions. The anode-electrolyte, cathode-electrolyte, and electrolyte decomposition reactions occurred intensely. Extreme changes led to an increased spontaneous heat accumulation, ultimately spiking the cell temperature, and indicating thermal runaway. Figure 8a,b illustrate the exothermic reactions for the 403.15 and 423.15 K oven tests, respectively. Figure 8b has a superimposed graph in which the full range of heat generation is presented.

OR PEER REVIEW 17 of 24
The cathode-electrolyte reaction is represented by , which is the conversion degree of the reaction. The has an initial value of 0.04. Figure 7 shows the evolution of the value of . The cathode-electrolyte reaction during the 403.15 K oven test started around t = 310 s and then increased slowly. A considerable increase was observed around t = 650 s, where the corresponding battery temperature was approximately 438.15 K. The conversion degree reached unity, indicating the completeness of the reaction during the 423.15 K oven test at t = 700 s, where the temperature peaked as seen in Figure 6. The cathode-electrolyte reaction accelerated as the temperature was increased, especially when a threshold temperature was met.
The electrolyte decomposition reaction, which occurs when the temperature of the cell is about 473.15 K, is represented by the dimensionless concentration of electrolytes, , with an initial value of unity. Electrolyte decomposition did not occur for the 403.15 K oven test as the cell temperature did not meet the threshold criterion. A sharp decline in was observed for the 423.15 K oven test at t = 700 s, where the temperature spike occurred. The elevated battery temperature triggered the electrolyte decomposition reaction and exhausted the reaction quickly for about 50 s.
The 403.15 K oven test demonstrated two immediate exothermic reactions, namely the SEI decomposition and anode-electrolyte reactions. On one hand, the cathodeelectrolyte reaction for the 403.15 K oven test was extremely sluggish and is negligible. On the other hand, the 423.15 K oven test exhibited all four thermal runaway exothermic reactions. The anode-electrolyte, cathode-electrolyte, and electrolyte decomposition reactions occurred intensely. Extreme changes led to an increased spontaneous heat accumulation, ultimately spiking the cell temperature, and indicating thermal runaway. Figure 8a,b illustrate the exothermic reactions for the 403.15 and 423.15 K oven tests, respectively. Figure 8b has a superimposed graph in which the full range of heat generation is presented. During the 403.15 K oven test, the SEI decomposition reaction and anode-electrolyte reaction mainly contributed to the heat generated in the cell seen in Figure 8a. The SEI decomposition proceeded and ended earlier than the anode-electrolyte reaction. The cathode-electrolyte reaction occurred at a sluggish pace of heat generation, while the electrolyte decomposition reaction did not proceed. The occurrence and severity of the exothermic reactions are analogous to the behavior of their respective dimensionless During the 403.15 K oven test, the SEI decomposition reaction and anode-electrolyte reaction mainly contributed to the heat generated in the cell seen in Figure 8a. The SEI decomposition proceeded and ended earlier than the anode-electrolyte reaction. The cathode-electrolyte reaction occurred at a sluggish pace of heat generation, while the electrolyte decomposition reaction did not proceed. The occurrence and severity of the exothermic reactions are analogous to the behavior of their respective dimensionless chemical parameters or reaction degree.
All four exothermic reactions occurred during the 423.15 K oven test, as seen in Figure 8b. Qualitatively, the SEI decomposition reaction was parallel with the results from the 403.15 K oven test with a slightly longer initial heat generation rise. The two electrode-electrolyte and electrolyte decomposition reactions occurred far more intensely. The cathode-electrolyte reaction gave off the maximum heat generation rate for all reactions at 5.7 × 10 8 kW m −3 , which was negligible in the 403.15 K oven test. Similar extreme heat generation was determined by Shelke et al., where a heat generation of about 10 7 kW m −3 was experienced for an NMC battery undergoing thermal runaway at around 573 K [63]. The same occurrence was exhibited by the model of Zhang et al., where maximum heat generation (3 × 10 5 kW m −3 ) due to cathode-electrolyte reaction was achieved at 500 K [20]. The generated model exhibited higher heat generation primarily due to the higher temperature of thermal runaway occurrence and due to the different battery chemistry investigated.
The cathode-electrolyte reaction was the primary heat generation source of the thermal runaway. Selecting the appropriate cathode material with excellent thermal stability is critical to reducing the risk of thermal runaway, thus improving battery safety.

Conclusions
A multiphysics model was successfully developed for an LCO-graphite-PVdF-HFP pouch LiPo battery capturing the electrochemical, thermal, and thermal runaway behavior using COMSOL Multiphysics. The simulated model consisted of a P2D electrochemical model and a 3D thermal and thermal runaway model and was compared using experimental results using a commercial GEB 585460 LiPo battery.
Simulated electrochemical behavior was investigated and validated, showing competitive accuracy with a low RMSE of 0.06 V for 1 C discharge. The electrochemical model can be quickly adapted to different battery chemistries by simply modifying the physical material properties. Additionally, electrochemical behavior from geometric changes, such as electrode thickness, can be rapidly simulated and utilized to optimize battery design. Simulation results showed a two-step discharge phenomenon for the battery but were not apparent in the experimental cell. Working voltage plateau regions were shortened as the C-rate increased, demonstrated in both the simulation and experiment, due to higher polarization at higher discharge currents.
Simulated thermal behavior during discharge exhibited a viable RMSE of 1.5 K at 1 C discharge with an average heat generation of 0.90 W m −3 . However, the model showed appreciable qualitative discrepancies from experimental results. A buffering of temperature rise was present in both the experiment and the simulation. The increase in the temperature gradient of the battery and sink increased the heat dissipation by convection and hindered temperature buildup. Another factor that affected the phenomenon was heat generation. The reversible and irreversible battery heat generation was successfully determined using the model. The simulated thermal profile showed a prominent decrease in average surface battery temperature, which was absent in the experiment. The temperature lowering by heat absorption was primarily due to the reversible heat generation, where the endothermic values were present in the simulation. The reversible heat is a vital function of the entropic change coefficient intrinsic in nature to the electrode material. As the electrochemical and thermal models were coupled by temperature and heat generation, both models' veracity heavily relies on the accuracy of parameters determined in the literature used by the study.
The simulated thermal runaway model indicated that thermal abuse at 403.15 K in the oven test does not initiate thermal runaway and heavily relies on heat generation by SEI decomposition and anode-electrolyte reaction. Thermal abuse during the 423.15 K oven test triggered runaway and exhibited all four exothermic reactions in which the cathode-electrolyte reaction is the primary heat source reaching a heat generation of 5.7 × 10 8 kW m −3 . This entails the careful selection of cathode material to minimize the likelihood of thermal runaway is necessary.
The developed multiphysics model can serve as an initial bird's-eye view of predicting battery behavior at various operational conditions. It is recommended to explore extended numerical simulations situated in thermal runaway, such as short circuit analysis, which explores the heat generation from cell puncture and separator breakdown, mechanical abuse that changes the volumetric configuration of the cell during thermal runaway, thus altering the heat transfer mechanism, a full 3D geometric model of the battery thus including layering effects in pouch-type and jelly-roll cells, a gas-pressure buildup that considers the venting mechanism leading to additional heat release conditions, and fire and explosion models that consider the magnitude of sudden pressure burst. Furthermore, developing a parametric study and a sensitivity analysis of the parameters used in the models is recommended to extend the optimization of battery design and enhance the reliability and quality of the model.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/en16062642/s1, Figure S1: Lithium-ion battery charge and discharge mechanism for an LCO-graphite battery cell illustrate the lithium ions' intercalation and deintercalation.; Figure S2: Schematic diagram of the P2D electrochemical model.; Figure S3: Schematic diagram of the study.; Figure S4: Schematic illustration of the thermocouple positioning at one side of the battery during thermal and thermal runaway tests.; Figure S5: An illustration of the geometrical meshes of the multiphysics models used in the study.; Figure S6: Simulated results of 10 C discharge rate.; Figure S7: Simultaneous analysis of simulated working potential and average surface battery temperature during 1.0 C rate discharge.; Figure S8: Simulation and experimental average battery surface temperature vs. time during 443.15 K oven test; Figure S9: Three-dimensional temperature distribution simulation of 423.15 K oven test.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.