Multi-Scale and Multi-Dimensional Thermal Modeling of Lithium-Ion Batteries ^†

Abstract

In this study, we present a three-dimensional (3-D), multi-scale, multi-physics lithium-ion battery (LIB) model wherein a microscale spherical particle model is applied to an electrode particle domain and a comprehensive 3-D continuum model is applied to a single cell domain consisting of current collectors, porous electrodes, and a separator. Particular emphasis is placed on capturing the phase transition process inside the lithium iron phosphate (LFP) particles that significantly influences the LIB charge and discharge behaviors. The model is first validated against the experimental data measured at various discharge rates. In general, the model predictions compare well with the experimental data and further highlight key electrochemical and transport phenomena occurring in LIBs. Besides elucidating the phase transition evolution inside LFP particles and location-specific heat generation mechanism, multi-dimensional contours of species concentration, temperature, and current density are analyzed under a 3-D cell configuration to provide valuable insight into the charge and discharge characteristics of LIBs. The present multi-scale LIB model can be applied to a realistic LIB geometry to search for the optimal design and operating conditions.

1. Introduction

Since Padhi et al. [1] first introduced lithium iron phosphate (LiFePO₄, LFP) as the positive electrode material for lithium-ion batteries (LIBs), LFP has been regarded as an excellent positive electrode material for large-scale LIBs owing to several attractive features, such as the high theoretical capacity (170 mAh/g), good thermal and electrochemical stabilities, reasonable cycle life, low cost, and low toxicity. Moreover, LFP exhibits a relatively low flat open-circuit voltage around 3.4 V vs. Li/Li+, which tends to alleviate side reactions and electrolyte decomposition. However, LFP suffers from a poor rate capability ascribed to poor transport properties inside the LFP particles, including low electronic conductivity, poor lithium ion diffusivity, and sluggish phase transformations between the lithium-deficient phase (α phase: LiaFePO₄) and lithium-rich phase (β phase: Li_1−bFePO₄). Various electrode designs and fabrication techniques have been attempted to improve the electron and/or lithium ion transport inside the LFP electrode [2,3,4,5,6,7,8].

To enhance the performance of the LFP electrode, a parallel effort has been made in LIB modelling and simulations [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]. Srinivasan et al. [9] developed a mathematical discharge model for the LFP electrode based on the idea of a shrinking core. They insisted that transport limitations led to loss in the utilization of the core in active materials, which could be improved by small particle sizes and a high diffusivity. However, they neglected to include the diffusion of the α phase, and their model was only effective for discharge. Furthermore, the model is not able to capture inconsistent phase transition characteristics in LFP electrodes such as charge/discharge asymmetry and path dependence. Later, Kasavajjula et al. [10] upgraded the shrinking core model by considering lithium ion diffusion in both the α and β solid phases and the interface mobility, which is a mixed-mode phase transformation. They found that the rate capability of the phase transformation electrode material mainly depends on lithium ion diffusion; thus, a high rate capability can be achieved by reducing the particle size. Although the shrinking core model is known to have limitations in predicting the effects of multiple cycles and inconsistent phase transition characteristics, such as charge/discharge asymmetry, path dependence, and tangential front propagation [11,12], recently, Khandelwal et al. [13] further improved a shrinking core model and presented a generalized core-shell model for LIBs with LFP that could capture the phenomena of path dependence, charge/discharge asymmetry, and tangential front propagation caused by phase transitions in the electrode active material. Later, the energy equation was additionally included in the model to consider non-isothermal effects during cell operation [14]. They compared individual heat generation terms, including the heat release due to the phase transition. Meanwhile, the geometrical characteristics of LIBs are important because the large size and extreme operating conditions yield non-uniform thermal and electrical behaviors.

Besides the shrinking core and core-shell models, there are electrochemical models for LFP, such as the phase field model [13,14,15] and resistive reactant model [16,17,18]. The phase field model can capture the more general bulk transport limitation compared to the shrinking core model, though it requires computational elaboration. The resistive reactant model has the advantages of describing multi-cycle operation and features of charge/discharge asymmetry, even if it neglects phase transition and requires several parameters. Andersson and Thomas [19] suggested two different LIB models with LFP electrode. One model is a radial model that was radius-dependent process and another model is a mosaic model that consider many particles. These models allow to consider the inactive LFP regions. Laffont et al. [20] observed the juxtaposition of the LFP by using high-resolution electron energy loss spectroscopy (HREELS) and indicated that the classical shrinking core model is not able to capture the effect of lithium reactivity mechanism in LFP. Meethong et al. [21] found that the miscibility gap between the lithium-deficient phase and lithium-rich phase in the LFP particle was reduced as particle size decreased. Delmas et al. [22] described the LFP reaction mechanism using a domino-cascade model that accounts for the coexistence of intercalated and deintercalated individual particles. They suggested that the minimization of the elastic energy improved the reaction process in LFP electrode. Ramana et al. [23] presented a cooperative phase boundary motion model, showing that the mechanism of phase transition in LFP particles is similar with the motion of nucleation fronts and dislocations in superplastic alloys. Liu et al. [24] proposed a dual-phase solid-solution reaction mechanism based on the observation of nonequilibrium lithium insertion and extraction in the LFP. Xu et al. [27] considered the effects of the tab placement on the commercial LP2770120 prismatic LFP/graphite battery by treating the tabs as three-dimensional (3-D) and local cell units as one-dimensional (1-D). Unfortunately, they neglected the phase transition in the active material mentioned above. Panchal et al. [28,29] analyzed thermal behaviors of large-scale LIBs with LFP electrodes, i.e. large sized 20 Ah-LFP prismatic battery cell and 18650-LFP cylindrical battery cell under various C-rates. Jiang et al. [30] developed a reduced low-temperature electro-thermal coupled model and proposed an efficient charging strategies in order to evaluate the condition of battery health. Lin et al. [31] investigated the electrochemical and thermal behavior of a prismatic LIB with LFP under a mixed charge-discharge cycle by using the calorimetric technique. They indicated that simulation results based on calorimetrically measured value was well matched with experimental data.

