Battery Thermal Runaway User Guide

The AcuSolve/SimLab battery solution for thermal runaway provides a virtual experimental platform to test pack designs for thermal safety, including the assessment of propagation characteristics, heat shield effectiveness, and general understanding of thermal abuse tolerance. Three models are available to simulate thermal runaway:

NREL abuse model
Arc reaction heat model
Heat rate model

All models provide a heat-source ( $S_{T R}$ ) that is applied to the energy equation. Further details of their physics and governing equations can be found in the Thermal Runaway Theory Manual.

NREL Abuse Model Formulation

The NREL model builds on previous models in literature ([1], [2], and [3]). The model solves four or five equations representing different material decomposition reactions, which are exothermic in nature, under thermal abuse conditions. These reactions can lead to uncontrollable heating and subsequently thermal runaway. The choice between four or five equations depends on whether the decomposition of the electrolyte is included as a decomposition reaction. Each of these decomposition reaction rates is governed by an Arrhenius equation. The heat generated is added as a source term in the energy equation based on the component mass and reaction enthalpy. To determine the type of model, the type command is set to NREL_abuse_model and equation_type is set to either four_equation or five_equation. For example:

BATTERY_THERMAL_RUNAWAY_MODEL( "NREL_model" ) {
    …
    equation_type = four_equations
    type = NREL_abuse_model
}

Details of the governing equations can be found in the Thermal Runaway Theory Manual.

The thermo-physical and chemical kinetic parameters required for each Arrhenius reaction equation and associated heat generation include:

Mass of the cell components (anode, cathode, SEI (Solid Electrolyte Interphase), and electrolyte)
Frequency factor
Activation energy
Reaction enthalpy

These parameters can be obtained from literature, experiments, for example, parameter fitting of DSC data, or by selecting pre-defined models (LCO (Lithium Cobalt Oxide, LiCoO2), NMC (Nickel Manganese Cobalt) and NCA (Nickel Cobalt Aluminium)) in SimLab. An example of LCO parameters for the cathode decomposition would be:

BATTERY_THERMAL_RUNAWAY_MODEL( "thermosphysical_chemical_kinetic" ) {
    …
    type = NREL_abuse_model
    cathode_decomposition_frequency_factor = 6.67000E+11
    cathode_decomposition_activation_energy = 2.03000E-19
    cathode_decomposition_reaction_enthalpy = 314000.0
    …
}

The mass of the individual cell components can be provided as either a total mass (kg) or a volume specific mass (kg/m³). Details of this calculation are provided in the section Mass calculation (NREL model). To set the mass type the mass_type command is set to total for total mass or specific for volume specific mass. The masses of each of the components: anode_mass, cathode_mass, electrolyte_mass, and sei_mass are then set in either kg (for total mass) or kg/m³ (for specific mass). For example:

BATTERY_THERMAL_RUNAWAY_MODEL( "total_mass_example" ) {
    …
    type = NREL_abuse_model
    mass_type = total
    anode_mass = 0.006
    cathode_mass = 0.012
    electrolyte_mass = 0.004
    sei_mass = 0.012    
    …
}

An additional sixth Arrhenius equation, based on the work by [3], can be included to model an internal short circuit event. The equation requires the following parameter definitions:

Internal short circuit temperature
Frequency factor
Activation energy
Reaction enthalpy

The internal short circuit temperature is the trigger temperature that enables the solution of the sixth equation, that is, prior to the trigger temperature the reaction is deactivated. To enable an internal short circuit the command internal_short_circuit is set to on and the internal_short_circuit_temperature is provided. As an example:

BATTERY_THERMAL_RUNAWAY_MODEL( "internal_short_circuit ) {
    …
    type = NREL_abuse_model
    internal_short_circuit = on
    internal_short_circuit_temperature = 430.0    
    …
}

A typical AcuSolve input for the NREL abuse model based on LCO chemistry for a five-equation model would be defined by

