Improved internal short circuit models for thermal runaway simulations in lithium-ion batteries Open Access

Dwindling conventional energy resources and increasing awareness regarding the environmental impact of their emissions have pushed for action to decarbonize the energy economy. Significant improvements have occurred on the supplier as well as on the consumer ends. According to a report by the United States Environmental Protection Agency (USEPA), the transportation sector contributes to 28% of the total greenhouse gas emissions.¹ Hence, it is no wonder that recent years have seen an unprecedented number of electric vehicles (EVs) on the road.² The boom in the EV market has been a result of the progress in research and widespread adoption of lithium-ion battery (LIB) technology. These batteries have higher operational voltage (3.6 V), higher energy density, and great cycling performance, which have enabled the manufacturers to design battery packs with fewer cells and greater reliability.³ Since a long time after their conception in the late 1800s,⁴ the EV is finally ready to compete with the gasoline-based automobile on an equal footing.

LIBs suffer from a critical handicap known as thermal runaway (TR).⁵ Above a critical temperature, they start self-heating due to the initiation of exothermic decomposition chain reactions.⁶ If left unchecked, the battery temperature can shoot up to 800–1000 °C in a matter of seconds.⁷ Hence, it is a critical safety concern for their use in regular transportation vehicles. Often, the primary trigger for thermal runaway is the internal short circuit (ISC) of these cells, which can be caused by many factors, ranging from mechanical and thermal to electrical malfunctions.⁶ 

A dedicated effort toward understanding ISC phenomena in LIBs was first put forth by Santhanagopalan et al.⁸ They coupled electrochemical and thermal models with thermal runaway reactions and studied the characteristics across multiple parameters. They identified that aluminum-anode short led to the highest cell temperature inside the cell among all possible short circuit scenarios. This was attributed to the low electrical resistivity of the anode material, low activation energy of the anode decomposition reactions, and poor heat dissipation mechanisms. In the same year, Maleki and Howard induced ISC using three methods: indentation, small nail penetration, and cell pinching.⁹ They compared their experimental results with the numerical results from a thermal model. For simplicity, they skipped the detailed modeling of electrochemistry and chemical decomposition inside the cell and used the experimental values for heat generated by electrical, electrochemical, and chemical decomposition reactions. Zavalis et al. developed a 2D coupled electrochemical thermal model in prismatic Li-ion cells to simulate three different short circuit scenarios.¹⁰ They highlighted that the mass transport of Li ion particles in electrolyte is the most important factor that determines the rate of temperature increase.

Chiu et al. developed an extensive model coupling 3D thermal and nail penetration models with a 1D electrochemical model.¹¹ Chen et al. presented a unique approach toward ISC modeling using a network of electrical resistances between multiple nodes in the cell.¹² Although the model makes several simplifying assumptions, the numerical results seem to agree well with the experiments conducted by the authors. Zhang et al. developed the first coupled mechanical–electrical–thermal model for simulating the mechanical abuse induced short circuit¹³ and followed up by developing a multiscale coupled model to predict the electrochemical–thermal response of the cell during an evolving mechanical failure with a simultaneous internal short circuit.¹⁴ As an example of another unique modeling approach, Liu et al. used equivalent circuit models to simulate the electrochemical characteristics of a large format 20 Ah Li-ion cell and coupled it with a thermal model.¹⁵ They used a relatively simplistic constant electrical resistance approach to simulate the ISC and validated its simulation results with experimental studies in the literature.

Fang et al. developed another popular model that has been widely adopted in the literature. They extended the existing general purpose coupled electrochemical–thermal model to replicate the anode-aluminum and cathode–anode ISC scenarios and validated it with experiments.¹⁶ Much research efforts were also conducted at NREL as indicated by the reports of Kim et al.¹⁷ and Keyser et al.¹⁸ They developed multi-physics ISC models and analyzed the thermal performance of the cell as a function of short resistance, cell capacity, etc. Zhao et al. conducted numerical investigations on a large format Li-ion cell subjected to nail penetration.¹⁹ They focused on a single ISC scenario, aluminum–copper, but used multiscale coupled models to simulate a cell having commercial dimensions. They simplified their computations by assuming a constant short resistance between the aluminum and copper layers and hence not solving the solid potential equation within the nail body. Their major finding was the identification of different heating scenarios, i.e., localized vs global, as a function of short resistance and cell resistance. Their conclusions matched the results of the NREL studies and have been accepted by the research community since then.^19,20

