An Electrochemical-Thermal Model for Lithium-Ion Battery Packs during Driving of Battery Electric Vehicles

In this study, Li-ion battery packs that are large enough to be installed in actual BEVs were targeted. Such a battery pack does not actually exist and is only a design assumed for this model. An overview of the battery pack is shown in Fig. 1, and its performance is presented in Table I. The battery pack consists of 288 laminate cells. Although the cooling systems in several BEVs are based on liquid cooling ^4,5 or cooling with phase change materials, in this model there is only air cooling to the outside. This is because the modeling does not focus on the impact of the cooling system, and this design could make the discussion of the results easier to understand. In addition, there are BEVs that do not have a cooling system, as in this case, and we do not believe that this design will cause any deviation from reality. The laminate cells are numbered No. 1–18, starting with the one at the bottom. The laminate cell is assumed to be a common NMC/graphite (Li_y (Ni_1/3 Co_1/3 Mn_1/3)O₂ /Li_x C₆) system.

Zoom In Zoom Out Reset image size

Most system simulations adopt models based on electrical ECMs as battery models because they have low calculation costs. Meanwhile, models based on electrochemical phenomena, such as the Newman pseudo two-dimensional model, ^16,22 have the advantage of being faithful to the electrochemical mechanism, which makes the understanding of the results obvious. In particular, they are used for cyclic degradation and internal short circuits. In the proposed model, the single particle (SP) model ²³ is adopted for each cell. This model can reproduce electrochemical phenomena with a relatively low computational cost.

During the charge/discharge process in a Li-ion battery, Li ions are transported from the positive/negative electrode to the opposite electrode through interfacial reactions on the active materials and diffusion and migration in the electrolyte. Moreover, electrochemical reactions occur at the interface of the active material in the positive/negative electrode. In the SP model, the processes are formulated as a diffusion equation and electrochemical equations, as presented in Table II. Note that these equations are adopted from a previous report ²³ and are not developed in this study. The model assumes that the effects of diffusion and migration in the electrolyte are negligible.

Table II. Governing equations of the single particle model. ²³ Subscript j indicates a positive or negative electrode. Note that these equations are adopted from a previous report ²³ and are not developed in this study.

Items	Equations
Active material	$\begin{array}{c}\displaystyle \frac{\partial {c}_{s,j}}{\partial t}-\displaystyle \frac{{D}_{s,j}}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial {c}_{s,j}}{\partial r}\right)=0\end{array}$  [1]
Active material/electrolyte interface	$\begin{array}{c}\displaystyle \frac{\partial {c}_{s,j}}{\partial r}=\displaystyle \frac{{j}_{\mathrm{Li},j}}{F}\end{array}$  [2]
	$\begin{array}{c}{j}_{\mathrm{Li},j}={i}_{0,j}\left\{\exp \left(\displaystyle \frac{{\rm{0.5}}F{\eta }_{j}}{{\rm{RT}}}\right)-\exp \left(\displaystyle \frac{-{\rm{0.5}}F{\eta }_{j}}{{\rm{RT}}}\right)\right\}\end{array}$  [3]
	$\begin{array}{c}{i}_{0,j}=F{k}_{j}{c}_{s,j}^{{\rm{0.5}}}{\left({c}_{s,\max ,j}-{c}_{s,j}\right)}^{{\rm{0.5}}}{c}_{l}^{{\rm{0.5}}}\end{array}$  [4]
	$\begin{array}{c}{\eta }_{j}={\phi }_{s,j}-{E}_{\mathrm{eq},j}\left({c}_{s,j}\right)\end{array}$  [5]

From the electric current conservation law in each battery cell, the relationship between the current applied to the battery cell ${I}_{{\rm{app}}}$ and the current density on the active material interface ${j}_{{\rm{Li}}}$ is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{j}_{\mathrm{Li},j}=\pm {I}_{{\rm{app}}}/{S}_{\mathrm{tot},j},\end{array}\end{eqnarray} \tag{ 6 }$$

where ${S}_{{\rm{tot}}}$ is the total interface area of the active materials, and subscript j indicates the positive or negative electrode. ${j}_{{\rm{Li}}}$ has opposite signs for the positive and negative electrodes. ${I}_{{\rm{app}}}$ has a positive charge value and a negative value on discharge. ${S}_{{\rm{tot}}}$ is calculated using the following equations:

$$\begin{eqnarray}&&\begin{array}{c}{S}_{\mathrm{tot},j}={V}_{s,j}{a}_{v,j}={S}_{{\rm{sheet}}}{t}_{s,j}{\varepsilon }_{j}\left(3/{r}_{p,j}\right),\end{array}\end{eqnarray} \tag{ 7 }$$

where ${V}_{s}$ is the volume of the active material in each cell, ${a}_{v}$ is the specific area of the active material, ${S}_{{\rm{sheet}}}$ is the area of the electrode sheet, ${t}_{s}$ is the thickness of the electrode sheet, $\varepsilon$ is the volume fraction of the active material, and ${r}_{p}$ is the radius of the active material particle.