The previous studies clearly demonstrate that the micro-scale phase transition phenomena inside LFP particles significantly affect the macro-scale LIB performance, particularly in terms of the rate capability, degree of charge/discharge asymmetry, and path dependent response under charge/discharge cycles. Therefore, the phase transition effect should be considered in macroscale LIB modeling and simulations. While a number of LIB models can predict many details, most of the phase transition modeling and simulation studies have been conducted based only on a 1-D mathematical formulation decoupled from the multi-dimensional LIB analysis. In this study, we present a multi-scale LIB model wherein a microscale spherical model for an electrode domain is coupled with a comprehensive 3-D macroscale LIB model to deal with a full cell domain comprising electrodes, a separator, current collectors, etc. To capture the phase transition processes inside the LFP particles, the numerical scheme of the Landau transformation [32] is adopted to effectively resolve the moving boundary problem in the fixed mesh. The multi-scale LIB model is first validated against experimental data measured at different discharging rates. Then, the detailed effects of the phase transition parameters, such as the transport properties and equilibrium concentrations for the α and β phases on macroscale cell performance and distributions are studied via multi-scale LIB simulations. This study clearly indicates that multi-scale and two-phase LIB simulations are indeed necessary to faithfully capture the discharge-charge characteristics of LIBs with LFP positive material.

2. Numerical Model

2.1. Macroscale LIB Model

A 3-D, macroscale, transient, non-isothermal, electrochemical transport-coupled LIB model was developed using the lithium-air battery model presented in our previous work [33]. The model that deals with the coupling of multi-dimensional transport phenomena with electrochemical kinetics can be applied to a single cell with the negative electrode, separator, positive electrode, current collectors, and tabs, and simulated to predict the macroscale dynamic behavior of LIBs during charging and discharging. In addition, we also developed a microscale electrode particle model (described in detail in Section 2.2), aiming to capture the microscale physical phenomena occurring inside the electrode particles, such as solid-state diffusion, ohmic/mass transport losses, and phase-transition processes. The microscale electrode model is incorporated into the macroscale full cell model to fully elucidate the micro-macroscale coupled phenomena during LIB charging and discharging.

Figure 1 schematically shows the multi-scale LIB model, main electrochemical reactions, and resultant species/charge transport during the charge/discharge of a LIB. The cell is composed of five components such as a negative/positive electrode and a copper/aluminum current collector, and a separator. The electrode consists of active material, electrolyte, polymer binder and conductive filler. The active materials of anode and cathode are graphite meso-carbon microbeads and LFP, respectively. The electrolyte is commonly lithium hexafluorophosphate (LiPF₆) in ethylene carbonate:dimethyl carbonate. The porous separator is made up of polymer matrix and liquid electrolyte. The lithium intercalation/deintercalation processes occurring in each electrode are as follows:

$$\mathrm{Negative}\mathrm{electrode}:{\mathrm{Li}}_x{\mathrm{C}}_6\underset{\mathrm{charge}}{\overset{\mathrm{discharge}}{\rightleftarrows }}{\mathrm{Li}}_{x-z}{\mathrm{C}}_6+z{\mathrm{Li}}^{+}+z{\mathrm{e}}^{-}$$

(1)

$$\mathrm{Positive}\mathrm{electrode}:{\mathrm{Li}}_{y-z}{\mathrm{FePO}}_4+z{\mathrm{Li}}^{+}+z{\mathrm{e}}^{-}\underset{\mathrm{charge}}{\overset{\mathrm{discharge}}{\rightleftarrows }}{\mathrm{Li}}_y{\mathrm{FePO}}_4$$

(2)

The assumptions made in developing the model are as follows:

(1)

The concentrated binary electrolyte is assumed

(2)

Side reactions that are significant only at relatively high temperatures are ignored

(3)

Lithium intercalation/deintercalation kinetics are assumed to be first order in the reactants and described by the standard Butler-Volmer equation

(4)

The active materials of the electrode are made up of spherical particles with a uniform size

(5)

The volume change during cell operation is insignificant, resulting in constant electrode porosities

(6)

Multi-particle behavior and anisotropic propagation of reaction interface in LFP particles are neglected although these behaviors have shown [34,35,36,37]

