Multiscale modeling and characterization for performance and safety of lithium-ion batteries

A key for diversifying and commercializing alternative energy sources is the development of advanced electrical energy storage technology with new materials, chemical systems, and manufacturing processes. Efficient energy storage would enable better utilization of intermittent sources such as wind and solar, increase market penetration of hybrid and all-electric vehicles, and consequently, reduce oil consumption and greenhouse gas emissions. However, particularly for vehicles, development of safe and economically viable electrical energy storage systems faces significant scientific and engineering challenges.¹ To address these challenges, we need an integrated approach that combines experiments, theory, and simulations.

In order to meet the transportation and storage needs,² electrical energy storage devices (batteries and supercapacitors) require energy and power densities that are at least a factor of 3–5 higher than those corresponding to the current state-of-the art Li-ion technology. The fundamental difference between the two is that batteries store energy by faradaic processes, generating charge in the battery electrodes, whereas the electrochemical capacitors store charge directly by forming electric double layer at the electrolyte/electrode interfaces.³ The batteries can store and release large amounts of energy, albeit slowly, therefore having comparably low power density (slow charge/discharge). They have limited lifetimes due to material degradation associated with repeated electrochemical reactions at the electrodes. Conversely, supercapacitors have modest energy density but have high power density with little loss of performance during charging/discharging cycles, since involves no chemical reactions at the electrodes. In this paper, we focus on Lithium-ion batteries (LIBs).

Most high-capacity rechargeable LIBs in mobile power applications (vehicles, tools, phones, etc.) are based on intercalation reactions in which the mobile cations are inserted into a host structure of electrodes. These batteries synchronize ion and electron transport between the cathode, electrolyte, anode, and current/heat collectors. In cases with temperature gradients during operation or short circuit, heat conduction through the layers and phonon interactions with electrons and ions will add another transport mechanism into the batteries. Any mismatch in the transport across the scales can contribute to local temperature excursions and heterogeneities of the chemical species that can result in catastrophic failure, such as explosion of this hazardous mix of fuel and oxidizer (analogous to solid propellants). A typical Li-ion cell sandwich consists of positive (cathode during discharge) and negative (anode during discharge) porous composite insertion electrodes, a porous separator, negative (copper) and positive (aluminum) current collectors. The negative electrode is typically a carbon based composite, while a lithium metal oxide (e.g., LiCoO₂, LiMn₂O₄, LiNiO₂, or mixed [Li(Ni,Co)O₂]) is used for the positive porous electrode. The amount of energy stored in a battery is directly proportional to the amount of active material that moves back and forth between the electrodes. Battery power is limited by (1) the rate of diffusion of Li cations into cathode or anode particles, (2) the electrical conductivity of cathode materials, and (3) the rate and distance of ionic transport in the electrolyte and between the electrode particles.

Repeated charging and discharging degrades a Li-ion battery, affecting its chemical and structural stabilities, capacity, voltage, power, and safety. Further improvements in the battery performance will require addressing both materials science and chemistry challenges.⁴ At the device level, scale-up and integration of LIB unit cells into battery packs for vehicle powertrains⁵ involve system power and thermal management, abuse tolerance, and safety. Temperature imbalance can lead to shortened life; and in the extreme case, temperatures in excess of 90 °C can trigger exothermic chemical reactions initiated by metastable Solid-Electrolyte Interphase (SEI) layers that can lead to cell thermal runaway.^6–8 Thermal excursion can be especially severe when a cell is unintentionally shorted or otherwise abused. To date, a wide range of cell chemistries and additives have been developed to optimize cell performance for different applications. New materials can improve a variety of performance characteristics: increase energy/power density,⁹ reduce impedance, reduce self-discharge, increase voltage, improve efficiency, reduce cost, increase cycle life, prevent corrosion and leakage, control polarization, and increase safety. At the cell and the pack levels, three-dimensional structures have been proposed for reducing the diffusion paths.^10–13 

Mechanical deformation¹⁴ of batteries can trigger thermal runaway events through short-circuit. The deformation (through abuse) can lead to rearrangement of electrode material, transfer of normal and shear stresses to the separator, and eventually the rupture of separator.¹⁵ This can on further impact lead to various types of internal shorts such as between the anode and aluminum collector, cathode and copper collector, or even aluminum and copper collectors.^16,17 The electrical contacts and the subsequent thermal runaway event will depend on the contact area, contact resistance, cell capacity, and the ability of the current flow through the short to generate sufficient heat to cause the local temperature to rise above ∼90 °C to trigger the thermal runaway events.

Given all these complex options for chemistry and materials, it is important to have a simulation capability and models that can predict performance, design, and support experiments; optimize material components and geometries; and estimate safety and durability in an integrated fashion.¹⁸ 

Figure 1 shows the computational methods available for modeling and simulation of the batteries, electrochemical systems, and materials at various length and time-scales, which can be summarized as