The input parameters of the model are presented in Table III. These values are the same as those of our prototype cell reported in a previous study, ¹⁷ although the area of the sheet is adjusted to match the capacity of 57.5 Ah. The parameters were also validated with the actual discharge curve of the alternating current impedance in a state of charge (SOC) of 100 %, ¹⁷ and it was confirmed that these values were similar to those reported in other studies. ²³ In detail, the capacity is defined as the amount of electrons applied to the battery in order for the open-circuit voltage (OCV) to be between 3.0V and 4.2V, and this value is 57.5Ah per cell. Since OCV assumes no dependence on temperature in this study, the capacity does not change with temperature. In addition, the improvement of Li diffusivity at high temperatures is taken into account shown in Table III, and the decrease in internal resistance is taken into account.

Table III. Input parameters for the single particle model. ¹⁷ T indicates the temperature.

Item	Parameter	Unit	Value
Positive electrode	Thickness ${t}_{s}$	μm	42
	Volume fraction of the active material $\varepsilon$	—	0.29
	Particle radius ${r}_{p}$	μm	4.85
	Diffusion coefficient ${D}_{s}$	m2s−1	2.2 × 10−12 exp(−9.2 × 102/T)
	Rate coefficient $k$	m2.5mol−0.5 s−1	1.8 × 10−8 exp(−1.1 × 103/T)
Negative electrode	Thickness ${t}_{s}$	μm	55
	Volume fraction of the active material $\varepsilon$	—	0.39
	Particle diameter ${r}_{p}$	μm	13
	Diffusion coefficient ${D}_{s}$	m2s−1	9.9 × 10−12 exp(−1.6 × 103/T)
	Rate coefficient $k$	m2.5mol−0.5 s−1	1.5 × 10−8 exp(−1.6 × 103/T)
Common	Area of the electrode sheet ${S}_{{\rm{sheet}}}$	m2	3.55

The assumptions made for the model are as follows.

• (1)  
  The Li concentration in the active material in the cell is described as the concentration of two representative particles.
• (2)  
  The impact of the electrolyte and current-collecting foil is negligible.

The capacity of a Li-ion battery decreases and its internal resistance increases with repeated charging and discharging, which is known as cyclic degradation. In each cell, the proposed model considers the impact of the solid electrolyte interface (SEI) layer growth in negative electrodes, phase transition in positive electrodes, and crack propagation. These models were proposed in a previous study. ¹⁷ Table IV lists the governing equations. Note that these equations are adopted from a previous report ¹⁷ and are not developed in this study.

Table IV. Governing equations of the cycle degradation model. ¹⁷ Note that these equations are adopted from a previous report ¹⁷ and are not developed in this study.

Items	Equations
SEI layer growth in negative electrodes	$\begin{array}{c}{j}_{s}=\displaystyle \frac{{\rm{0.5}}{i}_{0,s}F{\eta }_{s}}{{\rm{RT}}}\end{array}$  [8]
	$\begin{array}{c}\displaystyle \frac{\partial {\delta }_{{\rm{sei}}}}{\partial t}=\displaystyle \frac{{M}_{{\rm{sei}}}}{F{\rho }_{{\rm{sei}}}}{j}_{s}\end{array}$  [9]
	$\begin{array}{c}{\eta }_{s}={\eta }_{s,0}-\displaystyle \frac{{\delta }_{{\rm{sei}}}}{{\kappa }_{{\rm{sei}}}}{j}_{s}\end{array}$  [10]
	$\begin{array}{c}{Q}_{{\rm{sei}}}=\displaystyle \frac{F{\rho }_{{\rm{sei}}}}{{M}_{{\rm{sei}}}}{\delta }_{{\rm{sei}}}{S}_{{\rm{tot}}}\end{array}$  [11]
Phase transition in positive electrodes	$\begin{array}{c}\displaystyle \frac{\partial {\theta }_{p}}{\partial t}={k}_{p}\left(1-{\theta }_{p}\right)\end{array}$  [12]
	$\begin{array}{c}{\rho }_{p}={\rho }_{{\rm{cubic}}}{\theta }_{p}\end{array}$  [13]
Crack propagating	$\begin{array}{c}{R}_{c}={R}_{c,0}\displaystyle \frac{{a}_{0}}{{a}_{0}-{a}_{c}}\end{array}$  [14]
	$\begin{array}{c}{a}_{c}={C}_{c}{Q}_{{\rm{tot}}}\end{array}$  [15]

As shown in Equation 15 in Table IV, in this model the crack propagation is assumed to be dependent on the total transported electric charge amount ${Q}_{{\rm{tot}}}.$ It is calculated by integrating the applied current to the battery cell. Therefore, in this model it is assumed that cracks in the cell are always generated by the application of current, and note that crack propagation is not predicted by stress evaluation. Furthermore, the resistance due to cracks is not estimated from the size of a single crack, but represents the frequency of occurrence of multiple cracks that resist the current flowing in the active material. It is assumed that there are no cracks in the initial state, and that if the current is sufficiently advanced, the resistance will be extremely high due to the numerous cracks. However, because of Equation 14 in Table IV, it is necessary to provide initial resistance due to cracking, and the value should be small enough to have no effect.