Kim et al. published fundamental work in understanding the transport processes during ISC.²⁰ They modeled a Cu–Al short circuit for varying the values of short circuit electrical resistances. They found out that diffusion in the solid active material plays a significant role at high shorting currents. The concentration depletes rapidly in the electrolyte and accumulates on the electrode surface. This leads to non-uniform transport properties, resulting in higher heat generation rates, and a shift in the heating regime from localized to global.

Recently, Zhang et al. also conducted multiscale simulations by coupling a 2D electrochemical–thermal model with a large-scale TR model for a Li-ion pouch cell.²¹ They provided a safety regime map plotting the volumetric heat source density with time and marking different zones on the chart to indicate the areas vulnerable to catastrophic TR. They achieve the combined analysis of ISC and TR in a “coupled-decoupled manner,” i.e., the thermal runaway kinetics equations are solved for the larger 3D cell model and the electrochemical equations are solved for the micro-scale model. The temperature and heat release densities are then exchanged between these models to couple them together.

It should be noted here that ISC is not the only cause of TR in LIBs. Studies have been reported where TR occurred without ISC.^22,23 The primary mechanism in such cases is the reaction of oxygen produced at the cathode with the anode. Ren et al. conducted accelerated rate calorimetry (ARC) tests on 24 Ah commercial pouch Li-ion cells.²⁴ They used thermally stable separators to avoid separator shrinkage and hence ISC. They observed that the battery electrical resistance increased sharply before TR due to the depletion of electrolyte, porosity reduction of electrodes, and shutdown of the separator. In that case, the heat generated by ISC had a minimal contribution as compared to the chemical degradation reactions. All simulations carried out herein are different from the above-mentioned scenario. We focus on a particular case where ISC is triggered by nail penetration. The current always has a low resistance path to travel between the electrodes, i.e., the penetrating filament. In these simulations, the ISC is the cause for TR and not the effect.

The primary aim of this article is to develop a low-cost alternative to the coupled electrochemical–thermal model for Li-ion batteries. We found a shortage of empirical ISC models having a sound physical foundation. There is no shortage of numerical models available in the literature, as discussed above, but almost all decouple the length scales involved in the problem to replicate large format commercial cells and save computation time in doing so. Although this approach leads to the modeling of geometries with practical applications, we suspect that it comes at the cost of losing fundamental information regarding the mechanisms of ISC and heat transfer happening at smaller length scales. We address this problem by dividing it into two sections. In the first section, we extensively simulate a small test cell, while keeping model-simplification assumptions at a minimum. Because the test cell is small enough, we can solve the coupled electrochemical–thermal equations along with the thermal runaway chemical kinetics equations at each mesh point in a reasonable computational time. Most studies found in the literature either assume a range of short resistances,^19,20,25 calculate the short resistances using experimental voltage and current profiles,²⁶ or skip coupling the TR reactions altogether.²⁶ To the best of our knowledge, this is the first study to analyze the current pathways in a completely coupled electrochemical–thermal model with an integrated TR chemical kinetics solver, in a 2D geometry. In the second section, we use these fundamental results to develop an empirical model. We attempt to bridge the gap between the fundamental simulation studies and the commercial application-based modeling approaches, by developing low cost ISC heat generation models, keeping the physics of the phenomena as intact as possible.

An axisymmetric disk-shaped cell is used as the test geometry for all parametric studies. The cell consists of five stacked layers: the positive current collector, the positive active material (cathode), the separator, the negative active material (anode), and the negative current collector. The individual material properties and geometric parameters are listed in Table I, where T is the local temperature, c is the species concentration, F is the Faraday constant, soc is the state of charge of the cell, and R is universal gas constant.

To induce the electrical short, a penetrating filament in the shape of a cylinder is modeled to connect the respective layers, depending on the type of short. The electrical conductivity of the filament material is initially set to 0. Within a very short time interval (from t = 0 s to t = 10⁻³ s), the electrical conductivity is ramped up to the actual physical value. This method is very commonly used in the literature to simulate the initiation of the internal short circuit.^21,26