• electronic structure methods for modeling the atomistic structures and properties¹⁸ that provide fundamental insight into the atomistic and electronic processes governing energy and power densities, as well as structural integrity;

• classical force-field methods for the dynamics and structure of materials,^19,20 which are used to study defect formations and evolution;

• kinetic Monte Carlo (KMC) for modeling electrochemical reactions²¹ in interfacial chemistry;

• phase field method^22–24 for modeling mesoscopic transport of ions and species in order to understand the various resistances at the secondary particle level and optimize the chemistry;

• coupled micro-macroscopic models to simulate spatiotemporally varying fields, such as ions and species^25–27 at the cell/electrode levels to optimize safety and performance;

• macroscopic and/or system-level models^28–39 for simulation of full cells and battery packs, including life prediction and abuse scenarios.

Till date, battery simulations relied heavily on simplified unit cell models (e.g., 1-D and pseudo-2D) and equivalent circuit models to account for cell performance variations,^7,40,41 thereby invariably ignoring the underlying cell-to-cell electrochemical and thermal coupling.

In this paper, we first describe the experimental techniques spanning the macro- and meso-scales illustrated in Fig. 1. In particular, we describe detailed thermal and mechanical characterization techniques followed by methods that provide spatial distribution of lithium. In Section III, we describe the underlying physics models, design, structure, and capabilities of the new simulation system VIBE—the Virtual Integrated Battery Environment. VIBE was developed under the U.S. Department of Energy Computer Aided Engineering for Batteries (CAEBAT) program and is described briefly in Section IV. This integrated simulation environment provides researchers with a tool to more accurately and efficiently explore problems ranging from the effects of crystallographic texture in electrodes on cell performance to modeling impact and internal short circuits. Several applications of VIBE are described in Section V, along with experimental validation where possible.

The most common cell measurements carried out for model development and validation are (a) charge and discharge cycling at different rates,⁴² (b) cyclic voltammetry,^43,44 (c) Electrochemical Impedance Spectroscopy (EIS),^43,45 and any of the above can be complemented by (d) temperature detection at few locations using thermocouples.⁴⁶ These measurements mostly reflect the global response of the cell, module, or pack. We have developed several thermo-mechanical testing systems for measurement of mechanical and thermal responses of battery cell components (cathode, anode, and separator), battery cells, and cell stacks under normal or abuse conditions. The experimental data are used to develop and validate the cell performance and safety models. The experiments are also used to provide input to the models, such as thermal properties of Li-ion cell components.⁴⁷ Typically, the models are validated using global response such as charge/discharge profiles, and thus there is a need to obtain spatiotemporal variation of Li within the cell—that is where the characterization tools such as micro-Raman and Neutron imaging are very useful and are described in the later part of this section.

Temperature is one of the important parameters of a working battery. For normal operation with long life, the battery needs to be maintained at a uniform thermal profile, preferably close to room temperature. However, the temperature can rise not only during normal operations but also under abuse scenarios. Thus, monitoring temperature and providing thermal management is important for safe and efficient operation of batteries. In addition, because the temperature is a reflection of electrochemical activity, it is a good measure to validate coupled electrochemical-electrical-thermal models. Heat sources arising during cycling reflect polarization and charge transfer impedance, in addition to bulk resistance in electrodes, and thus, thermal measurements connect observed temperature with underlying electrochemical processes.

In our case, multiple-channel cycler with 25 A current was used to study the thermal responses of Li-ion cells under various charging/discharging conditions. Infrared (IR) imaging cameras and contact temperature sensors were used to measure the surface temperatures of cells during the charging/discharging and mechanically induced short circuit as described below.

Up to 200 A discharging current was used to study temperature responses of large format cells (1C to 5C for 40 Ah cells or cell strings). The surface temperature was measured using an IR camera, and several thermocouples were mounted on the surface.⁴⁸ The IR camera reading was calibrated using the surface thermocouples. IR camera FLIR A325 micro-bolometer focal plane array (FPA) camera was used in this work. It has 320 × 160 pixels and operates at 8–12 μm wavelength with the maximum frame rate of 60 Hz.

Figure 2 shows IR imaging during the internal short circuit test. For cells that undergo thermal runaway after short circuit, a little material is left to be examined after cell combustion. However, with IR imaging, the histories of temperature distribution, current, and deformation provide practical information set that is used for safety model development and validation.