BATTERY_THERMAL_RUNAWAY_MODEL( "five_equation " ) {
    equation_type = five_equations
    formulation_type = lumped
    ecm_model = no
    type = NREL_abuse_model
    internal_short_circuit = off
    internal_short_circuit_temperature = 430.0
    mass_type = total
    anode_mass = 0.006
    cathode_mass = 0.012
    electrolyte_mass = 0.004
    sei_mass = 0.012
    anode_decomposition_frequency_factor = 2.50000E+13
    anode_decomposition_activation_energy = 2.24000E-19
    anode_decomposition_reaction_enthalpy = 1.714000E+6
    cathode_decomposition_frequency_factor = 6.67000E+11
    cathode_decomposition_activation_energy = 2.03000E-19
    cathode_decomposition_reaction_enthalpy = 314000.0
    sei_decomposition_frequency_factor = 1.67000E+15
    sei_decomposition_activation_energy = 2.24000E-19
    sei_decomposition_reaction_enthalpy = 257000.0
    electrolyte_decomposition_frequency_factor = 5.14000E+25
    electrolyte_decomposition_activation_energy = 4.54000E-19
    electrolyte_decomposition_reaction_enthalpy = 155000.0
    initial_anode_lithium_fraction = 0.75
    initial_sei_lithium_fraction = 0.15
    initial_sei_thickness = 0.033
    initial_degree_of_cathode_conversion = 0.04
    initial_electrolyte_fraction = 1.0
    sub_iteration_minimum_time_step_size = 0.0001
    sub_iteration_time_step_size_tolerance = 0.0001
}

The solution of the ODE for the Arrhenius reactions are stiff and require time-step adaptation for an accurate and efficient solution. The exact details of the adaptive time-stepping scheme are given in Adaptive time stepping (NREL and multistage model).

ARC Reaction Model

The ARC reaction model is a staged Arrhenius-based kinetic model that can be directly fitted to Accelerating Rate Calorimetry (ARC) data. An example of typical ARC data is shown in Figure 1.

Figure 1. Heat-Rate Versus Temperature for a Typical ARC Test

To enable the arc reaction, model the type is set to arc_reaction_model.

BATTERY_THERMAL_RUNAWAY_MODEL( "ARC reaction model" ) {
    …
    type = arc_reaction_model
    …
}

This model requires the ARC data to be divided into stages or bands that loosely describe physical processes that occur during thermal runaway.

Note: Since it is a data fitting approach some of the chemical processes can be lumped together. More details are given later.

The number of stages is determined by you during the data fitting phase. Figure 2 illustrates an example of dividing the ARC data into four stages based on temperature, where Ts represents the start temperature and T₁, T₂, and T₃ represent the end temperatures of each stage one, two, and three, respectively. The actual stages are defined as: Stage 1: T₁-T_s; Stage 2: T₂-T₁; Stage 3: T₃-T₂; Stage 4: T_max-T₃. T_max is the maximum temperature recorded from the ARC test as is determined automatically from the data.

Figure 2. Four Stage Fit Showing T_s, T₁, T₂, T₃ and the Division of Stages

Although an ARC test does not individually identify the decomposition of the reacting components, it can still be represented by Arrhenius kinetics with physical meaning depending on the number of stages, for example, two stages may group together SEI, anode decomposition into the first reaction kinetic and cathode conversion and electrolyte decomposition into the second reaction kinetic. Using optimization methods, the unknown parameters can be determined for the following system of equations:

\frac{{dα}_{i}}{dt} = α_{i}^{n} {(1 - α_{i})}^{m} \cdot A_{(a, i)} \cdot e x p (- E_{(a, i)} / (k_{b} T))

m c_{p} \frac{d T}{d t} = \sum m_{i} c_{p, i} Δ T_{i} η_{i} \cdot \frac{{dα}_{i}}{dt}

The unknown parameters for each stage (

i

) are:

$A_{(a, i)}$ : The frequency factor (units: 1/s).
$E_{(a, i)}$ : The activation energy (units: J).
$n,$ and $m$ : The parameters that determine the reaction order and type, for example, auto-catalytic reaction when m>0.
$η_{i}$ : The reaction enthalpy adjustment factor.

The fitting is performed on a 0D model that ignores conduction. However, in the final simulation the right-hand side of the above equation is used as the heat-source in AcuSolve.

S_{T R} = \sum_{i} Q_{i} = \sum_{i} m_{i} c_{p, i} Δ T_{i} η_{i} \frac{d α_{i}}{d t}

The input for the multi-stage ARC reaction model is given as an array of stages and Arrhenius kinetic/material parameters.

As an example, a two-stage model would have the following structure.

Table 1.
Stages	Enthalpy ( $h_{i}$ )	Activation energy ( $E_{a, i}$ )	Frequency factor ( $A_{a, i}$ )	Reaction order ( $n$ )	Reaction order ( $m$ )	Initial condition ( $α_{i}$ )
1	176620.88	1.5E-19	-2.2551e6	1	0	1
2	1521147.329	2.9445e-19	2.31741e15	0	7.5	0

