Boosting Rechargeable Batteries R&D by Multiscale Modeling: Myth or Reality?

Abstract

This review addresses concepts, approaches, tools, and outcomes of multiscale modeling used to design and optimize the current and next generation rechargeable battery cells. Different kinds of multiscale models are discussed and demystified with a particular emphasis on methodological aspects. The outcome is compared both to results of other modeling strategies as well as to the vast pool of experimental data available. Finally, the main challenges remaining and future developments are discussed.

Figure 1

Schematic representation of a LIB cell.

Figure 2

Multiscale character of a battery cell.

Figure 3

Workflows of (a) stand-alone models, (b) multiscale models based on sequential linking (MSMSL), (c) multiscale models based on iterative coupling (MSMIC), and (d) multiscale models based on tight coupling (MSMTC). “PE” refers to “physical equation” (mathematical equation based on a fundamental physics theory which defines the relations between physics quantities of an entity) and “MR” to material relation (materials specific equation providing a value for a parameter in the physics equation).

Figure 4

A variety of kinetic phenomena, including Li diffusion and first-order phase transformations involving nucleation and interface migration, occur within individual electrode particles during each charge and discharge cycle of a Li-ion battery. Reproduced from ref (93). Copyright 2013 American Chemical Society.

Figure 5

Computation approach integrated by first-principles calculation, cluster expansion, and MC simulation. Reproduced from ref (101). Copyright from IOP Publishing.

Figure 6

Important crystal structures and Li hop mechanisms in common intercalation compounds. Many intercalation compound chemistries have either (a) a layered crystal structure (with an ABAB or ABC stacking of a close-packed anion sublattice) or (b) a spinel crystal structure characterized by a three-dimensional interstitial network for Li ions. (c and d) Diffusion in these crystal structures is often mediated by vacancy clusters (divacancies in the layered form and triple and divacancies in the spinel form) if Li occupies octahedral sites. Reproduced from ref (93). Copyright 2013 American Chemical Society.

Figure 7

(i) Schematic diagram of phase transition during deintercalation from (a) O3-LiCoO₂ and (b) O2-LiCoO₂. Blue octahedrons represent CoO₆, red balls indicate oxide ions, and yellow balls represent lithium ions. (ii) Ordered (a) O6-Li_1/3CoO₂ and (b, c) O2-Li_1/4CoO₂ phases found to be stable at room temperature by first-principles calculations. The lattice denotes the lithium sites within a Li plane, and the filled circles correspond to Li ions. For the (b, c) O2 host, the unfilled small circle denotes the projection of the Li sites of an adjacent Li plane. For the (a) O6 host, the lithium positions of the two adjacent layers are different: the projection of one is represented by the small unfilled circle, and the projection of the other is represented by the large unfilled circle. Reproduced from ref (104). Copyright 2003 American Chemical Society. Reproduced with permission from ref (105). Copyright 2012 Royal Society of Chemistry.

Figure 8

Flowchart of the theoretical approaches used in ref (114). Reproduced from ref (114). Copyright 2012 American Chemical Society.

Figure 9

Time snapshots of lithium ion composition during the lithiation–delithiation cycles for the cubic spinel Li_1+xTi₂O₄ electrode particles. Shown are the composition fields at time t = 0.3 ms, t = 0.45 ms, t = 0.6 ms, and t = 0.9 ms for the problems with no nuclei, two nuclei, and five nuclei. The blue regions with composition close to x = 0 are in the α-phase, and the red regions with composition close to x = 1 are in the β-phase. Intermediate values of the composition correspond to the two-phase regions. Reproduced with permission from ref (119). Copyright 2016 American Chemical Society.

Figure 10

(a) Trajectories of the three Li atoms originally situated in the crystalline LiFePO₄ structure along the b axis obtained from the MD simulation at 2000 K. The three specific Li atoms are highlighted with different colors, and all other Li atoms are omitted for clarity. A zigzag diffusion pathway can be clearly identified (shown with the curved arrows) along the b direction and confined in the (ab) plane. (b) Snapshots from the MD simulation of the fully lithiated LiFePO₄ at 2000 K showing the second diffusion mechanism involving the formation of Li_Fe antisite. Large (light green) and small (brown) spheres represent Li and Fe atoms, respectively. A pair of Li (pink) and Fe (blue) atoms are highlighted to illustrate this mechanism. Dashed lines indicate the size of the used supercell. Reproduced with permission from ref (123). Copyright 2011 American Chemical Society.