For mechanical abuse, one needs controlled application of load on the cells to determine the displacement, stress, voltage response, temperature distribution, etc. For studying these scenarios, a mechanical load frame was built at Oak Ridge National Laboratory to apply lateral load to stacked cell components (anode, cathode, separator, etc.) or fully assembled batteries using different loading setups, including single-side indentation, double-side pinch, pinch-torsion, and high speed impact. The loading was applied through a Kollmorgen servomotor, which has a rotary actuator that allows precise control of angular position. It uses a screw to translate rotation into linear motion with the speed of 0.01 mm/s to 5 mm/s. The spherical indenter radii varied from $14$ in. to 2 in. in diameter. Below, we describe various tests developed to study cell components, internal short-circuit triggered by pinching, high-speed crush, and impact testing. Some of these tests are suitable to extract material properties that can be used as inputs to models, while others provide valuable data to validate the detailed models described in Section III D for various strain rates.

A servomotor system was used to apply compressive force to the battery components (stacks of layers) or full size Li-ion cells. A single stack of anode/separate/cathode was sealed in a pouch and partially charged. The open circuit voltage was monitored during the mechanical loading. At high data acquisition rates, a small voltage drop (e.g., 0.025 to 0.1 V compared to cell voltage of ∼4 V) due to mechanical failure (rupture) of the separator and contact of the electrodes can be detected. For cells without electrolyte (sometimes referred to as “dry” cells), a circuit consisting of an external battery (1.5 V to 9.0 V) and current-limiting resistor is used to monitor the onset of short circuit. Thermocouples and infrared camera described in Section II A are used to monitor the surface temperatures to quantify the effect of short-circuit.

An important aspect of Li-ion cell safety research is to develop an understanding of the origin of the internal short circuit and its relationship to thermal runaway. Many techniques have been proposed to study the onset and evolution of internal short circuit, such as nail penetration, small indentation, and embedded instigators.^49–51 Our method⁴⁸ uses two spheres to laterally pinch the battery from two sides. The strain rate and indenter geometry can be modified to cause single layer failure deep inside the Li-ion cells.⁴⁸ The test results are used to model the deformation and failure of the cell components by finite element analysis.⁵² 

By controlling the loading speed and using small diameter spheres (1/2 in. in diameter), small (1 mm–2 mm) punctures are created in a single layer of separator. The resulted short circuit can cause local discharge of energy and heating. The amount of energy and the temperature rise depend on many factors, such as cell total capacity, state of charge (SOC), electrode chemistry, separator, and electrolyte. Conditions that trigger thermal runaway are also studied in Ref. 48.

The polymer separators can deform quite extensively under compressive load without breaking. In our experiments,⁵³ we discovered that the separator failed faster when additional shear force was applied, so that a new torsion-compression test was developed.¹⁵ Using a bi-axial load frame, a small twisting is added during the compression of the cell. When torsion is added to the pinch test, it results in consistent and early onset of short circuit in small pouch cells.

The same pinch test system was applied to large format cells intended for electric vehicle applications. A scaled up test chamber is needed to handle the possible smoke and fire. In order to induce a minimum size internal short circuit, two small indenters (1/2 in. in diameter) were used in the pinch test and very slow loading speed (25–50 μm/s) was applied. A torsion component was added by rotating the cell in the chamber 5° after about 50% of the pinch-only indentation depth in order to induce an internal short at smaller compressive strain. There was a clear improvement in the consistency of short circuit conditions and higher capacity cells or higher state of change can be tested without thermal runaway.¹⁵ 

In order to simulate response of battery cells during a vehicle crash, a high speed loading setup was developed. The electric-driven test frame was redesigned to reach a maximum speed of 12 in./min or 5.08 mm/s. Although this is still below typical vehicle collision speeds, this improvement provides some mechanical deformation data. Most large format cells are about 6 mm thick, so that the entire test is over in about 1 s.

To achieve even higher impact speed, two tests were developed:

A large diameter ball (3 in. in diameter) is dropped from a guide tube from 2 to 6 ft. A small, stainless steel sphere is placed on top of the cell stack or single cell. The large ball transfers the impact impulse to the small sphere, which causes deformation to the stack or cell.

A high speed air gun is used to send a projectile at high speed into the cell or cell strings.

Thermal and mechanical characterization techniques described above provide a macroscopic description of the battery or battery components under different operating parameters. The macroscopic validation using that data assures that the models agree with the experiments in terms of an average response, but still there is a need to validate the results at smaller length scales to be more dependable. In this section, we provide a brief description of the spectroscopic and imaging methods applied to battery electrodes for studying local state of lithiation or local lithium content at the micron scale. The method is based on correlation between the observed spectroscopic or imaging signal at a given SOC or voltage of any given battery electrodes and the local inhomogeneity at a relatively large length scale between tens of micrometers to millimeters. The measurements on cycled electrodes are compared to pristine (uncycled) electrodes to look for possible electrode degradation. The degradation can be the result of: (i) loss of structural integrity of the electrode materials, (ii) unwanted electrode side reaction at specific sites or regions of the electrodes, and (iii) electrochemically inactive regions. All these factors contribute to capacity degradation when the battery is continually cycled.