The aim of this study is to understand the short circuit phenomena during TR from an applied thermal engineering perspective. Hence, the following assumptions are made to ease the computational complexity:

1. The cell is axisymmetric and the individual material properties are homogeneous.

2. Complete interfacial contact is assumed throughout TR and short circuit phenomena.

3. Electrode film resistances and contact resistances have been ignored.

4. Thermophysical and electrical properties are constant during decomposition reactions.

5. Many studies in the literature use a constant terminal voltage boundary condition. This is a good approximation to get the initial value of the short circuit current. However, to simulate the extended characteristics, we use a zero current boundary condition at the terminals. This is assuming that the short provides a path of much lower electrical resistance for the current, and hence, very negligible current flows in the external circuit.

6. No venting or combustion occurs throughout the course of the short circuit event.

It is coupled with the thermal runaway and electrochemical module by the terms Q_ISC and Q_TR, which represent the volumetric heat generation due to ISC and TR, respectively. ρ, C_P, and κ are the local density, specific heat, and thermal conductivity of the material, respectively. A constant heat transfer coefficient (HTC = 5 W/(m²K)) boundary condition, with an ambient temperature of 298 K, is applied to the radial edge of the cell, as shown in Fig. 1.

The electrochemical model is based on the porous-electrode theory developed by Newman and Tiedemann.³⁶ This is a popular model, which has been implemented extensively in the literature in 1D,^10,19 2D,²¹ and 3D¹⁶ forms. It has also been validated at a range of working temperatures and C-rates by various authors.^37,38 Hence, this section is only aimed to be a brief overview of the model. Readers are suggested to refer to the source texts for specific details. The validation of the replicated model developed in this study with data available in the literature³⁷ can be found in Fig. S3 of the supplementary material. At the highest discharge rate used for validation, the error between the numerical model developed in this study and the experimental data is less than 8%. The main reason for the discrepancies is that the calibration parameters used in Ref. 37, such as electrode film resistance, contact resistance, and cathode exchange current density, have either been neglected or not fine-tuned at all. This approach has been adopted to keep the model as general as possible and avoid overfitting to a particular set of experimental data. Another possible reason for the error can be attributed to the differences in source references used for physical properties.