Figure 11

MC simulation of the galvanostatic discharge process performed at room temperature for the Li_xFePO₄ olivine nanocrystals. The gray points represent the lithium atoms in the active particle. The superficial current density is 0.5 Å m^–2. The far-field flux of Li ions is perpendicular to the b direction of the cell. Microstructure at (a) 0 s (initial random solid solution), (b) 0.0001 s (two Li-enriched phases have formed and joined together), (c) 0.00056 s (growth of the Li-enriched phase), and (d) 10.81 s (the Li-poor phase is almost consumed). (e) Cell voltage as a function of Li concentration (mol) in the active material. Reproduced with permission from ref (126). Copyright 2011 John Wiley and Sons.

Figure 12

(Left) Schematic of the cell as modeled by Newman et al. Cell consists of a Li-metal negative and a LiFePO₄ positive with a separator between them. The cell is filled with electrolyte. The LiFePO₄ porous electrode is attached to a carbon-coated aluminum current collector. The electrode is assumed to consist of spherical particles at the surface of which an electrochemical reaction occurs. (Right) Illustration of the shrinking-core model with the juxtaposition of the two phases and the movement of the phase boundary. Both the single-phase and the two-phase regions are shown. Reproduced with permission from ref (129). Copyright 2004 The Electrochemical Society.

Figure 13

(a−h) Simulations of reaction-limited phase separation of a 500 nm single crystal of Li_0.5FePO₄ into Li-rich (black) and Li-poor (white) phases in an electrolyte bath at zero current and zero pressure, consistent with ex situ experiments. (a) Coherent phase separation, where lithium insertion causes contraction along the (c) [001] axis and expansion along the (a) [100] axis, leading to tilted interfaces aligned with {101} planes. (b) Loss of [001] coherency (e.g., due to microcracks) causes the phase boundaries to rotate to align with the [100] planes. Reproduced from ref (130). Copyright 2012 American Chemical Society.

Figure 14

(top) Faradaic reaction O + ne^– → R in concentrated solutions. Each state explores a landscape of excess chemical potential μ^ex. Charge transfer occurs where the curves overlap, or just below, by quantum tunneling (dashed curves). (bottom) Example of ion intercalation into a solid electrode from a liquid electrolyte. Reproduced from ref (131). Copyright 2013 American Chemical Society.

Figure 15

Li-graphite system. (a) In-plane Li ordering in fully lithiated stage I and stage II (left) and (b) 2 × 2 Li ordering [found in stage II′ (right) seen from above on the carbon honeycomb lattice]. (c) First-principles phase diagram obtained from MC simulations based on stage I and stage II separate cluster expansions. The phase regions are denoted in the following way: G graphite, II (stage II), IID (disordered stage II), I (stage I), ID (disordered stage I), and II′ (stage II with 2 × 2 Li ordering). Different symbols refer to cooling and heating runs, respectively. Reproduced with permission from ref (139). Copyright 2010 American Physical Society.

Figure 16

(a) Density plot of Li (green) atoms in the simulation box during an NVT simulation at 1000 K obtained as an average over a time period of 250 ps in an OLC. The green color represents all grids (regions) that registered a Li count during the simulation. (b) Free energy landscape for Li diffusing out from the OLC through a 12-membered ring. Reproduced from ref (142). Copyright 2015 American Chemical Society.

Figure 17

C=A–A=C, C=A–B=C and CH–X ∼ X–CH systems: (a) High- and (b) low-potential molecular systems (with A/B equal to N or CH and X equal to NH/O (X–X bond), CH (X=X bond) or C (X≡X bond)). The spin density distribution and Bader’s atomic spin density populations (in hundredths of electrons) of two representative systems are shown to highlight their difference (a) A=CH and (b) X=CH. (c) Calculated reduction potential (in V vs Li⁺/Li) as a function of the capacity for a one-electron process. Reproduced with permission from ref (158). Copyright 2015 Royal Society of Chemistry.

Figure 18