Note: Activation energy can be negative to account for a negative sign in the Arrhenius reaction equation.

The above in the .inp file would be given as follows:

BATTERY_THERMAL_RUNAWAY_MODEL( "two_stage_fit"){
    type = arc
    arc_type = multi_stage_fit
    multi_stage_model_parameters = { 176620.88, 1.5E-19, -2.2551E6, 1, 0, 1;
                                    1521147.329, 2.9445E-19, 2.31741E15, 0, 7.5, 0; }
}

Additionally, for a higher number of stages (three and four) to improve the accuracy of the final fit the heat produced by the final stage is not applied to the energy equation until above a specified temperature. The ODE solution is still solved throughout. The typical input, as an example for a four-stage fit, would be given by

BATTERY_THERMAL_RUNAWAY_MODEL( "four_stage_fit " ) {
    ecm_model = no
    type = arc_reaction_model
    formulation_type = lumped
    arc_type = multi_stage_fit
    parameter_fit_type = arrhenius_trigger
    multi_stage_model_parameters = Read( "ms.txt" )
    heat_rate_trigger_temperature = 500
}

heat_rate_trigger_temperature is the temperature at which the final stage heat is added to the energy equation. In this approach the reaction rate is still solved for the entire temperature range for the final stage. This is identified from the ARC data directly as the approximate temperature where the temperature jumps considerably. See temperature T₃ in Figure 2.

If a trigger temperature is included in the fitting process, it must be included in the simulation and vice versa.

The multi_stage_model_parameters can be read directly from a file, in which the file is formatted, as an example for ms.txt, as follows:

2.2244e+03 3.0693e-19 -1.2420e+20 1.0 0.0 1.0
2.1507e+03 3.8971e-19 -1.0e+25 1.0 0.0 1.0
…

Where each line of the file represents the parameters for a specific stage.

The ARC reaction model has three parameters related to the adaptive time-stepping scheme used to solve the governing thermal runaway equations. The exact details of the adaptive time-stepping scheme are given in Adaptive time stepping (NREL and multistage model).

Heat-Rate Model (Direct Reading of ARC Data)

A direct reading of thermal runaway ARC data is enabled via the following command:

BATTERY_THERMAL_RUNAWAY_MODEL( "Heat rate model" ) {
    …
    type = heat_rate
    …
}

In this approach a source term is calculated from the ARC data and applied directly in the energy equation given by:

S_{T R} = ρ_{c e l l} c_{p, c e l l} \frac{d T}{dt}

Where

ρ_{c e l l}

is the effective density of the cell,

c_{p, c e l l}

is the effective specific heat of the cell. These two terms come from the MATERIAL_MODEL defined for the cell.

d T / d t

is the heat rate data in K/s from the ARC test and is read from a file. A typical input for direct reading ARC data would be given as follows:

BATTERY_THERMAL_RUNAWAY_MODEL( "my_arc_model"){
    type = heat_rate
    formulation_type = lumped
    heat_rate_type = linear
    heat_rate_curve_fit_values = Read( "dTdt.txt" )
    heat_rate_curve_fit_variable = temperature
    heat_rate_max_temperature = 740
    heat_rate_trigger_temperature = 490
    heat_rate_onset_temperature = 400
}

In the above example, heat_rate_curve_fit_values points to a file: dTdt.txt. This file is a two-column array of temperature and heat-rate values, for example:

…
402.752 0.006483333
402.978 0.0075
403.165 0.006233333
403.366 0.006633333
403.531 0.005516667
403.692 0.00535
403.853 0.005433333
403.993 0.004616667
404.149 0.005183333
…

The direct-read heat-rate model defines specific temperatures to improve accuracy and enable or disable the model. These parameters are: heat_rate_max_temperature, heat_rate_trigger_temperature, and heat_rate_onset_temperature.

heat_rate_max_temperature: This parameter provides the maximum temperature recorded during the ARC test.
heat_rate_trigger_temperature: This parameter provides the temperature at or around the point when thermal runaway occurs, for example, T₃ from the ARC data in Figure 2, and allows the ODE to sub-step at a time-step smaller than the PDE time-step of the energy equation.
heat_rate_onset_temperature: This parameter is the onset temperature for self-heating and can be directly identified from the ARC data.

Additional Model Information