There are two primary phases in a Li-ion cell: the electrolyte phase and the solid phase. Charge transfer and mass transfer, due to the transport of electrons and Li-ions, are continuously occurring across both phases. Hence, charge conservation and mass conservation across each phase are the primary modeling equations. Figure 1(c) illustrates all charge and mass transfer interactions considered in the electrochemical models. Apart from these four equations, there is a fifth equation that defines the electrochemical kinetics relating the electrode current density to the activation over potential. Hence, there are five unknown variables—electrical potential in the solid phase (⁠$ϕs$⁠, electrical potential in the electrolyte phase (⁠$ϕl$⁠, insertion particle concentration in the solid phase (⁠$cs$⁠, salt concentration in the electrolyte phase (⁠$cl$⁠, and electrode current density (⁠$is$⁠.

Table III lists the values of parameters used for each reaction. A MATLAB code implementing this model for a lumped cell is included in the supplementary material. The reader is also referred to the original publication for elaborate understanding. The thermal runaway model is validated by replicating the lumped system model used in the parent publication, in COMSOL (see Fig. S4 of the supplementary material).

COMSOL Multiphysics is used for solving the coupled equations. The solution is computed in two steps: (1) current distribution initialization and (2) time dependent solution. For the stationary step, the solver choice is set to automatic in the COMSOL toolbox. For the second step, the variables are segregated into groups to improve numerical convergence. For all these different segregated groups, MUMPS is used as the linear solver, and the nonlinear method is set to “automatic highly nonlinear (Newton).” Because the coupled numerical equations are extremely difficult to solve even with advanced nonlinear solvers available in the toolbox, the relative tolerance for the time dependent step is set to 0.005. Depending on a case-to-case basis, the minimum damping factor is also adjusted to achieve numerical convergence within the specified tolerance limits. Figure S5 of the supplementary material illustrates the mesh independence of the current model for a particular case.

The effect of three parameters was analyzed on the electrical and thermal characteristics of the cell undergoing ISC. The parameters include the type of ISC, nail radius, and cell radius. While carrying out the studies on a particular parameter, the others are fixed at nominal values specified in the relevant sections.

In this study, we have simulated three different ISC scenarios: cathode–anode, aluminum-anode, and copper cathode ISC. Figure 2 illustrates these scenarios graphically. The cell radius was fixed at 0.72 mm, and the nail diameter was fixed at 2 μm. As shown in Fig. 2(a), the maximum temperature attained in the Al-anode ISC (∼1180°C) is more than one order of magnitude higher than the other two scenarios (∼70°C). Cathode–anode and Cu-cathode ISC scenarios behave like each other with less than 2°C difference between their peak temperatures. The severity of Al-anode ISC can also be gauged by the fact that it takes only 100 s to discharge the cell completely, whereas in the other two scenarios, the cell takes more than 2000 s to exhaust its stored electrochemical energy. The TR phenomenon involves the release of stored chemical energy inside the constituent materials of the cell. This stored chemical energy should be the same for the cells with the same amount of active materials, electrolytes, etc. However, the share of this energy that gets converted into heat during TR depends on the mechanism via which the decomposition reactions, mentioned in Table III, are triggered. The severity of the decomposition reactions increases immensely if the triggering mechanism imparts a higher temperature increase rate to the system because the rate of reactions is a function of temperature itself.

For case (a), the trigger temperature is T = 100 °C, while for case (b), the trigger temperature is T = 40 °C. Once the cell hits the trigger temperature, the decomposition reaction will start and provide additional self-heating, which typically is so high that it overpowers the heat dissipation rate and the temperature shoots up, resulting in TR. In case (b), the temperature never exceeds 40 °C and the decomposition reaction never starts. This is illustrated in Fig. 3 for better comprehension.

Hence, when comparing the severity of an ISC scenario, it is beneficial to look at the heat generated solely due to the electrochemical processes [Fig. 2(b)]. As expected, the Al-anode ISC has a heat generation rate, which is almost 35 times higher than the other two scenarios. If one integrates the three power curves shown in Fig. 2(b) with respect to time, to calculate the total electrochemical heat generated during ISC, one would get the same answer for all three cases because the cell is completely discharged in all the three cases. However, only one of them has a heat generation rate high enough to increase the cell temperature beyond the trigger temperature of the first few decomposition reactions.

However, why does the Al-anode ISC succeed in generating catastrophic heat generation rates where other two scenarios fail? The answer lies in the electrical pathways that the short circuit current has available at its disposal. Figures 2(c)–2(e) show the current, voltage, and resistance of the three ISC scenarios, respectively. It is not surprising that the electrical resistance in a particular scenario remains almost unchanged with time, given the assumption that material properties are not treated independent of temperature and reaction progression. The electrical conductivity of cathode, i.e., NCM111, is almost 25 times less than that of anode, i.e., lithiated graphite. Hence, the pathway in which the current can skip this bottleneck should have exceptionally lower electrical resistance. This can be seen in Fig. 2(e), where the Al-anode ISC has an electrical resistance of 2.4 Ω as compared to the 87 Ω in the Cu-cathode ISC and 95 Ω in the cathode–anode ISC. Hence, the lower electrical resistance of the Al-anode pathway leads to a much higher electrochemical heat generation rate, which is the primary mechanism for triggering thermal runaway.

Note that there is a fourth type of ISC, which we have intentionally omitted, where the two current collectors are connected via a penetrating filament directly (hence the name Cu–Al ISC).⁸ If one were to follow the electrical pathway reasoning elucidated in the previous paragraphs, it would turn out to be the pathway with the least electrical resistance. Hence, it leads to extremely high ISC current, which discharges the cell much faster than the other scenarios. For example, the total numerical electrical resistance for the cathode–anode ISC with a 2 µm nail is 94.9 kΩ. In a Cu–Al ISC scenario, the entire short circuit current will pass through the nail, with a minimal amount of current passing through any electrode layer. The total electrical resistance in that scenario will be 0.92 kΩ. Hence, with the same initial voltage, the short circuit current in the Cu–Al ISC scenario will be ∼103 times higher than that in the cathode–anode ISC scenario. This high amount of current in the electrochemical thermal model coupled with the local thermal runaway kinetics equations makes it an extremely hard problem to solve numerically. The simulation requires very small time step size, which was beyond the computational power at our disposal. Hence, this scenario is not considered here. According to the given reasoning, the Cu–Al ISC should be the most severe TR event out of the four scenarios. However, the existing literature⁸ suggests otherwise. This apparent contradiction can be resolved by considering the spatial temperature distribution, which has been conveniently skipped in the past literature, and is considered herein.

Figure 4 shows the local power and volumetric power densities as a function of time in two ISC scenarios. During Al-anode ISC, the electrochemical power generated with a 2 μm radius nail is ∼4 times larger than that of the cathode–anode ISC with a larger 18 μm radius nail. Furthermore, the regions with the largest heat generation contributions are different in both cases. The cathode layer is the main contributor during cathode–anode ISC, whereas the anode is the dominant contributor during Al-anode ISC. This is obvious given that the cathode region is bypassed during the discharge in the latter scenario.

This explains the three orders of magnitude difference in the volumetric heat generation in the nail region between the two cases. In addition to the variance among different layers, there are considerable spatial gradients within a particular layer. Figure 5 show the contours of volumetric power density and temperature at the time of maximum temperature attained in the cell (the dotted line shows the corresponding time, t = 65 s for the Al-anode ICS and t = 135 s for the cathode–anode ISC). The spatial gradients are observed in the layer with the largest power generation. Although the order of magnitude of power densities is similar for the two scenarios, the thermal gradients are significantly higher in the case of the Al-anode ISC. The highest power density spots are typically observed in the active material layers immediately adjacent to the penetrating filament. Hence, these spots are the most susceptible to high temperatures. However, in the case of Al-anode ISC, this spot is in the anode region, whose thermal conductivity is 66% lower than that of the cathode material (Table I). Coupled with the fact that the power generated is also ∼4 times higher, the result is a higher peak temperature and higher thermal gradients.

During Cu–Al ISC, these hot spots are virtually non-existent, due to the high thermal conductivity of Al and Cu, and the minimal share of power contribution from the current collectors. Hence, despite the high ISC current, Cu–Al ISCs lead lower cell temperatures than their counterparts.

In the subsection titled Effect of ISC type, we saw that the cell does not undergo TR for cathode–anode and Cu-cathode ISC scenarios when the penetrating filament radius is 2 μm. In this section, we choose an ISC type, i.e., cathode–anode ISC, fix the cell outer radius at 0.72 mm, and vary the nail radius. Figure 6 shows the thermal and electrical characteristics of the cases where the nail radius is increased from 5 to 18 μm. According to the results, there is a critical nail radius, only beyond which the cell undergoes TR. The temperature increase rate imparted by the electrochemical heat generated when a 5 μm nail is pierced into the cell is not enough to increase the cell temperature to the critical self-heating temperature before the heat dissipation rates plateau. We know from the previous parametric study that nails smaller than this radius generate even lower heating rates.

A larger nail radius leads to a larger temperature increase rate. A 10 μm nail successfully triggers the catastrophic event, and hence, it can be concluded that for this cell configuration and geometry, the critical electrochemical heat generation rate lies between 0.4 and 0.9 mW.

The total electrochemical heat released in each case is still the same, but the higher initial power generation discharges the cell more quickly as the nail radius increases. Figures 6(d)–6(f) show the electrical characteristics for this study. A nail with a larger cross section leads to a lower electrical resistance and hence increases the initial electrical current.

An astute reader might try to fit the resistance vs nail radius curve [Fig. 6(f)] using an inverse parabolic function $(Relectrical=Lnail/σπrnail2)$⁠. However, the curves do not show a good regression fit. This is because the current follows a more complicated two-dimensional path, which cannot be captured using the simplistic series resistance assumption. This will be discussed in more detail in the subsection titled Electrical pathways and resistances.

Changing the cell radius is equivalent to keeping the same cell configuration but changing its capacity. In this section, we analyze the susceptibility of cells with different capacities toward thermal runaway. Like the subsection titled Effect of penetrating filament radius, the type of cathode–anode ISC scenario is simulated with a 15 μm radius nail.

It can be seen from Fig. 7(a) that smaller cells are more vulnerable to TR as compared to larger cells in the same configuration. This goes against conventional wisdom, given that a larger cell contains a larger amount of electrochemical heat to dissipate and raise the cell temperature. The highlighted numbers in Fig. 7(c) represent the total heat, contributed by both TR and ISC, during the discharge. A cell of 1.775 mm radius generates almost 14.5 times more heat than a cell of 0.355 mm radius and still fails to undergo TR. This can be explained using the discussion of Fig. 3. Even though the total electrochemical heat generated in a larger cell is higher, the initial rate of heat generation is still the same if the nail radius and the ISC type are kept fixed, as shown in Fig. 7(b). The smaller capacity cells also have smaller surface areas for heat dissipation. Hence, according to Eq. (15), smaller cells attain the threshold temperature more easily and undergo TR. Figure 8 illustrates this point by plotting the total heat generation rate vs the heat dissipation rate. In a larger cell, the temperature increase rate is damped due to the higher heat dissipation surface area, because the rate of electrochemical heat generation does not increase with the increase in cell capacity. A larger cell only ensures that the discharge takes longer, hence reducing the susceptibility of the cell toward higher temperature rates and hence TR.

The key idea from this parametric study is that electrochemical heat generation is independent of cell capacity. Figure 7(f) shows that this can be attributed to the independence of electrical resistance of ISC from the cell capacity. Intuitively, it makes sense, as the penetrating filament seems to be the bottleneck for the ISC current. Hence, electrical resistance being a strong function of nail radius and not the cell radius is justifiable. In the subsection titled Electrical pathways and resistances, we will discuss why this understanding of the phenomenon is incorrect and how can we better justify the numerical observations.

Looking at an example cathode–anode ISC and assuming that the current is traveling from the anode to the cathode via the penetrating filament, the simplest form of electrical network would be to assume the five components of the cell to be arranged in a 1D series network. This simplification assumes that there is a homogeneous current density in the z direction in all layers as shown in Fig. 9(a).

The electrical resistance of each layer and the total electrical resistance in this simplified model, for the two cathode–anode ISC parametric studies, are shown in Tables IV and V, respectively. It can be readily seen that neither the absolute value nor the trend matches the calculated and numerical values of electrical resistances in the two parametric studies. Table I lists the electrical conductivity of the individual materials that are involved in the simulation. Based on those values and the individual layer resistances in Tables IV and V, the bottleneck electrical resistance is the cathode layer as it has the highest electrical resistance. Hence, it should ideally dominate the overall electrical characteristics. However, from the parametric studies, we observe a completely opposite trend [Figs. 6(f) and 7(f)]. When we double the nail radius from 5 to 10 μm, the theoretical change in overall network resistance should be 0.3%, whereas simulations indicate a resistance drop of 56%. Similarly, when the cell radius is doubled from 0.72 to 1.44 mm, the theoretical change in overall network resistance should be 74%, but the simulations indicate a change of 1%. It is obvious that our simplified series electrical resistance network is not appropriate to model the given problem.

Figure 9(b) shows the contours of the z-component of electrical current density for two different nail radii. The most striking observation from the contours is that a large fraction of the cathode and anode layers do not contribute to the vertical conduction of the ISC current. Most of the current traveling vertically in the active layers is constrained within a small region whose radius is slightly larger than the radius of the penetrating filament. Table VI lists the calculated values of effective radius of this region for the NCM111 layer using the numerically computed overall electrical resistance values. It turns out to be slightly larger than the nail radius and tends to approach it as the nail becomes larger. It should be noted that the nail is not the dominant electrical resistance, but its dimensions affect the cross section area of the cathode and anode layers that get effectively “utilized” for vertical current conduction. Modeling the 2D electrochemistry during the ISC event helps us identify the key electrical pathways and resistance involved during the process and enables us to create more accurate reduced order electrical resistance networks.

The MATLAB code used for modeling is included in the supplementary material. We observe that our ISC model shows a Mean Absolute Percentage Error (MAPE) of 45.7% as compared to the past ISC model, which has an MAPE of 105.9% when validated against their experiments. Hence, our empirical ISC model reduces the MAPE by ∼60%. This comparison also acts as a qualitative experimental validation for the coupled numerical simulations reported in this article.

The best use case of the developed model is to simulate the ISC scenarios in lumped system models. In this subsection, we have calibrated the model for the nail penetration scenario for a 24 Ah NCM111-graphite pouch cell. It is an empirical model, which captures the general trend of heat generation during an ISC even in LIBs. Hence, a model calibrated for a particular battery cannot be applied directly to a different one. Any change in the following parameters demands a recalibration of the model: cell capacity, cell geometry (pouch, cylindrical, prismatic, etc.,), cell chemistry (change in anode, cathode, electrolyte, and separator materials). The procedure to calibrate the model is relatively simple and can be done using the genetic algorithm toolkit available in the MATLAB. The MAE between the model outcome and the experimental data for the test cell should be used as the objective function. A vector consisting of four model parameters (K₁, m, t_rel, and K₂) is the optimization variable. Once calibrated, the ISC model can be integrated with three-dimensional coupled models for advanced use cases. Readers should note that the current study is specifically catered to ISCs in LIBs and hence should not be extended to other battery types. The mechanism for thermal runaway and internal short circuit in Solid State Batteries (SSBs) is an active field of research.^40,41 Future work on this study will focus on the development of the ISC models for the same.

In this study, we have tried to bridge the gap between the fundamental understanding of ISC in Li-ion batteries and its practical implementation in battery thermal management applications. We have developed a coupled electrochemical–thermal–TR chemical kinetics model in COMSOL to understand the effect of different parameters on the electrochemical heat generation profile during an ISC plus TR event. The key takeaways are as follows:

1. Among all the ISC scenarios, the Al-anode ISC scenario generates the largest power and leads to the highest cell temperatures. This is due to the elimination of the electrical bottleneck, cathode, from the current pathway.

2. The severity of TR increases as the penetrating filament increases in diameter. This is because of the decrease in the overall electrical resistance of the ISC event, which leads to higher temperature increase rates.

3. Smaller cells are more susceptible to TR as compared to larger cells due to a smaller surface area for heat dissipation. If the penetrating filament radius and the type of ISC are kept same, the electrochemical power generation magnitude remains independent of the cell capacity. Hence, a larger cell easily dissipates the heat generated and never reaches the critical temperature required to initiate TR.

4. The ISC current is not homogenously distributed in the entire cell. There is only a small region in both active material layers that contributes to the total electrical resistance of the ISC event.

5. The proposed empirical model developed using the results from the numerical simulation reduces the MAPE between the numerical and experimental temperature rates for lumped system models by 60% as compared to the state-of-the-art ISC models.

This study provides an improved and effective low-cost alternative for electrochemical heat generation models. Researchers can use our model in the thermal modeling of larger and more complex systems to accurately capture the physics of ISC, without developing coupled models. Hence, our work can save immense computational costs for design optimization studies in battery thermal management systems.

Additional figures including interpolation functions for the physical properties of materials, experimental validation of electrochemical and thermal runaway models, and mesh independence studies can be found in the supplementary material. A reduced order MATLAB code for simulating thermal runaway in a lumped system LIB model has also been included in the same.

The authors acknowledge the funding support from the Power Optimization of Electro-Thermal Systems (POETS) National Science Foundation Engineering Research Center with Cooperative Agreement No. EEC-1449548 and the Advanced Research Project Agency-Energy (ARPA-E) with Cooperative Agreement No. DE-AR0000895. N.M. acknowledges the funding support from the International Institute for Carbon Neutral Energy Research (No. WPI–I2CNER), sponsored by the Japanese Ministry of Education, Culture, Sports, Science and Technology.

The authors have no conflicts to disclose.

Bakhshish Preet Singh: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead). Yashraj Gurumukhi: Methodology (equal); Resources (equal); Software (equal); Writing – review & editing (equal). Engin Tekinalp: Data curation (equal); Investigation (equal); Methodology (equal). Hao Wu: Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Myung Ki Sung: Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). Nenad Miljkovic: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.