Absolute oxidation potentials (left axis) and estimated voltages (right axes, blue for Li and red for Na) for (●) PANI and (◆) CN-functionalized PANI at different degrees of oxidation. ○ indicate the potentials obtained from pristine PANI at 24% oxidation with different spatial positions of ClO₄⁻ anions. Green △ indicate the potential obtained from PANI at 24% oxidation using a larger unit cell.

Figure 19

Comparison of the Li⁺ trajectory in pure EC and the 3:7 EC:EMC mixture showing stronger caging effects in the mixture. Reproduced from ref (178). Copyright 2017 American Chemical Society.

Figure 20

Representation of polymer chains using a coarse-grained model in 3 different perspectives. Reproduced with permission from ref (283). Copyright 2014 The Royal Society of Chemistry.

Figure 21

Schematics of simulated polymers with uncharged polymer beads in green, charged polymer beads in blue, and counterions in red: (a) and (b) periodic pendants, (c) random block pendants, (d) periodic ionenes, and (e) fully random pendants. Reproduced from ref (286). Copyright 2012 American Chemical Society.

Figure 22

Schematic of charge conduction through ionic aggregates. Blue-filled circles are cations and red/black circles are sulfonate groups. Selected cations involved in the movement are highlighted with green, purple, and orange circles. Reproduced with permission from ref (288). Copyright 2016 The Royal Society of Chemistry.

Figure 23

Backbone polymer (gray) in the atomistic chain is completely removed through coarse-graining, leaving only anions (pink) bonded together with weak harmonic springs, neutralized by sodium cations (cyan). Reproduced with permission from ref (288). Copyright 2016 The Royal Society of Chemistry.

Figure 24

Schematic illustration of coarse-grained model for salt-doped diblock copolymers. The density fields as required by the Hamiltonian are constructed by mapping bead positions onto an underlying grid. The mass and charge of beads are mapped onto the eight nearest-neighbor grid sites (a, b, c, d, etc.) and the contributions of different sites to the free energy are determined by the relative distance of each bead to a nearby grid site. Reproduced from ref (291). Copyright 2016 American Chemical Society.

Figure 25

Timescales of different events: τ₁ is the time for intrachain ionic motion, τ₂ is the relaxation time of the polymer chain, and τ₃ is the waiting time of an ion between two interchain jumps. Reproduced with permission from ref (294). Copyright 2007 American Physical Society.

Figure 26

Macroscopic Li⁺ diffusion coefficients originating from various contributions to the Li⁺ transport and the total for EO:Li = 20 concentration. Reproduced from ref (258). Copyright 2006 American Chemical Society.

Figure 27

Alternative strategies for obtaining ion transport properties in polymers. The dashed arrow indicates the conventional brute-force approach, in which transport properties are obtained by running long and computationally expensive MD trajectories. The black arrows indicate the approach of the CS-DBP model, in which short MD trajectories are used to obtain parameters for kMC simulations that predict transport properties at a reduced computational cost. Reproduced from ref (298). Copyright 2015 American Chemical Society.

Figure 28

Classification of the types of adsorption orientations for EC/PC. Reproduced with permission from ref (306). Copyright 2017 Elsevier Ltd.

Figure 29

Decomposition of EC. Initially one electron is transferred from the surface to (a) Li-EC, causing a (b) C_carbonylO_ring bond to break. Subsequently, a second electron is transferred to (c) the EC– radical anion, triggering the breaking of a second C_carbonyl–O_ring bond and (d) generating the C₂H₄ + CO₃^2– pair. The net charges of the EC molecule and the CO₃^2–/C₂H₄ products are shown. Reproduced from ref (309). Copyright 2013 American Chemical Society.

Figure 30

SEI after 16 ps of simulation using 4 M (a) LiFSI or (b) LiTFSI in DME electrolytes vs a lithium surface. Reproduced from ref (315). Copyright 2017 American Chemical Society.

Figure 31

Energy profile for breaking C₂H₄ from o-EC^–/Li⁺, as obtained from QC calculations and gas phase simulations using ReaxFF at 10 K. Reproduced from ref (320). Copyright 2012 American Chemical Society.

Figure 32

(a) Schematic of a Li metal electrode dipped in an electrolyte. (b) Initial configuration of the cell. Reproduced with permission from ref (324). Copyright 2011 Elsevier Ltd.

Figure 33