Owing to the SEI layer growth, the capacity loss of the negative electrode is taken into account by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{c}_{s,n}^{* }={c}_{s,n}-\displaystyle \frac{{Q}_{{\rm{sei}}}}{F{V}_{s,n}},\end{array}\end{eqnarray} \tag{ 16 }$$

where ${c}_{s}^{* }$ is the Li-ion concentration modified for the capacity loss, ${c}_{s}$ is the Li concentration, and ${Q}_{{\rm{sei}}}$ is the capacity loss. In the SP model, ${c}_{s,n}^{* }$ was used instead of ${c}_{s,n}.$ The variation in the overpotential of the main reactions due to SEI is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{\eta }_{n}={\phi }_{s,n}-{E}_{\mathrm{eq},n}\left({c}_{s,n}^{* }\right)+\displaystyle \frac{{\delta }_{{\rm{sei}}}}{{\kappa }_{{\rm{sei}}}}{j}_{\mathrm{Li},n},\end{array}\end{eqnarray} \tag{ 17 }$$

where $\eta$ is the overpotential, ${E}_{{\rm{eq}}}$ is the open-circuit potential, ${\delta }_{{\rm{sei}}}$ is the thickness of the SEI, and ${\kappa }_{{\rm{sei}}}$ is the electric conductivity of the SEI layer. The subscript n indicates the negative electrode.

The variation in overpotential due to the phase transition in positive electrodes is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{\eta }_{p}={\phi }_{s,p}-{E}_{\mathrm{eq},p}\left({c}_{s,p}\right)-{\rho }_{p}{j}_{\mathrm{Li},p},\end{array}\end{eqnarray} \tag{ 18 }$$

where ${\rho }_{p}$ is the resistance per unit cross-sectional area owing to phase transition. The subscript p indicates a positive electrode.

The voltage of each cell ${V}_{{\rm{cell}}}$ is calculated by the following equations:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{V}_{{\rm{cell}}}={\phi }_{s,p}-{\phi }_{s,n}+{R}_{c}{I}_{{\rm{app}}}=\left({E}_{\mathrm{eq},p}+\,{\eta }_{p}+{\rho }_{p}{j}_{\mathrm{Li},p}\right)\\ \,\,-\left({E}_{\mathrm{eq},n}+\,{\eta }_{n}-\displaystyle \frac{{\delta }_{{\rm{sei}}}}{{\kappa }_{{\rm{sei}}}}{j}_{\mathrm{Li},n}\right)+{R}_{c}{I}_{{\rm{app}}},\end{array}\end{array}\end{eqnarray} \tag{ 19 }$$

where ${\phi }_{s}$ is the difference in the electrical potential between the active materials and the electrolyte, and ${R}_{c}$ is the resistance due to crack propagation. In this model, crack propagation considers the impact of the cell as a whole cell, rather than the impact of the positive and negative electrodes individually. ${I}_{{\rm{app}}}$ has a positive charge value and a negative value on discharge. The ${I}_{{\rm{app}}}$ in each cell is adjusted so that the cell voltage ${V}_{{\rm{cell}}}$ is the same among the three cells connected in parallel.

Heat generation ${Q}_{r}$ is calculated by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{Q}_{r}={\eta }_{p}{j}_{\mathrm{Li},p}{S}_{\mathrm{tot},p}+{\eta }_{n}{j}_{\mathrm{Li},n}{S}_{\mathrm{tot},n}+{R}_{c}{I}_{{\rm{app}}}^{2}.\end{array}\end{eqnarray} \tag{ 20 }$$

The input parameters of the model are listed in Table V. These values are the same as those of our prototype cell reported in a previous study ¹⁷ and were determined using cycle degradation test data at 25 °C and 70 °C.

Table V. Input parameters for the cycle degradation model. ¹⁷ T indicates the temperature.

Item	Parameter	Unit	Value
SEI growth in negative electrodes	Initial overpotential ${\eta }_{s,0}$	V	0.3
	Exchange current density ${i}_{0,s}$	A m−2	5.4 × 1011 exp(−7.6 × 103/T)
	Molar weight of SEI ${M}_{{\rm{sei}}}$	kg mol−1	0.01
	Density of SEI ${\rho }_{{\rm{sei}}}$	kg m−3	1000
	Conductivity of SEI ${\kappa }_{{\rm{sei}}}$	S m−1	3.41
Phase transition in positive electrodes	Cover rate ${k}_{p}$	s−1	1.8 × 1010 exp(−1.0 × 104/T)
	Resistance per unit cross-sectional area ${\rho }_{{\rm{cubic}}}$	Ω m2	0.1
Crack propagating	Initial resistance ${R}_{c,0}$	Ω	1.0 × 10−4
	Constant ${C}_{c}$	s−1	7.6 × 10−5

The assumptions made for the model are as follows.

• (1)  
  The side reactions that cause SEI occur uniformly throughout the negative active material in the cell.