Under the aforementioned assumptions, the proposed model for the LIB cell is governed by the conservation of mass, charge, and energy equations:

Mass conservation:

$$\mathrm{Lithium} \mathrm{ion} \mathrm{in} \mathrm{the} \mathrm{electrolyte} \mathrm{phase}: \frac{\partial \left({\epsilon }_e{C}_e\right)}{\partial t}=\nabla \cdot \left(D_e^{eff}\nabla C_e\right)+\frac{j}{F}-\frac{t_{+}}{F}j-\frac{i_e}{F}\nabla \cdot \left(t_{+}\right)$$

(3)

In Equation (3), j denotes the local charge transfer current density that represents the reaction rates of the lithium intercalation and deintercalation processes in Equations (1) and (2). These kinetic expressions are described by the Butler-Volmer equation as follows:

$$j=a_s{i}_0\left(\exp \left(\frac{{\alpha }_{a}F}{RT}\eta \right)-\exp \left(-\frac{{\alpha }_{c}F}{RT}\eta \right)\right)$$

(4)

where α_a and α_c are the anodic and cathodic transfer coefficients of the electrode reaction. The exchange current density, i₀, depends strongly on the local lithium concentrations (in both the electrolyte and solid phases) and the temperature,

$$i_0=k{\left(C_e\right)}^{{\alpha }_a}{\left(C_{s,\max }-C_s\right)}^{{\alpha }_a}{\left(C_s\right)}^{{\alpha }_c}$$

(5)

In Equation (5), C_e and C_s refer to the bulk concentration of lithium ions in the electrolyte within the porous electrode and the solid-state lithium concentration on the interfacial surface between the solid electrode and liquid electrolyte, respectively. The local value of C_s is calculated by a microscopic solid-state lithium diffusion model that will be described in Section 2.2. For an electrode composed of spherical particles of radius r_s, the specific interfacial area, a_s, can be written as:

$$a_s=\frac{3{\epsilon }_s}{R_p}$$

(6)

where ε_s is the porosity of the electrode. The local surface overpotential, η, is defined as:

$$\eta ={\varphi }_s-{\varphi }_e-U$$

(7)

$$U=-\frac{\Delta G}{nF}=U_{ref}+\left(T-T_{ref}\right)\frac{\partial U}{\partial T}$$

(8)

where U is the open circuit potential (OCP) of the lithium intercalation and deintercalation reactions, determined from the state of charge (SOC) and temperature. The detailed expressions of OCP for graphite and the LFP electrode materials are given in Figure A1 and Figure A2, respectively.

In the Equation (3), t₊ is the transference number of the lithium ion with respect to the velocity of the solvent. The effective diffusion coefficient, $D_e^{eff}$ is calculated by the Bruggeman relation:

$$D_e^{eff}=D_e{\epsilon }_e^{1.5}$$

(9)

$$D_e=5.34\times {10}^{-10}\exp \left(-1.565\times {10}^{-4}{C}_e\right)$$

(10)

Charge conservation:

$$\mathrm{In} \mathrm{the} \mathrm{electrolyte} \mathrm{phase}: \nabla \cdot \left({\kappa }^{eff}\nabla {\varphi }_e\right)+\nabla \cdot \left(\frac{2RT{\kappa }^{eff}}{F}\left(t_{+}^0-1\right)\left(1+\frac{d\ln f_{\pm }}{d\ln C_e}\right)\nabla {\mathrm{lnC}}_e\right)+j=0$$

(11)

$$\mathrm{In} \mathrm{the} \mathrm{solid} \mathrm{phase}: \nabla \cdot \left({\sigma }^{eff}\nabla {\varphi }_s\right)-j=0$$

(12)

where ϕ_s and ϕ_e denote the electronic and electrolyte potentials, respectively and f is the mean molar activity coefficient of the electrolyte. The intrinsic conductivity of the lithium ion through the electrolyte consisting of LiPF₆ in ethylene carbonate:dimethyl carbonate is a strong function of the lithium ion concentration [9]:

$$\kappa =0.0911+1.901\times {10}^{-3}{C}_e-1.052\times {10}^{-6}{C}_e+0.1554\times {10}^{-9}{C}_e$$

(13)

For the porous components of LIBs, the above expression is further modified by the Bruggeman factor based on the electrolyte fraction, ε_e, as follows:

$${\kappa }^{eff}={\epsilon }_e^{1.5}\kappa$$

(14)

Energy conservation:

$$\frac{\partial \left(\overline{\rho }{\overline{c}}_{p}T\right)}{\partial t}=\nabla \cdot \left(\varpi \nabla T\right)+Q_{heat}$$

(15)

where the heat capacity, ρc_p, and the effective thermal conductivity, ω^eff, can be defined from the corresponding component values:

$$\overline{\rho }{\overline{c}}_p={\displaystyle \sum _i{\epsilon }_i{\rho }_i{c}_{p,k}}$$

(16)

$${\omega }^{eff}={\displaystyle \sum _i{\epsilon }_i{\omega }_i}$$

(17)

The source term of the energy conservation equation, Q_heat, accounts for the heat generation during LIB charging and discharging, and is expressed as:

$$Q_{heat}=\underset{\mathit{Irreversible}\mathit{heat}}{\underbrace{\left(j\cdot \eta \right)+\left({\sigma }^{eff}\nabla {\varphi }_s\cdot \nabla {\varphi }_s\right)+\left({\kappa }^{eff}\nabla {\varphi }_e\cdot \nabla {\varphi }_e+{\kappa }_D^{eff}\nabla \ln c_e\cdot \nabla {\varphi }_e\right)}}\underset{\mathit{Reversible}\mathit{heat}(-\mathit{discharging},+\mathit{charging})}{\underbrace{\pm j\cdot T\frac{\partial U}{\partial T}}}$$

(18)

In Equation (18), the first, second, and third terms that represent the irreversible heat are due to lithium intercalation/deintercalation processes, the Joule heating in the solid phase, and the Joule heating in the electrolyte phase, respectively. The fourth term results from the reversible heat release associated with the entropy change of the lithium intercalation/deintercalation processes.

In addition, we consider Arrhenius type temperature dependence terms for the ionic conductivity, species diffusivity, and exchange current density by using $\Omega$ indicates a general parameter:

$$\Omega ={\Omega }_{ref}\exp \left(\frac{E_{act,\Omega }}{R}\left(\frac{1}{T_{ref}}-\frac{1}{T}\right)\right)$$

(19)

where the term E_act,_Ω denotes the corresponding activation energy for the general parameter.

2.2. Microscale Model for Electrode Particles

The solid-state lithium diffusion equation is described in a spherical coordinate system for an electrode composed of spherical particles to account for the effects of microscale lithium transport inside the electrode particles:

$$\mathrm{Lithium} \mathrm{ion} \mathrm{in} \mathrm{the} \mathrm{solid} \mathrm{phase}: \frac{\partial C_s}{\partial t}=\frac{D_s}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial C_s}{\partial r}\right)$$

(20)

where D_s is the solid phase diffusion coefficient for either the negative or positive electrode material. The relevant boundary condition for individual electrode particles can be written as follows:

$${D_s\frac{\partial C_s}{\partial r}|}_{r=0}=0,{D_s\frac{\partial C_s}{\partial r}|}_{r=R_p}=-\frac{i}{F}$$

(21)

where i denotes the local current density, which varies as a function of the local state of charge or discharge and is calculated by the macroscopic LIB model.

As seen in Equation (1), graphite (Li_xC₆) is used as the active material of negative electrodes; thus, Equation (20) is simply solved without any further numerical treatment. However, the LFP particles in the positive electrode undergo a phase-transition between the α and β phases during the lithium intercalation/deintercalation processes, which induces the numerical issue of moving the phase boundary within the LFP particle. To effectively resolve the moving boundary problem in a fixed mesh, we adopted the numerical scheme of the Landau transformation [32]. A detailed derivation is given below.

Equation (20) can be rewritten in terms of the two-phase interfacial position, s, as a function of charge or discharge time, t.

In the core region:

$$r^{\lambda -1}\frac{\partial C\left(r,t\right)}{\partial t}=\frac{\partial }{\partial r}\left(r^{\lambda -1}{D}_{\alpha }\left(C\left(r,t\right)\right)\frac{\partial C\left(r,t\right)}{\partial r}\right),0\le r\le s\left(t\right)$$

(22)

In the shell region:

$$r^{\lambda -1}\frac{\partial C\left(r,t\right)}{\partial t}=\frac{\partial }{\partial r}\left(r^{\lambda -1}{D}_{\beta }\left(C\left(r,t\right)\right)\frac{\partial C\left(r,t\right)}{\partial r}\right),s\left(t\right)\le r\le R_p$$

(23)

Moving boundary condition at the interface:

$${D_{\alpha }\left(C\left(r,t\right)\right)\frac{\partial C\left(r,t\right)}{\partial r}|}_{r=s{\left(t\right)}^{-}}-{D_{\beta }\left(C\left(r,t\right)\right)\frac{\partial C\left(r,t\right)}{\partial r}|}_{r=s{\left(t\right)}^{+}}=\left[C_{\beta }^{eq}-C_{\alpha }^{eq}\right]\frac{ds\left(t\right)}{dt}$$

(24)

where λ is a parameter that describes the geometry of the domain (planar (λ = 0), cylindrical (λ = 1), or spherical (λ = 2)). Equation (24) can be derived by assuming that there is a local equilibrium between the α and β phases under the principle of solute conservation. Equations (22)–(24) can be reformulated using the Landau transformation wherein the radial distance, r, and a solid phase lithium concentration, C, are transformed into a proportional position (u for the core region and v for the shell region) and location-specific lithium concentration (p for the core region and q for the shell region), respectively:

In the core region:

$$\begin{array}{l}\frac{\partial \left(ps{\left[us\right]}^2\right)}{\partial t}=\frac{\partial }{\partial u}\left({\left[us\right]}^2\left\{\dot{s}pu+\frac{D_{\alpha }}{s}\frac{\partial p}{\partial u}\right\}\right),0\le u\le 1\\ u=\frac{r}{s},p=p\left(u,t\right)=C\left(r,t\right)\\ B.C.\\ {D_{\alpha }\frac{\partial p}{\partial u}|}_{u=0}=0,{p|}_{u=1}=C_{\alpha }^{eq}\end{array}$$