Distribution of the SEI components for different electrolytes (left). Atomic configurations from MD simulations; the components of the SEI are identified (right). Reproduced with permission from ref (324). Copyright 2011 Elsevier Ltd.

Figure 34

Potential energy profile for the reduction of EC/Li⁺ and the radical termination reactions according to various pathways. ΔE_R and ΔE_B denote the reaction energy and reaction barrier, respectively. Color scheme: cyan, carbon; white, hydrogen; red, oxygen; purple, Li⁺; large blue sphere, electron. Reproduced from ref (321). Copyright 2016 American Chemical Society.

Figure 35

Li⁺ cation diffusion coefficients for ordered and amorphous SEI models consisting of Li₂BDC and Li₂EDC. Reproduced from ref (335). Copyright 2017 American Chemical Society.

Figure 36

Snapshots highlighting the Li⁺ distribution in ordered Li₂BDC and Li₂EDC, and amorphous Li₂BDC, all at 393 K. Reproduced from ref (335). Copyright 2017 American Chemical Society.

Figure 37

Density distributions for an electrochemical cell composed of an SEI-layer consisting of Li₂EDC with varying thickness. The thickness and voltages are indicated in each panel, with the density of Li (red lines) and carbonate groups (by center of mass, green lines) in the SEI and the EC (by center of mass, blue lines). Reproduced from ref (340). Copyright 2013 American Chemical Society.

Figure 38

Representative Li⁺-anion coordination in the double layer at (a) negative electrode potentials and (b) positive electrode potentials. Red and blue bars indicate the (x,y) plane of the anode and cathode surfaces, respectively. Reproduced from ref (243). Copyright 2016 American Chemical Society.

Figure 39

Model system and reaction scheme (cyan, carbon; red, oxygen; white, hydrogen; orange, phosphorus; green, fluorine; blue, lithium). Reproduced from ref (341). Copyright 2014 American Chemical Society.

Figure 40

Schematic representation of the 4 types of events selected to occur in the kMC simulations. Reproduced with permission from ref (343). Copyright 2011 The Electrochemical Society.

Figure 41

Temperature effect on the total charging time and SEI thickness during a full charging process. Reproduced from ref (344). Copyright 2017 American Chemical Society.

Figure 42

Effect of SEI solvent activation energy on the total charging time and SEI thickness during a full charging process. Reproduced from ref (344). Copyright 2017 American Chemical Society.

Figure 43

Calculated electrode potential during the first charge of the formation process for electrodes with particle radius R1 = 3 × 10^–6 m and R2 = 10 × 10^–6 m. Dashed gray lines correspond to the potentials 0.55, 0.525, and 0.5 V from top to bottom. Reproduced with permission from ref (345). Copyright 2017 The Electrochemical Society

Figure 44

Schematic of the single-particle model used to simulate SEI film growth by kMC. Reproduced with permission from ref (346). Copyright 2017 The Electrochemical Society.

Figure 45

SEI film thickness under OCV as a function of time for various initial SOC. The simulations assumed an initial SEI thickness of 1 nm. Reproduced with permision from ref (347). Copyright 2011 Elsevier Ltd.

Figure 46

Experimental and simulated Nyquist plots for Si-anode-based half cells. Reproduced with permission from ref (350). Copyright 2017 Elsevier Ltd.

Figure 47

Different structural models for LIB electrode proposed in literature. (a) Assumed main geometrical features. Reproduced from ref (359). Copyright 2010 The Electrochemical Society. (b) Reconstructed by tomography (in blue). Reproduced with permission from ref (361). Copyright 2014 Elsevier Ltd.

Figure 48

Representation of the volume averaging method applied to a porous medium made of a solid phase, σ, filled with a liquidphase, β. The scale of the entire domain is L, while the pore-scale characteristic length is l_β. Reproduced with permission from ref (363). Copyright 2013 Elsevier Ltd.

Figure 49

Schematic representation of the agglomerate scale (a scale between the crystal and electrode size). Reproduced with permission from ref (366). Copyright 2015 The Electrochemical Society.

Figure 50

Scheme of the MSM of a Li-S cell proposed by Thangavel et al. Reproduced with permission from ref (372). Copyright 2016 The Electrochemical Society.

Figure 51