• (2)  
  The capacity shift occurs only at the negative electrode.
• (3)  
  The phase transition that causes an increase in internal resistance occurs uniformly throughout the positive active material in the cell.
• (4)  
  The crack propagation that causes an increase in internal resistance occurs uniformly throughout the cell.

The thermal ECM is adopted for modeling the thermal phenomena in the battery pack, as in most system simulations. In this model, the thermal resistance, heat capacity, and heat generation are linked. The temperature of each component is calculated using the energy-balance equation.

$$\begin{eqnarray}&&\begin{array}{c}{C}_{T}\partial {\rm{T}}/\partial {\rm{t}}={Q}_{{\rm{link}}}+{{\rm{Q}}}_{r},\end{array}\end{eqnarray} \tag{ 21 }$$

where C_T is the heat capacity of a cell and is calculated using the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{C}_{T}={\rho }_{{\rm{cell}}}{{\rm{Vol}}}_{{\rm{cell}}}{C}_{p,\mathrm{cell}},\end{array}\end{eqnarray} \tag{ 22 }$$

where ${\rho }_{{\rm{cell}}}$ is the cell density, ${{\rm{Vol}}}_{{\rm{cell}}}$ is the cell volume equal to 4.6 × 10⁻⁴ m³, and ${C}_{p,\text{cell}}$ is the cell-specific heat at a constant pressure. Q_link is the heat transfer into or from the linked components, which can be determined from the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{Q}_{{\rm{link}}}={\rm{\Delta }}{\rm{T}}/{R}_{T},\end{array}\end{eqnarray} \tag{ 23 }$$

where ${\rm{\Delta }}T$ is the temperature difference between the targeted component and the linked components, and R_T is the thermal resistance of a cell, calculated by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{R}_{T}=\displaystyle \frac{{l}_{t}}{{l}_{w}{l}_{d}{k}_{t}},\end{array}\end{eqnarray} \tag{ 24 }$$

where ${l}_{t}$ is the cell thickness equal to 8.0 mm, ${l}_{w}$ is the cell width equal to 260 mm, ${l}_{d}$ is the cell depth equal to 220 mm, and ${k}_{t}$ is the cell thermal conductivity in the thickness direction.

The thermal resistance between the pack case and the outside air R_T,ht , is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{R}_{T,\mathrm{ht}}=1/h{S}_{{\rm{case}}},\end{array}\end{eqnarray} \tag{ 25 }$$

where h is the heat transfer coefficient and ${S}_{{\rm{case}}}$ is the cross-sectional area of the air. The heat transfer coefficient under turbulence flow can be calculated by the following equation: ²⁴

$$\begin{eqnarray}&&\begin{array}{c}Nu=\displaystyle \frac{{\rm{hL}}}{{k}_{{\rm{air}}}}={\rm{0}}{{\rm{.037Re}}}^{4/5}{{\rm{\Pr }}}^{1/3},\end{array}\end{eqnarray} \tag{ 26 }$$

where Nu is the Nusselt number, L is the width of the case, ${k}_{{\rm{air}}}$ is the thermal conductivity of air, Re is the Reynolds number, and Pr is the Prandtl number. For the sake of simplicity, $h$ does not take into account the time dependence and is calculated from the average speed of the BEV.

A schematic of the electrochemical–thermal coupled model for the battery pack is shown in Fig. 2. Each module is connected with the bottom of the pack case whose thermal resistance is negligible. Note that the modules are not connected to the rest of the case and are insulated. This is because this battery pack is mounted on the bottom of the car as shown in Fig. 1a, and is assumed to be in contact with the car's parts intricately in all areas except the bottom. In reality, there is some heat conduction passes between the parts in contact with each other, but for the sake of simplicity, this is not taken into account and is assumed to be insulation. In each cell, the electrochemical models and the thermal ECMs are coupled for temperature and heat generation. The model is constructed by using MATLAB Simulink® R2020a/Simscape. The input parameters of the model are listed in Table VI. Furthermore, the time step was fixed at 1.0 s. This was not the minimum time for convergence, but was set to sufficiently account for the time variation of the heat generation.

Zoom In Zoom Out Reset image size

Table VI. Input parameters for the thermal equivalent circuit model.

Parameter	Unit	Value
Cell density ${\rho }_{{\rm{cell}}}$	kg m−3	2003 20
Cell specific heat ${C}_{p,\text{cell}}$	J kg−1K−1	793 20
Cell thermal conductivity in the thickness direction ${k}_{t}$	W m−1K−1	1 20
Heat capacity of the pack case ${C}_{T,\mathrm{case}}$	J K−1	9.05 × 104a)
Average Reynold number of air Re	—	7.22 × 105b)
Prandtl number of air Pr	—	0.717 25
Heat transfer coefficient $h$	W m−2K−1	40.5 b)

^a): Assumed value. ^b): Calculated with the average speed of the worldwide-harmonized light vehicles test cycles (WLTC) class 3b. ^26,27

The assumptions made for the model are as follows.