Raman microscopy is a tool for investigating battery cathode and anode materials and their evolution under electrochemical cycling.^54,55 The Raman spectra of battery chemistry show changes against lithiation-delithiation and can be used as spectroscopic fingerprint to observe variation of local structural changes under continuous electrochemical transport. Nanda et al.⁵⁶ have analyzed the SOC dependence of LiNiCoAlO₂ (NCA) electrodes and mapped the regions of degraded region of cycled electrodes that were subjected to automotive duty cycles. A similar method has been applied to the high voltage lithium-rich cathode material, Li_1+xMn_{1−x−y−z}Ni_yCo_zO₂, where x is nominally in the range of 0.2.⁵⁷ For anodes, where carbon is still the most common material for lithium-ion applications, Raman spectroscopy provides important insights into degradation. For pristine or un-cycled anodes, the Raman signal from carbon primarily has two main bands, commonly referred to as D and G bands. The G band originates from 2 degenerate E_2g modes that correspond to in-plane vibrations of carbon atoms. The D1and D2 bands appear in disordered graphite samples, where the local bonding deviates from ideal sp² geometry. The extent of disorder can be described by the ratio of their areal intensities, given by A_D1/A_G. When the electrodes are cycled continuously under a duty or load cycle, the structure or order of the graphitic regions could be compromised by repeated intercalation and de-intercalation. Figure 3 demonstrates the changes observed in a cycled carbon electrode. We notice non-uniformity in the carbon regions in terms of their degree of disorder as shown by the change in the color contrast over the electrode spatial region. The brighter colors as shown in the color bar in Fig. 4 indicate higher degree of disorder (larger D/G ratio). The region indicated in blue as shown by the corresponding Raman spectrum is highly graphitic followed by green and yellow areas. The abundance of yellow region does indicate possible exfoliation of graphitic domains under extensive electrochemical cycling. The region in black does not correspond to carbon and could be either electrolyte salts or organic residues. Raman microscopy is also a powerful technique for probing materials phase evolution for non-intercalation type chemistries such as alloying reactions for high capacity silicon anodes.⁵⁴ 

Raman imaging (microscopy) discussed in Section II C 1 has a spatial resolution between 0.5 and 1 μm and probes between 50 and 100 nm in the bulk depending on the skin depth (conductivity) of the materials. In our goal towards probing the bulk microstructure of the electrodes in 3D, we have recently applied neutron imaging methods for battery electrodes. Given the sensitivity and contrast change between Li⁶ versus Li⁷, neutron imaging can be a powerful tool for resolving the 3D microstructure of the electrodes during electrochemical transport. Imaging methods combined with isotopic substitution especially varying the Li⁶/Li⁷ ratio and deuterating the electrolyte solvents to prevent strong absorption from hydrogen could further increase the contrast. We recently applied neutron imaging technique on discharged lithium-air electrodes to monitor the spatial lithium distribution across the air-cathode thickness.⁵⁸ The current spatial resolution of this method is on the order of 40 μm. Our results show that for thicker air cathode, we notice anisotropic distribution of discharge products (mainly lithium salts decomposed during the discharge process). 3D electrochemical modeling simulated using actual experimental parameter predicts kinetic limitations due to both lithium transport across the thick cathode (1 mm) and limited diffusion of oxygen at the other end of the cathode. This combined effect leads to spike in lithium concentration at the edges of the cathode, and the lithium concentration drops towards the center.

In Section II, we described experiments that provide thermo-electrochemistry-electrical and mechanical response. In addition, we have techniques that can probe Li distribution as a function of space and time at the mesoscale. In this section, we describe models that can simulate the phenomena measured by these experimental techniques and represent the four major physical phenomena in Li-ion batteries: (a) electro-chemistry, (b) electronic transport, (c) heat transport, and (d) mechanics. The electro-chemistry model acting over the electro-active components provides heat source and electrical resistance but requires temperature, geometry, stress, current, etc., from other components. The electrical component acting over the entire device needs resistance and geometry but provides current, heat source, etc. The thermal component acting over the entire domain needs heat sources, geometry, but provides temperature. The mechanical component takes temperature and provided deformed geometry/stress etc. If the processes are separated in spatiotemporal scales, one can perform loose coupling and obtain accurate solutions. However, in reality, the scales overlap, and one needs the ability to couple these models with various degrees of coupling to obtain accurate simulation. The software developed to enable such integration of components with different coupling is described in Section IV.

Several models applicable to redox-based energy storage systems have been developed with different degrees of complexity and idealization. In the most general form, the mass and charge conservation in 3D can be described with four differential equations applicable to host material (subscript s) and electrolyte (subscript e) phases^25,59

where $cs$ and $ce$ are solid and electrolyte lithium concentrations, $φs$ and $φe$ are solid and electrolyte potentials, F is the Faraday constant, and $t+0$ is the cation transference number in the electrolyte. Conservation of species (Eq. (1)) and conservation of charge (Eq. (2)) are connected via flux of charge associated with lithium ions $jLi$⁠. The latter is described by the Butler-Volmer equation of cathodic (anodic) reaction at the interface

