Benchmarking Between COMSOL and GPYRO in Predicting Self-Heating Ignition of Lithium-Ion Batteries

Recent studies have shown that self-heating ignition is a possible cause of fires when Lithium-ion batteries (LIBs) are stacked in large numbers, for example, during storage. The understanding of this ignition type is limited, and most current studies are based on numerical modelling. The different modelling tools found in the literature differ in their assumptions, capabilities, and resources needed, and may provide significantly different predictions. This study presents a benchmarking between COMSOL Multiphysics, which is one of the most prevailing tools used in modelling thermal-electrochemical behaviour of LIBs, and Gpyro, which is widely used in modelling ignition of solid fuels. Four case studies are designed with increasing levels of complexity: (1) just chemical kinetics at the microscale, (2) just heat transfer at the mesoscale, (3) self-heating behaviour at the mesoscale for coupled chemical reactions and heat transfer of a single cell, and (4) four-cell ensemble for multiphysics at a larger scale. The results of scenarios #3 and #4 are also compared to experiments. The results show that although COMSOL and Gpyro have significant differences in their assumptions and resources needed, both tools can accurately predict the critical conditions for ignition for self-heating, which validates their use to study the safety of LIBs.

W Material content kg m -3 C Electrochemical capacity Ah T Temperature K c p Specific heat capacity J kg -1 K -1 c Concentration (dimensionless) f Mass fraction function n Reaction order h Specific enthalpy J kg -1 h Convective heat transfer coefficient W m -2 K -1 m 00 Mass flux kg m -2 s -1 p Pressure Pa k Effective thermal conductivity W m -1 K -1 q 000 Volumetric heat generation W m -3 q 00 Heat flux W m -2 t Times Y Mass fraction z SEI thickness (dimensionless) Greek Symbols