(25)

In the shell region:

$$\begin{array}{l}\frac{\partial \left(q\left(R_p-s\right){\left[v\left(R_p-s\right)+s\right]}^2\right)}{\partial t}\\ =\frac{\partial }{\partial v}\left({\left[v\left(R_p-s\right)+s\right]}^2\left\{\dot{s}q\left(1-v\right)+\frac{D_{\beta }}{\left[R_p-s\right]}\frac{\partial q}{\partial v}\right\}\right),0\le v\le 1\\ v=\frac{r-s}{R_p-s},q=q\left(v,t\right)=C\left(r,t\right)\\ B.C.\\ {q|}_{v=0}=C_{\beta }^{eq},{D_{\beta }\frac{\partial q}{\partial v}|}_{v=1}=\left(R_p-s\right)\frac{j_n}{F}\end{array}$$

(26)

Moving boundary condition at the interface:

$${\frac{D_{\alpha }}{s}\frac{\partial p}{\partial u}|}_{u=1}-{\frac{D_{\beta }}{R_p-s}\frac{\partial q}{\partial v}|}_{v=0}=\left[C_{\beta }^{eq}-C_{\alpha }^{eq}\right]\frac{ds}{dt},u=1;v=0$$

(27)

In Equations (25)–(27), ṡ is an interface position velocity representing the interface position change with respect to time. Although Equations (25)–(27) appear more complex than Equations (22)–(24), they have simplified the problem in that all of the boundaries are now fixed. Consequently, we can solve the problem by applying implicit schemes that are not subject to any numerical restrictions.

2.3. Multi-Scale Modeling and Numerical Implementation

Figure 2 summarizes the multi-scale modeling and simulation approach used in this study wherein two computational domains at the electrode particle and cell levels are built on different coordinate systems. The 3-D transient, non-isothermal, macro-scale LIB model described in Section 2.1 is applied to the 3-D LIB cell domain for electrochemical-transport-thermal-coupled multi-dimensional simulations, while the solid-phase lithium transport in the electrode particles is predicted by the micro-scale electrode model within the particle domain. In this study, the 3-D macro-scale LIB model based on Cartesian coordinate system is coupled with the micro-scale electrode model based on spherical coordinate system. Then, both models defined in the different type and scale domains were implemented numerically into a commercial computational fluid dynamics program, ANSYS Fluent based on its user-defined functions and solved simultaneously. The model variable solutions in individual domains and the coupling variables exchanged during multiscale simulations are illustrated in Figure 2. The conservation equations are discretized by the finite volume method with the semi-implicit pressure linked equation (SIMPLE) algorithm and implicit time integration. The discretized finite volume equations are solved by the algebraic multi-grid (AMG) method.

The convergence criteria were set to 10⁻⁶ for the equation residuals. The grid and time step sizes are inversely proportional to the spatial and temporal gradients of variables calculated by the LIB simulations, and thus were properly chosen based on several evaluation profiles during discharge. The maximum time step size to ensure adequate time resolution is 1.0 s.

According to our grid-independent test, the 70.35 cm² computational cells require around 6.7 million grid points, which is computationally challenging for transient LIB simulations. To reduce the computational turnaround time, the 3-D LIB code was massively parallelized for parallel computing.

2.4. Initial and Boundary Conditions

Figure 3 shows the computational domain generated for the 1 cm² experimental test cell along with relevant dimensions and boundary conditions to validate the multiscale LIB model against the experimental data of Srinivasan et al. [9]. An orthogonal but non-uniform grid is written on the Cartesian coordinate system for macroscale LIB simulations. Meanwhile, the microscale electrode particle model is applied to every grid point of the macroscale cell domain and internally solved based on the macroscale variable solutions predicted by the 3-D macroscale LIB model. Detailed cell dimensions, initial conditions, and operating conditions are listed in

Table 1. Geometric, initial, and operating conditions.

Description	Value	Ref.
Thickness of the separator, δsep	20 μm	[38]
Thickness of the negative electrode, δNE	59 μm	[38]
Thickness of the positive electrode, δPE	62 μm	[9]
Thickness of the Cu current collector in the negative side, δCu	9 μm	[38]
Thickness of the Al current collector in the positive side, δAl	16 μm	[38]
Spherical particle size for negative electrode (RNE,p)	14.75 μm	[38]
Spherical particle size for positive electrode (RPE,p)	52 nm	[39]
Initial concentration of the solid phase lithium concentration in the negative electrode for discharge, C0s,NE,discharge	25,096 mol∙m−3	[14]
Initial concentration of the solid phase lithium concentration in the positive electrode for discharge, C0s,PE,discharge	230 mol∙m−3	[14]
Initial concentration of the electrolyte phase lithium concentration, C0e	1200 mol∙m−3	[14]
Initial concentration of the solid phase lithium concentration in the negative electrode for charge, C0s,NE,charge	12,548 mol∙m−3	Assumed
Initial concentration of the solid phase lithium concentration in the positive electrode for charge, C0s,PE,charge	19,902.5 mol∙m−3	Assumed
Initial temperature (T0)	298.15 K
Operating pressure (P)	1 atm
Volume fraction for the negative electrode, εe,NE, εs,NE, εf,NE, εNE	0.017, 0.62, 0.013, 0.35	[40]
Volume fraction for the positive electrode, εe,PE, εs,PE, εf,PE, εPE	0.065, 0.56, 0.025, 0.35	[40]
Volume fraction for the separator, εe,sep, εsep	0.5, 0.5	[40]