which is driven by the local overpotential $η=φs−φe−U$⁠, with $U$ being the open circuit potential (OCP) of a particular half-cell reaction under consideration. The exchange current density $i0$ is a function of lithium concentration

where $cs,max$ is the maximum stoichiometric content of Li in the host electrode material.

Effective solid and electrolyte transport properties, i.e., diffusivities (⁠$Dseff$ and $Deeff$ correspondingly) and conductivities (⁠$σeff$ and $keff$⁠), are described in terms of corresponding volume fractions $εs$ and $εe$ as exponential relationships of a general form $Θeff=εpΘ$⁠. Exponent p could be taken as Bruggeman coefficient of 1.5 or used as an adjustable parameter to fit the discharge curve data.⁶⁰ The diffusional conductivity $κDeff$ is determined by using the concentrated solution theory⁶¹ as

where $f±$ is the mean molar activity coefficient of the electrolyte.

Equations (1) and (2) describe the mass and charge balances in electrode and electrolyte domains of a redox cell in three dimensions. Several simplifications to this general model have been proposed. One of the more well-known approaches is the pseudo-2D (P2D) model developed by Doyle, Fuller, and Newman.^60–62 The model, sometimes referred to as the DFN model, is based on porous electrode theory and casts the equations of diffusion and charge transfer onto a simplified geometry. This geometry is described by two coordinates—one through the electrode thickness and the other being radial coordinate of an electrode particle idealized as a sphere. This additional coordinate related to the solid state diffusion in particle gives the model its name “pseudo 2D” because it is a one-dimensional approach where transport of species occurs along the thickness direction. Transport through the electrolyte is modeled by using the concentrated solution theory, and the pore-wall flux of lithium ions across the electrode-electrolyte interface is naturally set to zero in the separator region of the cell sandwich, while at the active material interface, it follows the Butler-Volmer equation. With such a simplification, solid state diffusion reduces to a 1D equation in radial coordinate r of a particle with radius $Rs$

In order to reduce compute time, Eq. (6) can be solved either by Duhamel superposition method²⁸ or by using approximations for the surface lithium concentration that come from the diffusion length formulas derived by Wang et al.²⁵ In the latter case, the approximation for surface concentration can be expressed as⁶² 

where $ls=Rs/5$ is the diffusion length corresponding to a spherical geometry with radius $Rs$⁠. With such a description, the initial surface concentration is equal to the average concentration of lithium in the solid and approaches a linear asymptote at long discharge/charge times of $cavg+jLilsFDs$⁠.

The second category of electrochemical models is based on a description of cell behavior from experimentally measured impedance characteristics. The model was developed by Newman, Tiedemann,⁶³ and Gu⁴² and accordingly is commonly referred as the NTG model. The cell current density (J) is linearly related to the cell potential (V) as

where Y and U represent the effective conductance and the OCP of the cell, respectively. They are expressed as polynomial functions of state of charge variable $θ$

The degree of polynomial (N) and the values of constants $ai$ and $bi$ are determined from the cell discharge curves for a number of C-rates. When the cell potentials are plotted versus applied current density, the cell voltages at zero current represent the OCPs of the cell, and the slopes of potentials represent the cell impedance $1/Y$⁠. Linear dependence of the latter on the current density is assumed, resulting in Y being a function of state of charge only. While the NTG model allows rapid simulation of large systems (modules or packs), it provides only the overall cell response and hence cannot be used in situations where local variations in species or effects are of interest.

We are also developing a 3D electrochemistry simulation capability known as AMPERES^8,64 (Advanced MultiPhysics for Electrochemical and Renewable Energy Storage) that preserves the full 3D treatment of Equations (1) and (2) and couples electrochemistry with electrical and thermal models on the same grid. This level of integration might not be critical for simulating all performance scenarios but can be significant in situations such as abuse, where there is less scale separation between processes. Detailed study of this model is beyond the scope of the present article.

The conservation of charge within the cell (or module or pack) is expressed by the Laplace equation

over the entire union of cell sandwich, current collector, and electrical interconnect component domains. Appropriate boundary conditions corresponding to charge or discharge behavior are applied to the current collector tabs. The values of conductivity $σ$ of the cell sandwich are determined as a result of the solution of the corresponding electrochemical set of equations.

In its general form, the three-dimensional heat conduction equation is expressed as