• (1)  
  The heat transfer in the direction perpendicular to the cell stacking is negligible.
• (2)  
  The thermal capacity of things, except the cell and case, is negligible.
• (3)  
  The heat is dissipated only to the outside air at the bottom of the case and is insulated from the rest.
• (4)  
  The temperature at the center of the cell is representative of the electrochemical phenomena of the whole cell.
• (5)  
  The heat transfer coefficient to the outside air can be estimated as the average of the driving speed of the BEV.

According to the statistics of EV fire incidents of electric vehicles between 2014 and 2019, internal short circuits (ISCs) accounted for 52% of the accident probability. ²⁸ Therefore, more accurate predictions of ISCs and the development of battery management systems based on them could be invaluable in designing safer battery packs. The main causes of ISCs are mechanical damage, ²⁹ manufacturing defects, ^30,31 overcharge, ³² and overdischarge. ^33,34 ISCs are expected to develop gradually. The early stages of ISCs, such as shorts of 100/10/1 Ω, have remaining time for them to run into thermal accidents because of their low heat generation power. ³⁵ Finally, the resistance of ISCs might become sufficiently small to lead to thermal runaway. In the proposed model, the early stage of ISCs in a cell is modeled to predict the impact of electrochemical and thermal phenomena in the battery pack.

Because the phenomena in an ISC cell are three-dimensionally complicated, it is not appropriate to use ECMs to predict them. Therefore, only the ISC cell is considered a three-dimensional phenomenon by using the finite element method. A schematic of the ISC cell model is shown in Fig. 3. This model is based on the ECM and includes the 3D detailed model of the ISC cell. In this model, the short-circuit object (SCO) is assumed to be in the center position of the cell. In the case of internal short-circuit phenomena, the resistance differs depending on the connection scenario between the Al/Cu foil and positive/negative electrode. ¹⁸ The early stage of ISCs is not assumed to be the scenario of the connection between the Al foil and Cu foil because it has the lowest resistance and is considered to be at risk of thermal runaway. Nevertheless, for the sake of simplicity, the early stage of ISCs is reproduced by adjusting the resistance of the SCO, and the scenario of the connection between the foils is assumed in the model shown in Fig. 3b. Moreover, as shown in Fig. 3c, the ISC cell is assumed to be in position 18. This is because the temperature is the highest in this cell, and the situation is prone to ISCs. Based on previous reports^,18,29 the shape of the SCO is modeled as a cylinder with a diameter of 3 mm.

Zoom In Zoom Out Reset image size

The ISC cell model was constructed based on the nail penetration model proposed by us in a previous report. ²⁰ The electric current through the SCO ${I}_{{\rm{ISC}}}$ is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{I}_{{\rm{ISC}}}=\displaystyle \frac{{\rm{\Delta }}{V}_{{\rm{ISC}}}}{{R}_{{\rm{ISC}}}},\end{array}\end{eqnarray} \tag{ 27 }$$

where ${\rm{\Delta }}{V}_{{\rm{ISC}}}$ is the difference in the electric potential between the Al and Cu foils, and ${R}_{{\rm{ISC}}}$ is the resistance of the SCO. In this model, assuming that the electrical resistance of the current collector foil and tabs is negligible, ${\rm{\Delta }}{V}_{{\rm{ISC}}}$ is calculated using the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}{V}_{{\rm{ISC}}}={V}_{{\rm{cell}}}.\end{array}\end{eqnarray} \tag{ 28 }$$

The above equations, the SP model, and the cycle degradation model are coupled to calculate the electrochemical state of the ISC cell and two parallel-connected cells, as shown in Fig. 3c. The heat generation of the SCO Q_ISC is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{Q}_{{\rm{ISC}}}={R}_{{\rm{ISC}}}{I}_{{\rm{ISC}}}^{2}.\end{array}\end{eqnarray} \tag{ 29 }$$

The thermal phenomena in the ISC cell are presented as a heat transfer equation, which is described by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}\rho {C}_{p}\displaystyle \frac{\partial {\rm{T}}}{\partial {\rm{t}}}+{\rm{\nabla }}{\rm{\Delta }}\left(-{\bf{k}}{\rm{\nabla }}T\right)={q}_{{\rm{gen}}},\end{array}\end{eqnarray} \tag{ 30 }$$

where $k$ is the thermal conductivity vector, and ${q}_{{\rm{gen}}}$ is the heat generation density. The volume integral of ${q}_{{\rm{gen}}}$ is ${Q}_{r}$ in the electrode domain, and ${Q}_{{\rm{ISC}}}$ in the SCO domain. The input parameters in the above equation are the same as those listed in Table VI. As shown in Fig. 3c, the thermal connection on the interface between the ISC cell and the adjacent cell is taken into account by exchanging the heat ${Q}_{{\rm{int}}}$ and the average temperature ${T}_{\mathrm{int},\mathrm{ave}}.$ In other words, the 3D detailed model of the ISC cell and the ECM of the other cells are solved simultaneously, taking into account the connections. The ISC cell model was constructed using COMSOL Multiphysics® 5.6, and the electric–thermal connection was achieved by LiveLink for Simulink.