. In addition, the kinetic, transport, and physiochemical properties used for the LIB simulations are given in

Table 2. Kinetics, physiochemical, and transport properties.

Description	Value	Ref.
Diffusion coefficient of the solid phase in the negative electrode, Ds	9.0 × 10−14 m2∙s−1	[13]
Diffusion coefficient of the α solid phase in the positive electrode, Dα	1.184 × 10−18 m2∙s−1	[13]
Diffusion coefficient of the β solid phase in the positive electrode, Dβ	3.907 × 10−19 m2∙s−1	[13]
Conductivity of the solid active materials in the negative electrode, σs,NE	2 S∙m−1	[27]
Conductivity of the solid active materials in the positive electrode, σs,PE	0.01 S∙m−1	[27]
Conductivity of the Cu current collector for negative, σs,Cu	6.33 × 107 S∙m−1	[27]
Conductivity of the Al current collector for positive, σs,Al	3.83 × 107 S∙m−1	[27]
Thermal conductivity of the separator, ksep	0.334 W∙m−1∙K−1	[27]
Thermal conductivity of the negative electrode, kNE	1.04 W∙m−1∙K−1	[27]
Thermal conductivity of the positive electrode, kPE	1.48 W∙m−1∙K−1	[27]
Thermal conductivity of the Cu current collector for negative, kCu	400 W∙m−1∙K−1	[27]
Thermal conductivity of the Al current collector for positive, kAl	160 W∙m−1∙K−1	[27]
Density of the separator, ρsep	1210 kg∙m−3	[27]
Density of the negative electrode, ρNE	2660 kg∙m−3	[27]
Density of the positive electrode, ρPE	1500 kg∙m−3	[27]
Density of the Cu current collector for negative, ρCu	8900 kg∙m−3	[27]
Density of the Al current collector for positive, ρAl	2700 kg∙m−3	[27]
Heat capacity of the separator, cp,sep	1518 J∙kg−1∙K−1	[27]
Heat capacity of the negative electrode, cp,NE	1437.4 J∙kg−1∙K−1	[27]
Heat capacity of the positive electrode, cp,PE	1260.2 J∙kg−1∙K−1	[27]
Heat capacity of the Cu current collector for negative, cp,Cu	385 J∙kg−1∙K−1	[27]
Heat capacity of the Al current collector for positive, cp,Al	903 J∙kg−1∙K−1	[27]
Activation energy for exchange current density for negative	4000 J∙mol−1	[40]
Activation energy for exchange current density for positive	20,000 J∙mol−1	[40]
Activation energy for solid phase Li diffusion coefficient	4000 J·mol−1	[40]
Activation energy for ionic conductivity of electrolyte	4000 J·mol−1	[40]
Anodic and cathodic transfer coefficients for an electrode reaction, αa, αc	0.5, 0.5	[21]
Transfer number, t+	0.363	[20]
Faraday constant, F	96,485 C∙mol−1
Universal gas constant, R	8.3143 kJ∙kmol−1∙K−1
Maximum solid phase concentration for negative, Cs,max,NE	31,370 mol∙m−3	[20]
Maximum solid phase concentration for positive, Cs,max,PE	20,950 mol∙m−3	[20]
Equilibrium concentration for the α phase, Ceqα	0.048Cs,max,NE	[19]
Equilibrium concentration for the β phase, Ceqβ	0.89Cs,max,PE	[10]

.

LIBs can be operated under either a constant current or constant voltage mode during the charging and discharging processes. Therefore, the constant current or constant voltage condition below is applied to each tab or sidewall of the current collectors:

$$\mathrm{Negative}\mathrm{current}\mathrm{collector}:{\varphi }_s=0$$

(28)

$$\begin{array}{lll}\mathrm{Positive}\mathrm{current}\mathrm{collector}:& {\varphi }_s=V_{cell}& \mathrm{for}\mathrm{potentiostatic}\mathrm{mode}\\ & {\sigma }^{eff}\nabla {\varphi }_s=\frac{{\overrightarrow{i}}_{cell}}{A_{cc}}& \mathrm{for}\mathrm{galvanostatic}\mathrm{mode}\end{array}$$

(29)

where ${\overrightarrow{i}}_{cell}$ denotes the current applied to the cell. The other external surfaces, including the top, bottom, front, and rear surfaces of the computational domain, are assumed to be electrically insulated. For thermal boundary conditions, the adiabatic boundary condition is applied to the top and bottom surfaces of the computational domain, whereas the convective boundary condition is applied to the side walls of the current collectors to approximate single-cell experimental conditions:

$$-{\omega }^{eff}\frac{\partial T}{\partial \overrightarrow{n}}=h\left(T_{cell}-T_{amb}\right)$$