where $ρ$ is the density, $Cp$ is the specific heat capacity, $k={kx,ky,kz}$ is the anisotropic thermal conductivity, and T is local temperature. In Li-ion battery cells, the anisotropy of thermal properties is a result of their layered cell structure. The cells' in-plane conductivity can be several times higher than out-of-plane thermal conductivity. Equation (11) is solved for the entire cell, module, and pack. Computation of the heat source, q, depends on the domain. In electrical interconnect components and current collectors, heat source is represented by a simple ohmic heating model. The heat generation model within the electrochemical dual-electrode insertion cells was derived by Bernardi et al.⁶⁵ The heating within the cell, in addition to electrical work, is attributed to electrochemical reactions, changes in the heat capacity of the system, phase changes, and heat of mixing. In its simplified form, which commonly employed in battery simulations, the heat generation consists of irreversible energy loss due to cell polarization, reversible entropy change due to particular half-cell reaction and ohmic heating within cell sandwich

where the summation in general occurs over all reactions $j=1…M$⁠, which in case of Li-ion intercalation systems simplifies to two half-cell reactions at each electrode with corresponding OCP vs. lithium being denoted as $Uj$⁠. When the NTG model is used, the set of model equations is imposed on the cell level and the local solid and electrolyte potentials are not known. Thus, in this case, Eq. (12) is applied in its integral form as

where $η$ is the cell overpotential $(Vp−Vn−U)$ and h is the cell sandwich thickness.

Considerations of mechanical stress-strain state have entered the battery simulations relatively recently. Generally, in battery modeling, the term “mechanics” can refer to (a) cell (module or pack) response to externally applied loads and (b) stress-strain state inside the cell induced by internal, “electrochemical” loading. The former case is typically considered as an instantaneous structural response of a battery having some effective response to external load.^66,67 In the latter case, deformation on the microscopic level is induced by solid state diffusion and phase boundary transition and thus is a time and rate dependent process (although the mechanical deformation is still considered as instantaneous at each time step). In this work, we are primarily concerned with mechanical response of the batteries to external loading, which is one of the major safety concerns. Large deformations of the battery pack that may arise due to crash loads can create internal short circuit that may result in thermal runaway.

Mechanical behavior of the foil materials, being that of copper and aluminum, can be described by the elastic-plastic formulation

In the above equations, the total strain is the sum of elastic and plastic strains $ε=εe+εp$ and $dp=dε¯¯p:dε¯¯p$ is the increment of equivalent plastic strain. Yield surface $f=0$ is described in terms of deviatoric part (⁠$S¯¯$⁠) of the stress tensor $σ¯¯$ with the backstress $α¯¯$ either following Armstrong-Frederick type hardening rule or set to zero for pure isotropic hardening. These relationships are well established for ductile metallic materials that undergo crystallographic slip and consequent plastic flow.

The mechanical behavior of polymer separators can be described by the viscoelastic equations typical for polymers. It should be noted, however, that separators have a porous microstructure and the properties of dry and wet separators could be significantly different.⁶⁸ In order to accurately include influence of liquid-filled pores on the macroscopic mechanical behavior of battery separators, poroelastic theory can be utilized. In this case, the effective stress includes the influence of the liquid pressure inside the pores as

where p is the pore pressure and b is the Biot's coefficient related to the ratio of moduli of porous and bulk forms of separator polymer

With this description, the constitutive equations of elasticity (first equation in Eq. (14)) are represented as

where E and $ν$ are elastic properties of dry separator.

Mechanical behavior of composite coatings on the current collectors is not well understood. Electrode coatings contain phases of very different properties. Active materials in cathodes are metal oxide particles (which behave like ceramics) connected by flexible binder material into an electronically conductive network. Exact elastic properties of the constituents are not known at the present. However, effective properties of the composite electrode layers can be estimated by experiments on cells and inverse modeling to match the experimental response.

The structure of developed computational simulation suite, named VIBE, is illustrated in Fig. 4.