Mass calculation (NREL model)

In the NREL thermal runaway model, the mass of individual cell components can be expressed either as total mass or specific mass. Below is a brief description of the theory.

The mass of the cathode, anode, and electrolyte can be computed using the formula:

m_{T} = m_{s} V

, where

m_{T}

is the total mass,

m_{S}

is the specific mass, and

V

the volume of the jellyroll. The jellyroll volume can be the whole cell in the case of a homogenized representation of the battery or the actual jellyroll geometry. For instance, a homogenized 18650 cell may have a jellyroll volume of

1.66 \times 10^{- 5} m^{3}

, while a more precise estimation based on studies by [2] and [4] yield volumes of

1.05 \times 10^{- 5} m^{3}

, and

1.29 \times 10^{- 5} m^{3}

, respectively. Alternatively, if a detailed and accurate jellyroll geometry is available, it can be directly utilized.

Homogenized 18650 LCO cell: If the cell volume is $1.66 \times 10^{- 5} m^{3}$ and the specific masses of the cathode and anode of the cell are 722.89 $k g / m^{3}$ and 351.45 $k g / m^{3}$ , respectively, then the total mass of the cathode is 0.012 $k g$ , and 0.006 $k g$ for the anode.
Detailed Jellyroll 18650 LCO cell: If the volume of the jellyroll is $1.29 \times 10^{- 5} m^{3}$ and the specific masses of the cathode and anode are again 722.89 $k g / m^{3}$ and 351.45 $k g / m^{3}$ , respectively, then the total mass of the cathode is 0.0093 $k g$ and 0.0047 $k g$ for the anode.

Adaptive time stepping (NREL and multistage model)

Adaptive time stepping is necessary for the ODE-based thermal runaway models to increase computational efficiency. In the adaptive time-stepping algorithm, the ODE is sub stepped based on how the solution evolves over time. The algorithm updates the ODE sub stepping size once per PDE solution update. To determine the ODE time-step size, an estimate of how rapidly the solution is evolving is calculated using the formula:

e_{n} = \frac{y_{n} - y_{n - 1}}{1 + min (y_{n}, y_{n - 1})}

where $y_{n}$ is solution at the current time step, $y_{n - 1}$ is solution at the previous time step, $e_{n}$ is calculated for each ODE solution variable for all cells, with the maximum value being selected among them. Based on control theory, it can be shown that the updated time step can be given by an integral controller:

{Δt}_{I} = \frac{T o l}{e_{n}} \cdot Δ t_{n - 1}

Here,

T o l

is a user-defined tolerance, which determines how sensitive the time step should be when solution is changing,

Δ t_{n - 1}

is the time step size at the previous time step. Once

{Δt}_{I}

is determined, it is bound between two user-defined limits:

Δ t_{min}

and

Δ t_{max}

Note:

Δ t_{max}

cannot exceed

Δ t_{P D E}

, which is the PDE time step size, that is:

{Δt}_{n} = min (max (Δ t_{min}, {Δt}_{I}), Δ t_{max})

The adaptive time stepping commands for

T o l

and

Δ t_{min}

are provided in the BATTERY_THERMAL_RUNAWAY_MODEL section of the .inp file and are given by the commands sub_iteration_minimum_time_step_size and sub_iteration_time_step_size_tolerance, respectively. An example input would be:

BATTERY_THERMAL_RUNAWAY_MODEL( "my_arc_model"){
    …
    sub_iteration_minimum_time_step_size = 0.0001
    sub_iteration_time_step_size_tolerance = 0.0001
}

References

[1] T. D. Hatchard, D. D. MacNeil, A. Basu and J. R. Dahn, "Thermal model of cylindrical and prismatic lithium-ion cells," Journal of The Electrochemical Society, p. A755, 2001.

[2] K. Gi-Heon, A. Pesaran and R. Spotnitz, "A three-dimensional thermal abuse model for lithium-ion cells," Journal of Power Sources, pp. 476-489, 2007.

[3] P. T. Coman, E. C. Darcy, C. T. Veje and R. E. White, "Modeling Li-ion cell thermal runaway triggered by an internal short circuit device using an efficiency factor and Arrhenius formulations, "Journal of The Electrochemical", p. A587, 2017.

[4] P. J. Bugryniec, J. N. Davidson and S. F. Brown, "Computational modeling of thermal runaway propagation potential in lithium iron phosphate battery packs," Energy Reports, pp. 189-197, 2020.