(a) Schematics of the LOB cell model. (b) One of the surface area distributions (related to a pore size distribution) and electronic tunneling function adopted in the model. (c) Simulation results for Super P and Ketjen Black carbon-based positive electrodes for two discharge current densities. Reproduced with permission from ref (381). Copyright 2014 The Electrochemical Society.

Figure 52

Macroscopic and microscopic scale models with multiscale coupling variables shown. Reproduced with permission from ref (385). Copyright 2014 The Electrochemical Society.

Figure 53

Simulated lithium concentration within the cathode at various states of charge (electrolyte-black, carbon-binder domain-dark gray). (a) SOC 0 (initial state). (b) SOC 0.1. (c) SOC 1 (fully charged). Reproduced with permission from ref (361). Copyright 2014 Elsevier Ltd.

Figure 54

Scheme of different physicochemical processes taking place in a Li–S positive electrode which may affect the cell electrochemical impedance spectra: (a) reactions at the pore-scale, (b) precipitate film and electrochemical double layer, (c) imperfect contact between the precipitate film and the carbon substrate, and (d) electrode microstructure. (e) Interfacial impedance contributions are represented in terms of an equivalent electrical circuit. (f) Schematic illustration of a simulated electrochemical impedance spectrum, (g) and example of its validation against the experimental one. Reproduced from ref (394). Copyright 2017 American Chemical Society.

Figure 55

Three-dimensional calculation of Li₂O₂ formation in a reconstructed carbon-based Li-O₂ battery. Reproduced with permission from ref (395). Copyright 2015 American Chemical Society.

Figure 56

Workflow of the pore network approach. Reproduced from ref (396). Copyright 2018 American Chemical Society.

Figure 57

(a) Electrode density effect on its capacity. Reproduced with permission from ref (398). Copyright 2003 Elsevier Ltd. (b) Schematic representation of the typical LIB electrode fabrication steps at the industrial level.

Figure 58

DEM model of a LIB cathode calendaring for different hypothesis of particles morphologies. Reproduced with permission from ref (405). Copyright 2015 Elsevier Ltd.

Figure 59

Scheme of the computational workflow proposed in ref (408) to simulate the LIB electrode fabrication and its impact on the cell performance. Coarse grained molecular dynamics is used to simulate slurries, the coating, and the effect of solvent evaporation, whereas the resulting predicted electrode mesostructure is used in a 3D-resolved continuum model to simulate the cell operation. Adapted from ref (408). Copyright 2017 American Chemical Society.

Figure 60

Predicted through thickness strain distribution at the end of discharge (2000 s) under a fixed-free (10 psi) boundary condition using a quarter-sized 2D RVE model. The areas with a voidlike appearance have a strain value beyond the range and hence are removed from the display. Particles 2 and 4 have a close-packed pattern, whereas particles 1, 3, and 5 have a loose-packed pattern. The maximum strain in the separator is found near a loose-packed particle of larger diameter. Reproduced with permission from ref (411). Copyright 2010 Elsevier Ltd.

Figure 61

Length scale dependent physics impacting battery modeling. Reproduced with permission from ref (9). Copyright 2016 Elsevier Ltd.

Figure 62

Comparison between temperature contour of the model result and the test. Reproduced with permission from ref (418). Copyright 2014 Elsevier Ltd.

Figure 63

(a) MSM architecture and (b) temperature distribution in an unrolled (left) and a rolled cylindrical (right) cell subject to 2C discharge current. Reproduced with permission from ref (422). Copyright 2014 Elsevier Ltd.

Figure 64

(a) Distribution of lithium concentration inside all particles. Horizontal axes represent particle location in the electrode by the x position along the electrode thickness direction and by the y position in the plane of the electrode. The vertical axis represents a point in the particle by the radial r position. (b) Lithium concentration in particles at y = 0. (c) Radial stress distribution in particles. (d) Tangential stress distribution in particles. Results shown are at the instant of 300 s. (e) Temporal profiles of the maximum tensile radial stress and the maximum compressive tangential stress for particles next to the separator (SP) and next to the current collector (CC). Reproduced with permission from ref (429). Copyright 2017 Elsevier Ltd.

Figure 65

Schematics of an imaginary computational approach combining MSM with machine learning. “X”, “Y”, and “Z” refer to some characteristics being used as an input in the study (e.g., atomic composition of an active material).