Introduction
Lithium-ion batteries (LIBs) are an important energy storage technique for mobility such as consumer electronics, electric vehicles (EVs), and even smart grids [1][2][3]. As an energy storage technology, fire safety is also important for the LIB industry, especially for large-scale applications, which may result in more severe consequences. While the energy density of LIB increases, safe storage and transport become big challenges to the industry [4]. Several major fires related to LIB storage and transport [5][6][7] have been reported in the past decade, while the understanding of the potential causes remains limited. The industry and battery community focuses more on chemistry and identify the cause to all kinds of abuse conditions and manufacture defects [1,8,9]. However, fires could also be initiated by another fundamental mechanism driven by heat transfer, self-heating ignition [10][11][12], which is barely discussed in the literature. Self-heating ignition is a spontaneous ignition of a reactive material resulting from internal exothermic reactions [10][11][12]. For a large ensemble of reactive materials, the heat generated inside by low-temperature exothermic processes is difficult to dissipate and may accumulate, causing self-heating and finally ignition at ambient temperature such as self-heating ignition of coal mines and biomass [10,13]. This kind of fundamental cause of fires is also possible for LIBs stacked in large sizes, but has received little attention in the literature.
Since self-heating ignition happens more often when reactive materials are stacked in large piles, it is difficult to study via large-scale experiments due to the heavy costs and large spaces needed. Current most feasible approach is through investigating the reaction kinetics by means of lab-scale experiments [14][15][16][17] and then applying it via physical theories or heat transfer models to make large-scale predictions [18]. According to chemical kinetics, models for self-heating ignition can be divided into two categories. The first type is one-step global models, which are based on classical self-heating theories (e.g., Frank-Kamenetskii's theory [19]) and assume an effective global one-step reaction to represent the major reaction studying self-heating behaviour. This simplification focuses on the dominating reactions and supports theoretical analysis and fast predictions. It has been widely used to study and predict the self-heating ignition of traditional reactive materials such as coal, shale, and biomass [13,19]. The second type is multi-step kinetics models. Some materials or a system with mixtures of different components (e.g., mixed sawdust and olive oil [19]) over a range of temperature (from ambient temperature to self-ignition) may undergo several reactions playing crucial but different roles, which requires multi-step kinetics to accurately model. Multi-step models contain more complex kinetics, which usually could not be solved analytically and thus require numerical methods to obtain simulation results.
Hatchard et al. [20] and Kim et al. [21] used three and four-step reactions to model the self-heating ignition of a single 18,650 type of LiCoO 2 (LCO) cell in oven experiments. The four-step kinetics have been used widely to model the thermal behaviour of LIBs in different shapes [22], different oven temperatures [23], and a stack of cells [24][25][26]. In previous work [18,27], we used the four-step kinetics to model the self-heating behaviour of large LIB ensembles. On the other hand, He et al. [28,29] used a one-step effective model based on Frank-Kamenetskii's theory to simulate self-heating ignition of battery ensembles with a different number of cells. Liu et al. [4,30] also used Frank-Kamenetskii's theory to model their experiments on piles of LIBs. The one-step effective models and multi-step models have been developed based on different simulation software, prevailing in fields. Different simulation platforms may have a different focus on particular phe-nomena (like heat transfer, chemical kinetics, or fires), therefore, have different simulation performances on different physics modelling. This study conducts a benchmarking analysis on two popular simulation platforms, which support numerical modelling of self-heating ignition of LIBs: COMSOL Multiphysics and Gpyro.
COMSOL is a numerical tool aiming to provide numerical solutions for a wide range of scientific and engineering problems, especially for problems cross disciplines. It contains various basic modules and interfaces to solve governing partial differential equations (PDE) based on the Finite Element Method (FEM). COM-SOL has been widely used to analyse the electrochemical and thermal behaviour of LIBs [31][32][33], and it contains basic modules for LIBs and heat transfer. Gpyro is an open-source tool designed for pyrolysis modelling of combustible solids based on the Finite Difference Method (FDM). The embedded coding covers mechanisms of mass transfer, heat transfer, and chemistry in both gas and condensed phases, and has been validated against many fire problems including selfheating ignition studies [10,28,34]. Focusing on self-heating ignition of LIBs, the two numerical tools have different computational performances since they are prevailing in a different area and may have different capabilities in simulating chemical kinetics, heat and mass transfer, and multi-physics coupling. Our group has used both COMSOL [18,26] and Gpyro [29,35] for this fire problem, and it is of interest to benchmark COMSOL and Gpyro for the predictions of self-heating ignition of LIBs, and analyse their advantages and limitations.
Four scenarios are analysed in this study with different levels of complexity: (1) kinetics at the microscale (pure chemical reactions), (2) heat transfer at the mesoscale (pure heat transfer), (3) self-heating behaviour at the mesoscale for coupled chemical reactions and heat transfer of a single cell and (4) a four-cell ensemble for Multiphysics at a larger scale. The performances of COMSOL and Gpyro for the different physics involved are compared and quantitatively analysed to assess their capability in modelling this fire problem.

Theoretical Frameworks
This section covers the basic setup for both COMSOL and Gpyro with specific consideration of self-heating ignition of LIBs.

Numerical Set Up for COMOSL
The detailed development of self-heating ignition model using COMSOL is introduced in our previous numerical studies [18,27]. This section summarizes the governing equations, chemical kinetics, and assumptions adopted by considering the special characteristics of self-heating ignition of LIBs.
Self-heating ignition focuses on the transition from a stable state to the onset of thermal runaway, with a special focus on early-stage reactions at temperatures below 200°C. The thermophysical properties (q, c p , k) are assumed to be constant in this temperature range. In this study, a battery cell is simplified as a bulk that all reactants are homogeneously distributed inside the cell. Previous studies have shown that a lump thermal model [20] already provides adequate accuracy for the temperature prediction for the single-cell compared to a 3D model [21] for the focused temperature range. Therefore, this lumped simplification is acceptable and at mean time, significantly reduces the computational costs. The governing equation could be simplified as: COMSOL supports simulations on complex geometries with 0D, 1D, 1D axisymmetric, 2D, 2D axisymmetric, and 3D modelling. It also contains various types of interfaces for boundary conditions. For self-heating ignition analysis, the commonly used thermal boundaries include: constant temperature, heat flux, convection, and thermal insulation. In regard to the chemical models used, COMSOL has a module to simulate multi-step reactions, the interaction between reactions, and species transport. The chemical reactions considered in this study include four major steps: solid electrolyte interphase decomposition (SEI decomposition), negative-electrolyte reaction, positive-electrolyte reaction, and electrolyte decomposition [36]. These four reactions are the sources considered for internal heat generation. The heat generated by the four reactions is given by (Tables 1 and 2): One of the most appealing features of COMSOL Multiphysics is its flexibility. While the software provides basic interfaces for most of the physical problems, it Table 1 The kinetic model for LiCoO 2 Lithium-Ion Battery [20,21] Chemical reactions also has easy access for users to make their own variable definitions and self-coding, even changing the governing equations. This allows users to adjust and develop models for the specific problems.

Numerical Set Up for Gpyro
Gpyro is a open-source tool designed for pyrolysis modelling for combustible solids. The embedded code in Gpyro considers a bulk of porous reactive material. The governing equations include the conservation equations [28,37] The original equations in Gpyro consider a 1-D problem. The later version of Gpyro supports direct simulation of 3D problems with simple shapes. However, the 3D simulations extensively increase its computational costs and have convergence issues. In this study, we adopted Pseudo 3D simplification by incorporating the heat transfer from the other dimensions into a source term in the 1D governing equation to consider 3D heat transfer. For example, Eq. (16), h vl [28,38] is the volumetric heat transfer coefficient considering the heat transfer from x and y directions. Gpyro also has various interfaces to set thermal boundaries such as constant temperature, heat flux, and, convection, and thermal insulation. The reactions scheme considered in Gpyro is similar to the classical self-heating theories which assume an effective overall Arrhenius type of reaction. For LIBs, the one-step effective reaction considered is: The input parameters for the Gpyro model are listed in Table 3. The most attractive feature of Gypro is that it is developed specifically for pyrolysis modelling for combustible solids, therefore already embedding many interfaces with consideration of most situations that are encountered in pyrolysis and fire modelling. Users can easily set up the parameters and simplify the problem. The limitation is that if the problem exceeds equations provided, it is not easy for users to make their new definitions, as it requires alters of the source code directly.

Set Up for Benchmarking
This study conducts four case studies to compare and analyse the numerical performances with increasing complexity for (1) kinetics at the microscale (pure chemical reactions), (2) heat transfer at the mesoscale (pure heat transfer), (3) selfheating ignition of a single cell, and (4) a four-cell ensemble. Figure 1 is the schematic of the four case studies and Table 4 shows the details. Case (1) and case (2) consider a ideal condition while case (3) and case (4) are compared against the oven experiments for Sanyo LiCoO 2 (LCO) prismatic cells [28,29,39]. The detailed experimental setup is introduced in our previous work [28,29]. The dimensions of these cells are 34 9 10 9 50 mm, with a nominal capacity and voltage of 1.88 Ah and 3.7 V, respectively. In case 1, an adiabatic condition is considered for a small battery (such as a coin cell) with a mass of 1 g at an initial temperature of T 0 . In this condition, all Table 3 Heat transfer and kinetic parameters used in Gpyro simulations and simulations by COMSOL with the one-step reaction model for Sanyo LiCoO 2 (LCO) prismatic cells at 100% SOC. [29] Parameters Description Value Unit Figure 1. Diagram of 4 cases of benchmarking between COMSOL and Gpyro. Case 1 benchmarks the chemical kinetics by considering selfheating of 1 g of the sphere shape LIB with an initial temperature of T 0 placed in an adiabatic condition. Case 2 considers the heat transfer between the bulk of solid sample with a constant internal heat source q 000 c and ambient at T a . Case 3 and case 4 benchmark the oven experiments of a single cell and a four-cell ensemble, respectively. the heat generated by internal chemical reactions heats up the battery cell. Therefore, the temperature information is directly related to the chemical kinetics, and can be used to benchmark chemical kinetics, as shown in Eq. (22): With such a small size of the battery cell, the internal temperature gradient can be ignored, and the cell can be assumed to be a lumped system, which eliminates the effects of heat transfer and focuses on chemistry. While Gpyro considers a onestep global reaction, COMSOL consider multi-step reactions. Therefore, in COM-SOL, we try both the one-step global reaction scheme (same as Gpyro) and fourstep reaction scheme to benchmark with Gpyro. One issue for the four-step reaction scheme is that the kinetic parameters are based on an old generation of E-One/Moli Energy ICR18650 (18 mm diameter, 65 mm length) 1.65 Ah LCO/graphite cells at 100% state of charge (SOC) by the work conducted by Hatchard et al. [20] in 2001. The battery techniques in terms of energy density and material safety have improved dramatically since then, and the parameters may no longer suit the current LCO batteries. Also, it is difficult for the industry to obtain all the detailed kinetic parameters for a specific battery of interest, as it will require fundamental experimental analysis of every separate component making up the battery.
In this study, the energy ratio b is used to consider the effects resulting from the different energy densities between batteries. The heat generated by the four-step reactions is assumed to proportionally change with the stored chemical energy (which is roughly proportional to the electrochemical capacity). Equation (23) is then modified as: b Table 4 The list of shape and dimensions, boundary conditions, and internal heat sources of four cases where, C is the electrochemical capacity of the targeted battery and, C ref is the electrochemical capacity of the reference battery, which is 1.65 Ah. This work tries to investigate if the adapted function based on the kinetics [20] developed 20 years ago could give a reasonable prediction of the self-heating ignition behaviour of the current LCO batteries. Case 2 considers a heat transfer scenario. A constant heat source is used in this case to eliminate the complex chemical reactions. The geometry considered in this case is the dimensions of the four-cell ensemble in the oven experiment, which is 34 9 40 9 50 mm. Gpyro solves equations in one dimension, and adopts a volumetric heat transfer coefficient, h vl [29] to consider the heat transfer from the other two dimensions. COMSOL supports simulations on a complex 3D geometry. Therefore, both 1D and 3D simulations are conducted for COMSOL to benchmark with Gpyro. The boundary conditions consider both convective and radiative heat transfer with the ambient.
Case 3 and case 4 simulate the oven experiments on a single cell and a four-cell ensemble, respectively. The two cases consider the coupling effects of both chemistry and heat transfer when predicting self-heating ignition of the batteries. Both a one-step reaction scheme and a four-step reaction scheme are simulated by COMSOL and compared with the results of Gpyro and the oven experiments [28,35]. Convective and radiative heat transfer with the surroundings is considered for the boundaries.
Since case 4 represents the most complex condition, we conducted mesh independence analysis on case 4 for both COMSOL and Gpyro to determine the suitable mesh size and time step. For 1-D simulations in the depth direction with a total length of 40 mm, both COMSOL and Gpyro have converged temperature profiles at 50 nodes, where further increasing node number would not increase accuracy. For the 3D simulations by COMSOL, the mesh with 560 elements already predicts almost the same results as the mesh with 1080 elements (less than 1% of relative errors). Therefore, the meshes with 50 nodes and 560 elements are selected for 1D and 3D simulations. For the time step, COMSOL has a smart algorithm to self-adapt the time step used in simulation, for example, BDF (backward differentiation formula), which will automatically increase the time step under the premise of the settled tolerance. The time step can be further constrained by setting a maximum step constraint, where the maximum time step could not exceed the settled value. In this study, simulations with a maximum time step of 6 s already converged, which has almost the same result with a maximum time step of 0.6 s. Gpyro usually uses a fixed time step pre-settle by the user, where if the simulation could not converge under the given time step it will automatically decrease the time size until the settled tolerance is satisfied. In this study, the time step used in Gpyro simulations is 0.01 s. Figure 2 shows the benchmarking results for case 1 considering microscale chemistry. The heat transfer and chemical kinetic parameters for Gpyro are listed in Table 3. The one-step reaction model considered in COMSOL uses the same parameters as Gpyro for comparison. The parameters for the four-step reaction model considered in COMSOL is the same as the parameters listed in Table 1, except for the thermophysical properties, q, c p , k, which uses the values listed in Table 3. The energy ratio, b, considers two values: b = 1, which represents the same chemical kinetics for the battery in Hatchard et al. [39] with 1.65 Ah, and b = 1.14, which is adapted to the battery with 1.88 Ah [28,29]. Two initial temperatures are considered: T 0 = 130°C, and T 0 = 140°C. Focusing on self-heating igniton phenomenon of LIB, we consider temperature below 200°C, that after this temperature, the battery involves a more complex combustion stage and is out of the scope of this research.

Microscale Chemistry, Case #1
As Fig. 2 shows, the predictions by the one-step model of COMSOL and Gpyro agree well for both ambient temperatures. At T 0 = 130°C, a small increase of temperature (less than 2.5°C) after 200 min is predicted for both onestep models in COMSOL and Gpyro. The internal reaction rates are limited, with low heat generation at this initial temperature. For T 0 = 140°C, both models predict the exponential increase of temperature after 75 min, while the temperature predicted by COMSOL one-step model is a little faster (3 min) ahead than Gpyro. Those small differences may result from the different numerical methods (FEM for COMSOL and FDM for Gpyro) and different criteria for numerical convergence by the two tools. Overall, COMSOL and Gpyro agree well with the one-step modelling of chemical dynamics.
However, the four-step model by COMSOL considering two b values predicts a fast temperature increase at T 0 = 130°C. For b = 1, which represents the same kinetic parameters as Table 1, the temperature is predicted to grow exponentially at t = 37 min, which is even earlier than the result from the one-step models at T 0 = 140°C. The temperature for b = 1.14, which is modified by considering the variation of energy density, are predicted to grow faster, around 10 min ahead of that for b = 1. These fast temperature increases by the four-step model are caused by the SEI decomposition reaction, which is predicted to take place at around 100°C and generates heat. The heat generated by the SEI decomposition reaction will raise the temperature quickly in the adiabatic condition and initiate other reactions, resulting in a rapid temperature increase in a short period. Since the four-step reactions model already predicts exponential temperature growth at T 0 = 130°C, the scenario for T 0 = 140°C is no longer needed to be compared. In regards to the simulation time required, it took COMSOL around 8 min to run a one-step simulation and 10 min for a four-step simulation, while Gpyro took around 20 min for a one-step simulation. COMSOL has higher computational efficiency.  boundary conditions as Gpyro 1-D model. A 3-D heat transfer is also considered for COMSOL. A constant heat source (q 000 c ¼ 10 kW m -3 ) is assumed for the solid bulk to eliminate chemical dynamics. Convective and radiative heat transfer is considered for the solid bulk and the ambient at T a = 140°C. The thermophysical properties and heat transfer coefficients are the same as Table 3.

Mesoscale Heat Transfer, Case #2
As Fig. 3 shows, all three curves plateau at around 150°C. The temperature prediction by the COMSOL 1-D model has the same trend as the results by Gpyro, but is a little bit (around 3°C) higher. The 3-D model by COMSOL gives a higher prediction of temperature in the first 50 min. This difference may result from the consideration of 3D heat transfer. Gpyro and the COMSOL 1-D model use a volumetric heat transfer coefficient h vl as a source term to consider the heat transfer from other dimensions. This treatment assumes the same temperature distribution for each cross-section perpendicular to the dimension analysed. While the COMSOL 3-D model directly solves equations in three dimensions, which can consider the temperature gradience in the cross-sections. Overall, the final equilibrium temperatures by all three models are in agreement, which shows the same capability of modelling mesoscale heat transfer for COMSOL and Gpyro. For the computational efficiency, COMSOL has a better performance and took around 10 min (5 min faster than Gpyro). Figure 4 shows the benchmarking results for case 3, which are compared against the oven experiments on a single prismatic LCO cell. The COMSOL one-step model uses the same governing equations, boundary condition, kinetic and heat transfer parameters as the model developed by Gpyro 1-D model. The COMSOL Figure 4. The benchmarking results for case 3 (a single prismatic LCO cell) using Gpyro, COMSOL 1-step, and COMSOL 4-step models with two energy ratios b = 1 and 1.14 against oven experiment. a represents subcritical conditions and b represents supercritical conditions. The results for the experiments, Gpyro, COMSOL 1-step, COMSOL 4-step (b = 1), and COMSOL 4-step (b = 1) are presented by curves in red, blue, black, orange, and purple respectively (Colour figure online). four-step model considers two energy ratios b = 1 and 1.14 for 3D heat transfer. The corresponding kinetic parameters for the four reactions are the same as Table 1, while the heat transfer parameters are the same as Table 3. While selfheating problems especially focus on the transition from a safe state (subcritical condition) to an unsafe state (supercritical condition), we simulate the LIB under different ambient temperatures until the transition takes place. Figure 4a shows the subcritical conditions and Fig. 4b shows the supercritical conditions. The oven experiment showed that the temperature of the single-cell stabilized at T sub = 147°C (subcritical condition) and undertook thermal runaway at T sup = 149°C (supercritical condition).

Single-Cell Comparison, Case #3
Comparing the one-step models, Gpyro predicts the same subcritical and supercritical conditions as the experiment, while the COMSOL one-step model predicts T sub = 146°C for the subcritical condition and T sup = 147°C for the supercritical condition. The temperature profiles for the subcritical condition by the two models agree well with the experiment, with an error below 3°C. While for the supercritical condition, both models underestimate the time needed for the cell to develop into thermal runaway. Using temperature equals 200°C as a criterion, it took around 150 min for the experiment to reach that temperature, while Gpyro predicts around 50 min and the COMSOL one-step model predicts around 100 min at T a = 147°C. These results show that the reactivity of the one-step global reaction is a little overestimated and results in an earlier prediction of thermal runway. However, both models have a good prediction of the transition from the subcritical condition to the supercritical condition, which is the key parameter to assess self-heating ignition. Therefore, both COMSOL and Gpyro are capable of predicting self-heating ignition of LIBs with the one-step global reaction assumption.
For the four-step model, the four-step model was developed based on the work by Hatchard et al. [39] for a cylindrical LCO cell, which differs in its energy density from the prismatic cell tested in the experiments. In order to consider its effects, we introduced an energy ratio b. For b = 1, it represents exactly the same kinetic parameters as Table 1, the battery cell is predicted to stabilize at T sub = 153°C and thermal runaway at T sup = 158°C. While for b = 1.14, which adjusts the heat generation by considering the energy difference between the cylindrical cell in Hatchard et al. [39] and prismatic cell in this work, the predicted subcritical and supercritical conditions are T sub = 150°C and T sup = 155°C, respectively. The 3D heat transfer model predicts a faster temperature increase for the first 30 min, and the early-stage reactions such as SEI decomposition further accelerate this temperature increase rate. These results indicate that the prismatic cell analysed in this study may have a large difference in early-stage reactions compared with the cylindrical cell analysed in Hatchard et al. [39]. This may also result from the improvement of the safety of LIB techniques. For the time to reach 200°C, it took around 40 min for the scenario of b = 1 and 60 min for the scenario of b = 1.14, which are all lower than the prediction by one-step models. This is because the four-step models all predict a higher supercritical condition (T sup = 158°C and 155°C for b = 1 and 1.14), which is 5-10°C higher than the real experiment, and therefore a much faster reaction rate. For the computational time, it took COMSOL around 15 min for the 1-D simulation and 30 min for the 3D simulation. While Gpyro took around 30 min for the 1D simulation. Figure 5 is the benchmarking results for case 4, which compare against the oven experiments on a four-cell ensemble. The numerical set-up and parameters for the COMSOL one-step model, COMSOL four-step model, and Gpyro are the same as Sect. 3.3, except the geometry considered is a four-cell ensemble. This case with a larger size is supposed to demonstrate the effects of heat transfer.

Four-Cell Comparison, Case #4
For the one-step models, both COMSOL one-step model and Gpyro predict the same subcritical (T sub = 143°C) and supercritical (T sup = 138°C) conditions as the results obtained by the oven experiment. Again, the temperature profiles predicted by both COMSOL one-step model and Gpyro agree well for both subcritical and supercritical conditions, except the temperature predicted by COMSOL grows a little bit faster than Gpyro for the supercritical condition. The temperature measured in the experiment grows faster than the prediction by the two models in the initial 100 min for both subcritical and supercritical conditions, which may result from the variation of geometry and heat transfer. Comparing times taken to reach 200°C, it took around 230 min for the experiment to reach that temperature, while Gpyro predicts around 200 min and the COMSOL one-step model predicts around 160 min. Again, the predicted time by both models is lower than those measured in the experiments. Nevertheless, the accurate predictions on the transition from the subcritical to the supercritical condition demonstrate that Figure 5. The benchmarking results for case 4 (an ensemble of 4 prismatic LCO cells) using Gpyro, COMSOL 1-step, and COMSOL 4-step models with two energy ratios b = 1 and 1.14 against oven experiment. a represents subcritical conditions and b represents supercritical conditions. The results for the experiments, Gpyro, COMSOL 1-step, COMSOL 4-step (b = 1), and COMSOL 4-step (b = 1) are presented by curves in red, blue, black, orange, and purple respectively (Colour figure online).
both COMSOL and Gpyro are capable to predict self-heating ignition of LIBs with a one-step global reaction.
The COMSOL four-step model still predicts a higher temperature for both subcritical (T sub = 145°C for b = 1, and T sub = 142°C for b = 1.14) and supercritical (T sup = 150°C for b = 1, and T sup = 147°C for b = 1.14) conditions, as shown in Fig. 5. The temperature at the initial stage is still overestimated (maximum 25°C) compared to the experiments. For the time to reach 200°C, it took around 80 min for the scenario of b = 1 and 100 min for the scenario of b = 1.14. These large differences show that the LIB contains more complex reactions, which significantly affect the internal temperature change, especially for the supercritical conditions. The simplified one-step and four-step models are insufficient to predict the time to thermal runaway. For the computational time, it took COMSOL around 18 min for the 1-D simulation and 35 min for the 3D simulation. While Gpyro took around 40 min for the 1D simulation. Table 5 lists experimental results and predictions of the subcritical, supercritical, and calculated critical conditions for case 3 and case 4 by the different models discussed. Overall, Gpyro and the COMSOL one-step model have almost the same prediction as the experimental results for both single cell and the four-cell ensemble. The critical ambient temperature to trigger self-heating ignition, T a,cr , predicted by the COMSOL four-step model is around 7°C (b = 1) to 4°C (b = 1.14) higher than the experimental results for both single cell and the fourcell ensemble. This difference is acceptable considering the large differences in geometry, size, and improved techniques between the old generation of cylindrical cells from Hatchard et al. [39] and the prismatic cell in this study. Moreover, the decline of T a,cr for the single cell to the four-cell ensemble is accurately predicted for all models, that the experiments showed a 7.5°C drop and numerical models predict a 6 to 8°C decline. These results demonstrate that both COMSOL and Gpyro with effective one-step reaction schemes and four-step reaction schemes are capable of predicting the self-heating ignition of LIBs. Table 5 The subcritical and supercritical conditions found by the oven experiments and numerical modelling by Gpyro, COMSOL 1-step model, COMSOL 4-step model with b = 1 and 1.14 for the single-cell and the 4-cell ensemble Overall, both Gpyro and COMSOL accurately predict the transition from a subcritical condition to a supercritical condition, which is the key evaluation parameter to detect self-heating ignition of LIBs. However, for the detailed modelling of chemistry and heat transfer, both platforms have their strengths and weaknesses. Gpyro embeds a mature set-up for self-heating problems, which supports fast setup and simulations. Based on one-step reaction simplification, the required chemical parameters are much easier to obtain by experiments, which is easier for engineering and industry applications. However, it only supports simulations in which the reactants are in regular simple shapes. For complex geometry, the internal complex heat transfer may result in large errors. COMSOL has a customer-friendly interface and support solutions for complex multi-disciplines with complex shapes, which could study the fundamental mechanisms underneath and make reasonable predictions. However, it requires the user to have mature knowledge of the problem and be familiar with the software and take more time to setup. It also requires more input parameters, which might need significant tests to obtain for engineering applications.

Conclusions
This study presents a benchmarking analysis for COMSOL Multiphysics and Gpyro in predicting self-heating ignition of lithium-ion batteries. Four case studies with different levels of complexity are considered: (1) just kinetics at the microscale, (2) just heat transfer at the mesoscale, (3) self-heating behaviour at the mesoscale for coupled chemical reactions and heat transfer of a single cell, and (4) a four-cell ensemble for multiphysics at a larger scale. COMSOL one-step model predicts almost the same (within 1% error) results as Gpyro, while the COMSOL four-step model predicts an earlier (temperature sharp increase at initial temperature 10°C lower) temperature increase at the initial stage due to the early-stage reactions. For the heat transfer at mesoscale, the 3-D heat transfer model by COMSOL predicts a faster (10 min ahead) temperature increase at the initial stage, compared to the results by the 1-D model by Gpyro and COMSOL, which agrees well. For the overall prediction of the self-heating behaviour of LIB at mesoscale, both COMSOL and Gpyro give good predictions of the transition from a subcritical to a supercritical condition, which are also validated against the oven experiments for a single cell and a four-cell ensemble.
In general, COMSOL and Gpyro have large differences in their physical modelling, mathematical modelling, and numerical modelling. In this benchmarking study, we consider the effects of energy difference for old (COMSOL) and new (Gpyro) LiCoO 2 battery cells, one-step (Gpryo) and four-step (COMSOL) reaction kinetics for chemistry, 3D modelling (COMSOL) and Pseudo 3D modelling (Gpyro) for heat transfer, and FEM (COMSOL) and FDM (Gpyro) for simulation performance. For physical modelling, Gpyro assumes a 1-D simplification for geometry while COMSOL supports direct 3D simulation on complex geometry. For mathematical modelling, Gpyro considers more complex physics which solve heat and mass conservation equations for both gas and solid phases but consider an overall global one-step reaction, while COMSOL adopts simplifications to solve conservation equations for solid-phase but support simulations with multistep reactions. For numerical modelling, Gpyro and COMSOL use different numerical discretization methods, numerical calculations, algorithms, and computational efficiency. However, this study proves that both tools can accurately predict the critical ambient temperature to trigger thermal runaway, which is essential for predicting self-heating ignition of LIBs. This paper validates their use to study the safety of LIBs. Although this study focuses on self-heating, the validation also indicates the tools could be used for other LIB ignition events.