The Open Architecture Software (OAS)^69–71 is a modular and extensible software infrastructure that supports multiple modeling formulations and computer codes. The computational models use the OAS to couple the physical phenomena. The coupling strategies can range from one-way coupling, where the results from one physics model are provided to another, to fully-implicit nonlinear coupling where all the physics models are solved simultaneously. The design of the environment allows integration of new models to replace or enhance existing components. The battery state stores the variables describing the state of the battery and is updated by each modeling component as the simulation advances in time from one state to the next. A user interacts with the system either by command line or through the Eclipse (http://www.eclipse.org/) based Integrated Computational Environment (ICE).^72,73 ICE (http://www.eclipse.org/ice) is a general purpose set of interactive graphical tools for four primary activities that are common across all modeling and simulation projects: creating input files, launching jobs, visualizing and analyzing simulation results, and managing data.

As an example of coupled simulation, we modeled performance of a 4.3 Ah pouch cell having dimensions 70 mm × 5 mm × 110 mm. The simulation included electrochemical, electrical, and thermal transport components launched in this sequence by the simulation driver at each time step. The electrochemical processes in the cell were modeled with NTG model using the experimentally measured discharge profiles at different C-rates. The details of the cell description and numerical values for the NTG parameter fits for impedance and OCP can be found in Ref. 69.

The cell model included all 17 anode and cathode electrodes, and the nylon/aluminum/polypropylene pouch. The latter was included in thermal simulation step only. An idealized weld with zero resistance connected the current collectors to the cell tabs where a current flux corresponding to the applied C-rate was imposed as a boundary condition. Cooling of the cell surfaces was simulated by convective boundary conditions with a heat transfer coefficient that was determined from experiments as 15 W/m² K.⁶⁹ This value was obtained experimentally from the surface measurements of decreasing temperatures of cells during post-discharge cool-down in ambient laboratory air in the absence of heat sources. Figure 5(a) shows the computed distribution of temperature in the cell at the end of 5C discharge. The results of temperature evolution with time for different applied currents were validated against the experimental IR measurements described earlier. As can be seen from comparison of experimental and simulated temperature profiles in Fig. 5(c), the modeling approach is capable of predicting the cell behavior with good accuracy. This is rather promising considering the compute time savings provided by NTG model for electrochemical behavior.

Four prismatic cells with validated model parameters (Fig. 5) were used to assemble battery modules. The cells were connected either in parallel (4P) or in series (4S). Each of the cells contained 34 charge transfer zones (regions over which the different physics are averaged and coupled) and 37 current collector zones (including tabs). The overall module system thus consisted of 308 zones for connection in parallel and 300 zones for connection in series. The resulting finite element model (FEM) has approximately 150 000 nodes. Temperature distribution in the module under uneven cooling and high discharge current was computed. The two surfaces of the module were cooled at different rates. One side was cooled imitating slow natural convection at room temperature with heat transfer coefficient of 15 W/m² K. The other side had a forced cooling applied with the same coefficient equal to 55 W/m² K, which imitates cooling with fast moving air. No cooling was applied to the tabs and interconnects. 5C discharge current was applied in both cases. The results are shown in Fig. 6. Similar temperature distributions can be observed in both 4P and 4S configurations. In the current scenario of an air-cooled module, uneven cooling can be considered a result of failure of the cooling system, which stops forced flow of air on one side of the module (air duct blockage or fan malfunction). This results in significant temperature gradients across the module.

It can also be seen from Fig. 6 that when only the sides (i.e., the pouch material) of the module (or a cell) are cooled, a high temperature occurs at the tabs because a significant amount of heat is transported via metal current collectors. Dissipation of heat through metal current collectors can be used for cooling strategies alternative to the standard practices of using cooling fins sandwiched in between every 3rd and 4th cell.

In order to demonstrate this concept, a thick cell containing 84 cell sandwiches across the thickness was modeled. The same electrochemical model as in case of prismatic cell discussed above was applied with the discharge current scaled correspondingly to produce a 5C discharge rate. Figure 7 shows the comparison of the coupled electro-chemical and thermal solution under various cooling boundary conditions. Idealized liquid cooling was simulated by holding the corresponding cell surfaces at 295 K (Dirichlet boundary conditions). The rest of the surfaces had convective air cooling applied with convective heat transfer coefficient equal to 15 W/m² K. As can be seen, the cell with the new design allowing core cooling via current collectors operates at controlled temperature range with improved uniformity (Fig. 7(c)).

The benefits of cell liquid cooling via current collectors are evident in the through-thickness temperature profiles in Fig. 8. The figure shows the same cooling scenarios as in Fig. 7. The lowest temperature across the cell is achieved with liquid cooling applied to current collectors (Fig. 8(c)). It is evident from these plots that the cooling of the battery from sides by extending the current collectors not only allows us to manufacture thicker cells but also reduces the thermal gradient by accessing the core of the battery.

The above examples demonstrate robustness of loosely coupled electrochemical-thermal simulations in which underlying cell performance is described by the NTG model. The latter is assumed independent of temperature, which may be applicable only within a limited range of temperatures and cell currents (C-rates). More detailed description of the underlying physics, based on the porous electrode theory, provides lithium ion pore-wall flux, solid diffusivity, and diffusional conductivity of electrolyte as functions of temperature (Eqs. (3), (5), and (6)), thus enabling posing a two-way tightly coupled problem. To investigate the effect of tight coupling between thermal and electrochemical models, we applied a fixed-point iteration scheme at each time step (Picard iteration). In addition to exchange of variables (heat sources and temperature) via the battery state, such iteration demands convergence with pre-selected tolerance. In our case, the difference in temperature within one iteration should not exceed 0.1 K. A simple case of unrolled cell sandwich strip cooled at one end and consisting of spinel LiMn₂O₄ positive electrode paired with carbon anode was chosen.²⁸ The solution is shown in Fig. 9 for the applied current densities of 17.5 A/m² and 52.5 A/m² (1C and 3C rates correspondingly). Heat sources determined by Eq. (12) are shown as functions of time in Fig. 9 together with the thermal solution corresponding to 52.5 A/m² discharge. As can be seen, introduction of tight coupling between thermal and electrochemical components influences heat sources only at high applied currents. Apparently, the heat sources in electrodes attain slightly lower values when the appropriate temperature feedback is introduced into electrochemical component. Higher temperatures apparently increase diffusional conductivity of electrolyte, thus decreasing concentration polarization. Within the pseudo-2D model, solid state diffusion in a spherical particle (Eq. (6)) is solved by using the Duhamel superposition approximation, which is applicable when the diffusivity is constant. Thus, the diffusion constants $Ds$ in cathode and anode are not functions of temperature in this solution, unlike electrolyte conductivity and diffusion as well as exchange current density at the positive and negative electrodes, which all have Arrhenius type temperature dependencies.

As can be seen, at normal cell performance, introduction of tight coupling does not influence the solution significantly compared to the case when exchange of variables between thermal and electrochemical components occurs at the end of time step. Due to such low dependence, the convergence of Picard iteration is rather fast and in most cases occurs within four iterations. Therefore, the time penalty in order to obtain more accurate solution is negligible. In order to achieve better description of underlying physics, however, the solid state diffusivity should be dependent on temperature with corresponding activation energy for lithium diffusion determined from experiments.⁷⁴ 

Safe performance of Li-ion cells, modules, and battery packs is of a particular concern for electric vehicle manufacturers. The volatile components of a battery present a significant fire hazard in the event of crash. If an excessive compressive load is imposed on the cell, it can lead to a failure of the separator—the protective barrier that prevents the positive and negative electrodes from contacting each other and causing short circuits. Once a short occurs, the local rise in temperature may lead to thermal runaway. Therefore, from the safety modeling perspective, the ability of a model to detect the location and geometry of a short circuit contact in response to external load is essential. In order to be predictive, a model must resolve the components of the cell—separator, electrode coatings, and current collectors. At present, however, most of the structural mechanics models of Li-ion batteries are based on averaged properties that have to be determined experimentally for each type of cell.⁷⁵ In what follows, we describe the simulation of deformation in each of the cell layers for the pouch cell (180 × 150 × 6.5 mm) pinch test.^15,48,52 The pinch test involves indentation of the cell with steel sphere, 1 in. in diameter. The area of the cell that undergoes the deformation is small compared to the overall cell size so that the model depicts only the vicinity of the indentation. The corresponding finite element mesh is shown in Fig. 10, with the resolved layers that form the cell following the mechanics model as arranged in Table I.

We have performed the computational simulation in LS-DYNA⁷⁶ with fully integrated solid elements. The lower surface of the cell rested on a rigid wall. The hemispherical punch was modeled as a contact entity with rigid material. The failure of separator in this preliminary simulation was prescribed based on the maximum strain criterion with the failure strain equals to 10%. The ability of the model to resolve each component of the cell is demonstrated in Fig. 11, where the distribution of Von Mises stress on deformed mesh is shown in different cell components. The application of failure criterion in separator allows detection of the onset of short circuit. This is illustrated in Fig. 12(a), where the vertical coordinate (z-coordinate) of the cathode and anode is shown as a function of time. It can be seen that once the failure strain is exceeded in separator, the two layers come in contact, thus initiating the short circuit. This event is evidenced by the load drop in load-displacement curve in Fig. 12(b). The failure location is strictly a function of the maximum strain criterion, and in future studies, we will perform sensitivity analysis with respect to this parameter and validate against the ongoing experiments described earlier.

Predictive modeling of lithium-ion batteries requires accurate simulation of coupled nonlinear physicochemical phenomena across a range of spatiotemporal scales. We have developed an open-source computational environment for battery simulation known as the VIBE. VIBE currently includes traditional electrochemical models such as Newman-Tiedemann-Gu and Doyle-Fuller-Newman coupled to electrical and thermal models, as well as coupled volume-averaged fully 3D models. It is designed to allow new models, alternative chemistries, and extensions to microscopic/atomistic scales. Existing models have been validated with experimental data, and high-fidelity simulations have been performed in order to explore innovative approaches to cooling. VIBE is now being used in conjunction with mechanics simulations to better understand mechanisms of mechanical abuse that can lead to internal short circuits and thermal runaway. An integrated approach, combining theory, modeling, and experiments, is being used to better understand both the normal operation of LIBs and behavior under adverse conditions.

The research was performed using the resources at Oak Ridge National Laboratory (ORNL), managed by UT-Battelle, LLC, for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. The authors acknowledge the support of the Vehicle Technologies Program in the Office of Energy Efficiency and Renewable Energy, the Advanced Research Projects Agency-Energy (ARPA-E), the National Highway Transportation Safety Agency (NHTSA) of the U.S. Department of Transportation, and the ORNL Laboratory Directed Research and Development (LDRD) program.