(30)

where $\overrightarrow{n}$ and T_amb denote the unit normal vector out of the LIB current collectors and ambient temperature, respectively.

Assuming that a LIB is initially in a state of thermodynamic equilibrium, the initial conditions for the charge or discharge processes are defined as:

$$C_{s,ne}=C_{s,ne}^0,C_{s,pe}=C_{s,pe}^0,C_e=C_e^0,T=T^0$$

(31)

3. Results and Discussion

The 3-D multi-scale LIB model described in Section 2 was first validated against the experimental data measured by Srinivasan et al. [9], where a half cell (LFP vs. Li) was discharged to 1.09% SOC from 100% SOC. It should be noted that the positive electrode potential was taken from the LIB full cell simulation results for a comparison purpose. As seen in Figure 4, the simulated discharge curves generally agree well with the measured curves at the various C-rates ranging from C/2 (0.26 mA/cm²) to 2C (2.6 mA/cm²). The C-rate is defined as the discharge current divided by the theoretical current driven under which the battery would deliver its nominal rated capacity in one hour (nominal capacity), i.e.,:

$$\mathrm{C}-\mathrm{rate} = \frac{\mathrm{discharge}\mathrm{current}}{\mathrm{nominal}\mathrm{capacity}}$$

(32)

The comparison indicates that the capacity fade behavior with increasing C-rate, i.e., the characteristics of the LFP electrode, is accurately captured by the present LIB model. The agreement is relatively poor at 2C, where the experimental data show the steeper voltage drop than the simulation results. This deviation indicates that the discharge reaction and phase transition in the LFP particles are possibly more abrupt and complicated at the higher C-rate that seems to be underestimated by the present model. The evolution of the two-phase interface between the α and β phases and the lithium ion concentration inside the LFP particles are plotted in Figure 5 as a function of discharge time. At a higher C-rate, the two-phase interface inside the LFP particles moves faster toward the core but ends earlier because lithium ion transport and the phase transition inside the LFP particles are rate-limiting mechanisms. For instance, Figure 5a,b show that the final location of the two-phase interface at 2C is near the outer region of r/r₀ = 0.8 inside the LFP particles, which implies that a high degree of utilization loss occurs at the high C-rate. In contrast, the LFP particles are almost fully utilized at the lower C-rate, C/2, wherein the two-phase interface almost reaches the particle center (Figure 5a), and most areas of the LFP particles are transformed to the β phase (Figure 5d).

The simulation results during charging are shown in Figure 6, where using a higher C-rate results in a higher voltage rise and produces the lower charge capacity. The dimensionless interfacial position and lithium ion concentration profiles inside the LFP particles in Figure 7 clearly show the detailed behaviors of lithium ion transport and two-phase transition inside the particles during charging. When the higher C-rate is applied, the movement of the two-phase interface toward the core becomes faster and finishes earlier because of the more severe limitation of lithium ion transport and phase transition inside the LFP particles (Figure 7a). Delithiation occurs in the LFP particles during charging; thus, the lithium ion concentration profiles in Figure 7b,d appear opposite to those during the discharging processes. The phase transition occurs from the β phase to the α phase; thus, the α phase starts to appear from the outer region of the LFP particles. As the charging process proceeds, the more severe limitation of lithium ion transport and phase transition is observed at the higher C-rate, showing that a larger inner area of particles is still present in the β phase at the latter charge stages.

Various heat generation sources, including irreversible heat of the lithium intercalation/deintercalation, reversible entropic heat, and Ohmic Joule heating, are plotted in Figure 8 and Figure 9 (for discharging and charging) wherein the volumetric and location-specific heat source terms of energy equation in Equation (18) are integrated along the cell thickness direction and then averaged over the entire cell area. While all the terms are proportional to the C-rate, the irreversible heat source from lithium intercalation/deintercalation is dominant among the three heat sources. In particular, the irreversible heat sources in the positive LFP electrode show an increasing behavior as charging or discharging progresses. This indicates that during the galvanostatic charge or discharge mode, the overpotential in the LFP electrode substantially increases with time because of its phase-change characteristics; thus, more heat is irreversibly released at the latter charge/discharge stages.

The entropic heat, i.e., the second largest heat source among the three, is purely reversible; thus, its sign changes for charging and discharging. The sign of the entropic heat terms in the graphite carbon and LFP electrodes is determined by the change of entropy during its reaction, which is equivalent to dU/dT plotted as a function of θ in Figure A3 and Figure A4, respectively. As shown in Figure A3, dU/dT for the graphite carbon negative electrode is negative over the entire range of θ. Therefore, the entropic heat in the negative electrode is endothermic (cooling) for discharging (Figure 8b) and exothermic (heating) for charging (Figure 9b). In contrast, dU/dT for the LFP electrode (Figure A4) is negative in the range of 0.2 ≤ θ ≤ 1.0 but becomes positive when a small amount of lithium (θ < 0.2) is present in the LFP particles. As a result, the entropic heat terms observed at very early discharge and charge stages are endothermic (negative) but rapidly become exothermic (positive), which is indicative of the fast phase transition on the external LFP surfaces, i.e., from the α to β phase during discharging and from the β to α phase during charging.