The calculation meshes of the ISC cell model are shown in Fig. 4. The prism meshes were adopted, and the total number of elements was 42,850. It is confirmed that the calculation results do not vary with more refined meshes. For the same reason, the relative tolerance of the progress for the next time step was set as 1.0 × 10⁻³. This model adopted the backward differentiation formula in COMSOL Multiphysics® and the time step is automatically determined for each step. When the calculation of the thermal ECM advances one step, the data of temperature and heat generations are exchanged if the calculation of above ISC model advances to the same time.

Zoom In Zoom Out Reset image size

The assumptions made for the model are as follows.

• (1)  
  The impact of the current in the current-collecting foils and tabs is negligible because of the lower current flow under the early stage ISCs.
• (2)  
  The current passing through the SCO is not limited to the discharge of the ISC cell, but also includes the discharges of two cells connected in parallel.
• (3)  
  The scenario of the connection between the foils was assumed in the model, and the resistance of the SCO was adjusted to reproduce the early stage ISCs.
• (4)  
  The ISC cell is assumed to be in position No. 18, as shown in Fig. 3c, because the temperature is the highest in this cell and the situation is prone to ISCs.

In this study, the proposed model for the battery pack was not validated based on experimental methods. Battery packs are of the scale actually installed in BEVs, and it is estimated that the experimental measurement of these packs will require a significant effort in terms of equipment. We recognize that this type of validation with experimental approaches is extremely important and should be carried out. Nevertheless, in this study, we focus only on modeling, and this validation remains as future work.

Despite of the above fact, at the cell level, each model that constitutes the proposed model has been validated using the experimental methods reported in previous studies. ^17,20 All input parameters of the proposed model were set based on the prototype cell. The results of the SP model were confirmed to be in good agreement with those of the actual discharge curve and the alternating current impedance in the state of charge (SOC) 100%. ¹⁷ Moreover, the results of the cooperation between the SP model and the cycle degradation model were confirmed to be in good agreement with those of cycle degradation tests at 25 °C and 70 °C. ¹⁷ Finally, the ISC cell model is based on the nail penetration model, and the results of the temperature variation on the cell surface were confirmed to be in good agreement with those of the actual nail penetration tests. ²⁰ Therefore, if the scale effects are negligible, the proposed model for the battery pack may be considered a reasonable model.