The simulation results discussed above clearly show that the performance of the LIB with the LFP electrode is mainly limited by solid-state lithium diffusion and phase transition phenomena inside the LFP particles. Therefore, the parametric study is conducted to assess the effect of the lithium diffusivities of the α and β phases on the discharge voltage response and heat generation behavior. Figure 10a shows that the discharge capacity is significantly improved when the diffusivities of the α and β phases are one order of magnitude higher and accordingly reduced when the diffusivities of the α and β phases are one order of magnitude lower. Furthermore, Figure 10b shows that the rate of heat generation during discharge is considerably affected by the magnitude of the α and β phase diffusivities, where the heat generation rate is higher with lower values for the α and β phase diffusivities and vice versa. Therefore, achieving efficient solid-state lithium diffusion and fast phase transition inside the LFP particles is of paramount importance for improving the rate capability of LIBs and alleviating their thermal management requirement.

Multi-dimensional distributions of potential, temperature, and current density are also obtainable from the 3-D multiscale LIB simulation. The larger LIB geometry (10.5 cm × 0.67 cm, described in Figure 11) is adopted to investigate the key multi-dimensional contours and provide greater insight into the large-scale LIB operating characteristics because the uneven distribution of key variables is more prominent in a larger scale cell. The current density contours over the separator are plotted at six different discharge elapsed times in Figure 12. The expression of the local current density, $\overrightarrow{i_e}$, is obtained based on the charge conservation equation within the electrolyte phase, Equation (11), as follows:

$$\overrightarrow{i_e}=-{\kappa }^{eff}\nabla {\varphi }_e-\frac{2RT{\kappa }^{eff}}{F}\left(1-t_{+}\right)\left(1+\frac{\partial \ln f}{\partial \ln C_{L{i}^{+}}}\right)\nabla \ln C_{L{i}^{+}}$$

(33)

The highest local current density is formed underneath the positive current collecting tab, and it decreases with increasing distance from the tab. This is because the separator region near the positive tab has the lowest internal resistance for lithium ion and electron transport, while the positive tab made of aluminum exhibits a lower electron conductivity than that of the copper negative tab. As discharging continues, the uniformity of the current density distribution slightly decreases. Figure 13 shows the electronic phase potential distributions in the positive and negative current collectors at the same discharge elapsed times. Owing to the high electronic conductivity of the current collector (3.83 × 10⁷ for the aluminum positive current collector and 6.33 × 10⁷ for the copper negative current collector), a high degree of uniformity in the potential distribution was achieved inside both negative and positive current collectors. The electronic phase potential for the positive current collector is lowered at the latter discharge stage because the boundary value of 0 V is applied to the negative tab for the solid potential equation (Equation (28)).

The temperature contours in the middle of the positive electrode are shown at the six time-instants in Figure 14. The average electrode temperature increases from 25 °C and reaches an approximate steady state temperature of 35 °C at 180 s where the heat generation rate and heat removal rate by convection are in balance. More importantly, as discussed in Figure 8 and Figure 9, the temperature profiles are similar to the current density profiles in Figure 12 because of the dominant irreversible reaction heat. Therefore, the higher temperature rise observed near the positive current collecting tab is due to the higher current density and heat generation rate. Figure 15 shows the contours of dimensionless two-phase interfacial location (r/R_p) in LFP particles over the middle of LFP electrode at six different discharge elapsed times. At t = 10 s, the r/R_p in the entire electrode region is unity, indicating that the β phase is not yet formed at the early discharging stage. However, as the discharge proceeds, the value of r/R_p decreases, implying that the β phase begins to exist from the outer wall and it’s region is continuously expanded toward the core of particle due to the lithium intercalation. In addition, the lower r/R_p is observed near underneath the positive current collecting tab where the local current density is higher and thus the rate of lithiation is faster.

4. Conclusions

In this study, a multi-scale LIB model that accounts for two computational domains at the electrode particle and cell levels, was developed and applied to LIB cell with LFP electrode. While the comprehensive 3-D cell model based on Cartesian coordinate system predicts overall species, charge, and heat transport occurring in LIBs, more detailed lithium intercalation/deintercalation and solid-state lithium diffusion phenomena inside electrode particles are rigorously taken into account by the micro-scale electrode models, i.e. defined in spherical coordinate system. Particularly, for the LFP electrode (Li_yFePO₄), the two-phase modeling approach based on a Landau transformation was employed for the positive electrode model to effectively handle the moving boundary between the α and β phases. The multi-scale LIB model was validated against the experimental discharge voltage curves measured under different C-rates, where the poor rate capability of the LIB due to sluggish lithium diffusion and phase transition inside LFP particles was successfully captured. Detailed simulation results highlight different thermal behaviors of LIBs for the charging and discharging processes. Based on the same C-rate, relatively less heat is released during discharge compared to charge, which is mainly due to endothermic contribution in the negative electrode during the discharge. In addition, the irreversible heat due to the lithium intercalation/deintercalation processes in the LFP electrode is a significant contributor to overall heat generation, particularly at high C-rates, which demands a more rigorous thermal management for fast charging and discharging.