The electrical power applied to the battery pack is estimated by the electric consumption model for BEVs published by the Environmental Partnership Council (Tokyo, Japan, https://epc.or.jp/sim_2019model_english, accessed on, December 12/29/2020). In this model, the systems that make up the BEV, such as the vehicle mechanism, motor, and direct current (DC)-DC converter, are built in MATLAB Simulink®. The system simulation was similar to that in several reports. ³⁶ As this study focuses on the battery pack system, the features are not presented.

The typical input parameters for the BEV system simulation are listed in Table VII. In this condition, the calculation results of the electrical power are shown in Fig. 5. These are the electrical power required for the BEV to run at the speed pattern of the worldwide-harmonized light vehicles test cycles (WLTC) class 3b, ^26,27 shown in Fig. 5. Note that positive values of the electrical power indicate the charge, and a negative value indicates the discharge. In this model, regenerative braking is considered, and charging occurs when the vehicle stops.

Zoom In Zoom Out Reset image size

As the electrical power is equal to the values in Fig. 5, the electric current applied to the battery pack ${I}_{\mathrm{app},\mathrm{pack}}$ is calculated by the following equation:

$$\begin{eqnarray}&&\begin{array}{c}{I}_{\mathrm{app},\mathrm{pack}}=\displaystyle \frac{{P}_{{\rm{pack}}}}{{V}_{{\rm{pack}}}},\end{array}\end{eqnarray} \tag{ 31 }$$

where ${P}_{{\rm{pack}}}$ is the electrical power, and ${V}_{{\rm{pack}}}$ is the voltage of the battery pack.

As the standard condition, the calculation is performed as follows:

• (1)  
  The environmental and initial temperatures were 20 °C.
• (2)  
  The initial Li concentrations in all the cells were set as values corresponding to SOC = 100%.
• (3)  
  Initially, all cells were completely free of cycle degradation.
• (4)  
  The power pattern shown in Fig. 5 was applied repeatedly until the calculation was stopped.
• (5)  
  The calculation stops when the battery pack voltage reaches the lower limit of 288 V.

The calculation time for this simulation was 1.7 h using a general computer (CPUs:Intel® Xeon® Gold 6130 CPU@2.10GHz, Memory size:128 GB). The calculation results are presented in Fig. 6. Under this condition, the temperature results for module M1 are shown as representative because all the modules are in the same state. The applied current varied between approximately ±100 A. The voltage of the entire battery pack decreases slowly over time with a variation of approximately ±5 V until it reaches 288 V in 9.34 h. If the operation time is converted into driving distance, it is equivalent to 430 km. This is a large value considering that the driving ranges of the popular EVs released in 2017 or 2018 were 220–430 km. ³⁷ However, 430 km may be a normal BEV specification in 2021 because the range has increased significantly in recent years. For example, the range of the Tesla Model S is estimated to be 647 km.

Zoom In Zoom Out Reset image size

As shown in Figs. 6c and 6d, it was confirmed that the temperature in all the cells increased during driving. Cell No. 18 was the farthest from the outside air, and thus the temperature rise was the largest. At the final time, the temperature of No. 1 is 20.7 °C, that of No. 18 is 22.2 °C, and the average value is 21.7 °C. Note that this calculation model does not consider the contact thermal resistance, and thus the temperature may be lower than the actual value.

Assuming that cell No. 18 of module M1 is the state of the ISCs, the calculations are performed under three conditions: the SCO electric resistance is 100 Ω, 10 Ω, or 5 Ω. The resistance values are in the early stages of ISC. ³⁵ The other conditions were the same as those in the standard condition. The results of the time histories are presented in Fig. 7. From Fig. 7b, at stop time the voltages of module M1 are 54.02 V (standard condition), 53.96 V (100 Ω), 53.20 V (10 Ω) and 51.52 V (5 Ω). These differences are due to self-discharge and are very small for the pack voltage, as shown in Fig. 7a. Therefore, there was no significant difference in the driving ranges presented in Table VIII.

Zoom In Zoom Out Reset image size

From Fig. 7c, average temperature of the shorted cell No. 18 increases during driving and at stop time it reaches 23.9 °C (100 Ω), 27.9 °C (10 Ω), and 33.0 °C (5 Ω). This may be due to the fact that the smaller the short-circuit resistance is, the greater the amount of current carried and the greater the amount of heat generation. As shown in Fig. 7d, the average temperature of the SCO increases immediately after the driving starts and decreases slowly until the stop time. The temperature reaches a maximum of 26.5 °C (100 Ω), 67.1 °C (10 Ω), and 114.2 °C (5 Ω). As the maximum temperature in the case of 5 Ω reaches 90 °C, where the SEI layer begins to decompose and self-heating is initiated, ³⁸ the risk of thermal runaway is high. The SCO resistance, which is the limit of thermal runaway, was considered to be between 5 and 10 Ω. Figure 8 shows the temperature distribution in the shorted cell (No. 18) at the stop time in the case of 5 Ω. The high temperature is near the short-circuit (within a few millimeters), and the distribution is relatively uniform in most of the cells. From the voltage results, it can be observed that an ISC in one cell has only a very small effect on the battery pack voltage as a whole. Nevertheless, the temperature in the vicinity of the SCO is very high, and if thermal runaway occurs, it could lead to a major accident. Therefore, the battery management system, including the detection of internal shorts, must be sufficiently accurate.

Zoom In Zoom Out Reset image size

The results for each cell in the modules at the stop time are shown in Fig. 9. Note that module M1 includes the shorted cell, M2–16 are not shorted, and M2–16 are assumed to be in the same state under these conditions. As shown in Fig. 9a, the temperature in module M1 varies approximately linearly from No. 1 to No. 18 in the cases of ISCs, and the larger the cell number is, the larger the temperature difference compared to the non-shorted case (standard condition). The average temperature rise of No. 18 in module M1 is 2.2 °C (not shorted), 3.9 °C (100 Ω), 7.9 °C (10 Ω), and 13.0 °C (5 Ω), and in the case of 5 Ω, the temperature rise is 6 times higher than that of the not shorted case.

Zoom In Zoom Out Reset image size

From Fig. 9b, it can be observed that the capacity loss during driving increases with a larger number. This is because the SEI layer thickness increases with increasing temperature. The loss capacity of No. 18 in module M1 is 9.3 mAh (standard condition), 9.9 mAh (100 Ω), 12.5 mAh (10 Ω), and 15.9 mAh (5 Ω), and in the case of 5 Ω, the loss capacity is 1.7 times higher than that of the standard condition. As the capacity difference becomes more pronounced with repeated driving, over-discharge is likely to occur. Over-discharge is considered to be one of the factors of ISCs, ^33,34 and if the frequency of this incident increases, the ISC stage may develop.

From Fig. 9c, it can be observed that the Li concentration of the negative electrode at cells No. 16–18 in module M1 is lower than that of the other cells, and it is almost exhausted in the case of 5 Ω. This is because these three cells are connected in parallel and are significantly affected by discharge to the ISCs. Therefore, these cells are prone to over-discharge, which induces ISCs. Thus, it can be observed that an ISC increases the risk of thermal runaway not only in the cell where it occurs, but also in parallel-connected cells.

Using the standard condition as a basis, numerical simulations are performed in which the driving from SOC = 100% to the stop condition is repeated 1000 times. It is assumed that there are no cells in which ISCs occur. The state of cycle degradation at the last time in the run was set as the initial state of the next run. The Li concentration was set as SOC = 100% in the initial state, and the effect of charging from SOC at the stop time to SOC = 100% is not considered.

In this part, comparisons of the internal resistance of each cell is shown. The internal resistance is defined in the following equations:

$$\begin{eqnarray}&&\begin{array}{c}{R}_{{\rm{int}}}=\left({V}_{{\rm{cell}}}-{E}_{\mathrm{eq},p}+{E}_{\mathrm{eq},n}\right)/{I}_{{\rm{app}}}\end{array}\end{eqnarray} \tag{ 32 }$$

where ${R}_{{\rm{int}}}$ is the internal resistance per cell. Since this parameter is time varying, the time average from the start time to 0.5 h, when the first speed pattern of WLTC class 3b of each cycle, is used for comparison.

Figure 10 shows the time histories of the pack voltage for driving at the 1st, 400th, 700th, and 1000th run. From the start of the driving up to 6 h, the behavior was almost the same. After 6 h, the larger the number of runs is, the larger the voltage drop. This is due to the increase in the SEI film and the internal resistance caused by cyclic degradation. As presented in Table IX, the driving range also decreases with the number of runs, and the range of the 1000th run is 15 % lower than that of the standard condition.

Zoom In Zoom Out Reset image size

Figure 11 shows the results for each cell in a module at the stop time for each number of runs. Note that all modules (M1–16) are assumed to be in the same state under these conditions. As shown in Fig. 11a, the larger the number of runs is, the larger the temperature difference. This is because of the increase in heat generation due to the increase in the internal resistance. The temperature difference of cell No. 18 from the initial state is 2.2 °C (1st run), 5.7 °C (400th), 6.8 °C (700th), and 7.2 °C (1000th), and in the case of the 1000th run, the temperature rise is 3.3 times higher than that of the first run. This large difference in the temperature of the cells in the module is considered to cause a difference in the degradation degree, which reduces the overall performance of the module.

Zoom In Zoom Out Reset image size

As shown in Fig. 11b, the larger the number of runs, the larger the loss capacities accumulated from the first run. The accumulated loss capacity of cell No. 18 is 3.8 Ah (400th run), 6.5 Ah (700th), and 9.0 Ah (1000th), and in the case of the 1000th run the capacity is 16 % lower than the initial capacity 57.5 Ah. Moreover, the difference in the accumulated loss capacitance between cells No. 1 and No. 18 is 0.37 (400th run), 0.74 (700th), and 1.10 Ah (1000th), and it increases owing to the increase in the number of runs.

As shown in Fig. 11c, the larger the number of runs is, the larger the internal resistance averaged from the start time to 0.5 h. The internal resistance of cell No. 18 is 0.8 mΩ (1st run), 1.8 mΩ (400th), 2.1 mΩ (700th), and 2.2 mΩ (1000th), and in the case of the 1000th run, the internal resistance is 2.8 times higher than that of the first run.

Using the driving conditions in the previous section as a basis, numerical simulations were performed when the environmental temperature and the initial temperature of the cells were set as 30 °C or 40 °C.

Figure 12 shows the time histories of the pack voltage for each environmental temperature. In the first run, there was almost no difference in voltage, but in the 300th run, the difference became larger, especially at the end of the run. This is because the degree of cycle degradation is more advanced at higher temperatures. As presented in Table X, in the first run, the driving range increases slightly with the environmental temperature because of the increase in the electrochemical reaction characteristics, which induces a decrease in the internal resistance. In this calculation, the electric power applied to the battery was assumed to be 20 °C. Therefore, note that the driving range is assumed to be larger than in the actual condition owing to differences in air conditioning settings. At the 300th run, the driving range decreases with the environmental temperature because of the decrease in the capacity and increase in the internal resistance. The range at 40 °C in the 300th run was 17 % lower than that in the first run. The loss of the driving range is 3.5 times larger than that at 20 °C, which indicates that the cycle degradation progresses approximately 3.5 times faster.

Zoom In Zoom Out Reset image size

Figure 13 shows the results for each cell in a module at the stop time for each environmental temperature. Note that all the modules (M1–16) are assumed to be in the same state under these conditions. As shown in Fig. 13a, in the first run, the larger the environmental temperature is, the smaller the temperature difference. This is because of the decrease in the internal resistance due to the temperature increase. Nevertheless, at the 300th run, this trend was reversed owing to the effects of cycle degradation.

Zoom In Zoom Out Reset image size

As shown in Fig. 13b, the higher the environmental temperature is, the larger the loss capacities accumulated from the first run. The accumulated loss capacity of No. 18 at the 300th run is 2.8 Ah (20 °C), 6.0 Ah (30 °C), and 11.3 Ah (40 °C), and in the case of 40 °C, the capacity is 20 % lower than the initial capacity of 57.5 Ah.

From Fig. 13c, the trend of the internal resistance is as described above. At the 300th run, the internal resistance of cell No. 18 was 1.3 mΩ (20 °C), 1.9 mΩ (30 °C), and 2.0 mΩ (40 °C), and in the case of 40 °C, the internal resistance was 4.4 times higher than that of the first run. Thus, it can be confirmed that the environmental temperature has a significant effect on capacity loss and internal short circuits.