PHYSICS-BASED SIMULATIONS OF ELECTROCHEMICAL IMPEDANCE SPECTRA FOR LITHIUM-ION BATTERY ELECTRODES By Danqi Qu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Materials Science and Engineering – Doctor of Philosophy Computational Mathematics, Science and Engineering – Dual Major 2022 ABSTRACT Electrochemical impedance spectroscopy (EIS) is a powerful and non-destructive character- ization technique widely used in the electrochemical research field. It can measure many macroscopic properties such as internal resistance, capacitance, and diffusivity by fitting the obtained impedance with equivalent circuits. Each of the acquired quantities reflects an electrochemical mechanism, e.g, charge-transfer reaction, double layer formation, and mass transport, taking place in the electrode. However, the obtained quantity is a total value for the whole electrode. The underlying connections between the macroscopic properties, intrinsic material parameters, and electrode microstructures are not well understood. This dissertation focuses on building a modeling framework to simulate EIS processes with given electrode microstructures and intrinsic material parameters. With this simulation tool, we provide a digital bridge between battery electrode material properties, electrode microstruc- tures, and their corresponding EIS impedance. Capacitance of an electrochemical device originates from double layer formation in the electrolyte. However, there is a huge spatial discrepancy between the dimensions of double layer and electrode particles (or interparticle space). Thus, smoothed boundary method and adaptive mesh refinement are used to handle the scale discrepancy and the complex geometries of electrode particles in solving the Nernst-Planck-Poisson equations in simulating the double layer formation under voltage loading. The obtained double-layer capacitance is incorporated into multiphysics electrochemical simulations. Cathode electrode made of Nickel-Manganese-Cobalt (NMC111) oxide, is examined with this simulation tool. As a solid solution material, lithium transport in the NMC111 electrode particles is described by Fick’s law. EIS curves for various conditions, including different states of charge, electrolyte salt concentration, electrode microstructures, are extracted from the simulations and analyzed. The simulations properly reflect the relationships between particle exchange current density, reactive surface area, and the total resistance of the electrode. Anodes made of graphite, a phase-transforming material upon lithiation/delithiation, are also examined using the simulation tool. The Cahn-Hilliard equation is employed to model the phase transformation processes in the particles. EIS simulations are conducted on single-phase and multi-phase graphite. For single-phase or core-shell phase-distributed graphite particles, the simulated EIS curves exhibit a typical semicircle with a Warburg part. Interestingly, if phase boundaries intersect particle surfaces, a low frequency inductive loop appears on the EIS curve. Lastly, the simulation tool is applied to simulate EIS processes of a full-cell battery of both cathode and anode microstructures. On each electrode, the total current is comprised of capacitance and reaction currents. It is observed that, depending on the loading frequency, the ratio of capacitance-to-reaction current on the two electrodes can be significantly differ- ent. The simulation tool allows us to examine the details of electrochemical processes during EIS measurements. Copyright by DANQI QU 2022 This thesis is dedicated to my parents. v ACKNOWLEDGEMENTS My time at Michigan State University is the most challenging and wonderful period of my life so far. Without the help and support of many individuals, it is impossible for me to reach where I am today. First, I would like to thank my academic advisor, Dr. Hui-Chia Yu. I feel very fortunate to be his first student, which allowed me to help him establish the research platform at MSU and explore new research fields. On academic aspect, he always keeps a high standard and constantly challengs my critical thinking. He leads me to be the best scientific researcher I can be. He never hesitates to spend time with me on discussion and advising. He is always patient and generous to give support when I touched my blind spots in knowledge. For every paper and proposal, he will read and revise carefully sentence by sentence. Other than research, he truly cares about his students’ personal development. He often encourages me to travel to conferences and to build my network. I appreciate Dr. Yu’s teaching and caring over this five years. The experience in his group will be one of the most memorable time of my life. Secondly, I would like to thank my doctoral committee. I would like to thank Dr. Wei Lai for his valuable suggestions in my research. He has some experience on electrochemistry before and is willing to provide all he knows when I am confused. I would like to thank Dr. Yue Qi who was my committee member before she moved to Brown University. She was also the one who sent me my formal Ph.D. offer from MSU. I still remember her words, “Don’t be shy to ask questions and help, we are your resources.” I would like to thank Dr. Swain Gerg for his valuable insights on my dissertation. He raised lots of questions regarding my work from his experimental view during my comprehensive exam. I would like to thank Dr. Michael Murillo, who is one of my committee members from CMSE department. I learned a lot from his classes on machine learning for physics, which helped me find my first internship. I would like to thank thank Dr. David Hickey, who was willing to join my committee after my comprehensive exam. With his expertise in electrochemistry, my dissertation is strengthened from different perspectives. vi I would also like to thank my colleagues, Robert Termuhlen, Affan Malik, for providing valuable insights from the discussions and the questions in our group meetings. I would like to thank the ChEMS and CMSE office staff, Tiffany Owen, Donna Fernandez, Brad Tobin, Jennifer Keddle, Heather Dainton, Jessica Gallegos, Heather Williamson, Melinda McEwan, and Lisa Roy for helping me with various issues from administration to conference travel. I would like to thank my friends Runze Su, Alex Mckim, Jiyun Park, Jiawei Lu, Yining He, Shengyuan Bai, Yun Cheng and Fang Wu for our good time beyond research and classes. Finally, I would like to express my love and gratitude to my family. My father influenced and motivated me to pursue science research ever since I was a child. He is proud of me for what I have accomplished. My parents have been generously supporting me in all of my important periods with their selfless love and care. vii TABLE OF CONTENTS CHAPTER 1 BACKGROUND . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 NUMERICAL METHOD AND TECHNIQUES FOR SOLVING DIF- FERENTIAL EQUATIONS IN SIMULATIONS . . . . . . . . . . . . 23 CHAPTER 3 DOUBLE LAYER FORMATION ON ELECTRODES WITH COM- PLEX MICROSTRUCTURE . . . . . . . . . . . . . . . . . . . . . . 46 CHAPTER 4 PHYSICAL-BASED SIMULATION OF ELECTROCHEMICAL IMPEDANCE SPECTRUM ON CATHODE MATERIAL . . . . . . . 68 CHAPTER 5 PHYSICAL-BASED SIMULATION OF ELECTROCHEMICAL IMPEDANCE SPECTRUM ON GRAPHITE . . . . . . . . . . . . . 100 CHAPTER 6 FULL CELL SIMULATION . . . . . . . . . . . . . . . . . . . . . . . 137 CHAPTER 7 SUMMARY AND OUTLOOK . . . . . . . . . . . . . . . . . . . . . . 160 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 viii CHAPTER 1 BACKGROUND 1 1.1 Motivation Nowadays, with pressing world environmental concerns, clean energy has received extended attention. Among them, a clear trend is that electric vehicles will gradually replace internal- combustion-engine vehicles. The market demands of advanced electric vehicles raise in- creasing requirements for electrochemical energy storage systems, especially for lithium-ion batteries. To improve the range, charging time, lifespan, and safety of electric vehicles, re- search efforts have been put on enhancing the capacity, rate performance, cyclability, and reliability of lithium-ion batteries. A lithium-ion battery consists of two electrodes: anode and cathode that are separated by an insulating layer to prevent internal shortcuts. Electrolyte fills the interparticle space in the electrodes and pores in the separator to allow lithium (Li) ions shuttling between the two electrodes. A battery’s macroscopic performance is determined by many microscopic pro- cesses simultaneously occurring in the electrodes. For example, deintercalation of cathode (also called charging process) involves Li outward diffusion in the particles, electron migra- tion through the particle network, electrochemical reaction on particle surfaces, and Li-ion diffusion in the electrolyte, all of which take place in a complex cathode microstructure. The Li ions leaving the cathode migrate through the electrolyte, reach the anode particle surfaces, and are then inserted into anode particles. As can be envisioned, microstructural features, such as porosity [1], tortuosity [2] of interparticle space, reactive surface area [3], and particle size [1, 4], all affect a battery’s performance. These microstructural features sometimes vary through the electrode, which can result in substantially non-uniform reaction current density and might deteriorate cell performance or even cause safety concerns [5]. Without a com- prehensive understanding of the relationship between electrode microstructure and battery performance, it is impossible to improve electrode designs or optimize battery performance. Electrochemical impedance spectroscopy (EIS) is a widely used technique to measure properties of electrochemical devices, such as batteries and fuel cells. The EIS measures the response current (or voltage) to an oscillating voltage (or current) loading [6]. The 2 obtained EIS curve is fitted with an equivalent circuit model (ECM) that comprises of re- sistors, capacitors, and constant-phase elements [7], such that the total values of resistance, capacitance, and Warburg impedance are evaluated from the fitting. These resistance, ca- pacitance, and Warburg impedance correspond to the physical mechanisms of charge transfer reaction, polarization, and diffusion (mass transport) involved in the whole electrochemical process. However, the connection between the acquired values and the microscopic informa- tion is absent in the ECM fitting. It cannot provide physical insights of microstructural-level electrochemical processes. To investigate the relationship between macroscopic electrode performance and microscopic processes, complex microstructures must be explicitly consid- ered, which can be realized using electrochemical simulations [8, 9, 10]. Thus, this work aims to create a simulation tool based on physical mechanisms taking place in electrode microstructures to unravel the electrochemical processes during EIS measurements. The key components of this simulation framework are electrochemical double layer formation and coupled multiphysics electrochemical dynamics in EIS processes. The relevant backgrounds are reviewed in the following sections. 3 1.2 Electrochemical Double Layer Double layer is the key part in electrochemical capacitor and has a large effect on the EIS behavior of liquid-electrolyte batteries. This section will provide a literature survey of the formation of double layer and the resulting electrochemical capacitor. The following para- graphs contains brief descriptions of terminologies used in this thesis and a summary of published research progress related to the content of this thesis. 1.2.1 Electrochemical Capacitor Electrochemical capacitor (EC), also known as electrochemical double layer capacitor, su- percapacitor and ultracapacitor, is one type of electrical energy storage devices. A simple EC can be represented by two electrodes connected by an external circuit and immersed in a electrolyte, for example two metal plates in salt water. When applying external voltage to these two electrodes, an electron current flows from one to the other via the external circuit. To counter the electropotential polarization, ionic charge separation takes place at each side of the liquid-solid interface since electrons can not enter the electrolyte and ions can not enter the electrodes. Once the external circuit is cut off, the voltage persists and the energy is stored. In this state, the ions in the electrolyte are attached to the surface of each electrode by an equal but opposite charge. This status is maintained by electrostatic attractive force. The thin region with charge separation near the surfaces of electrodes is named “diffuse double layer”, which is the working mechanism of EC’s energy storage. Thus, understanding the double layer formation process is of great importance in developing physics-based EIS simulations. Conway described the difference between supercapacitor and battery behavior in mecha- nism and performance of electrochemical energy storage [11]. While the structure of ECs is somehow similar to batteries, their working mechanisms for storing energy are significantly different. In batteries, electrochemical energy is stored via chemical reactants by redox re- 4 actions taking place on the interfaces between electrode and electrolyte. A charge transfer process, or called Faraday process, occurs on the surface of electrode particles during charging or discharging. In many cases phase transformations occur in the bulk of electrode particles. Those processes usually lead to side reactions, and thus battery performance decays and their cycle live are limited. Additionally, the amount of energy that a battery can store depends on the quantity of electroactive materials that can accommodate the charge-transfer reactions. In contrast, in an electrochemical capacitor, electrostatic energy storage is achieved by the amount of charge separation, Q, and the stored energy is ∆G = CV 2 /2 = QV /2, where C is the capacitance and V is the electrostatic potential difference between the two electrodes. The cell performance curve of an EC shows a linear (or almost linear) function of the de- gree of charge, which is essentially different from batteries that shows a hyperbolic-sine-type curve determined by the thermodynamic behavior of single-phase battery reactants. Such a difference reflects the different work mechanisms of these two devices. Note that the working voltage of an EC or a battery must be lower than the breakdown potential of the electrolyte (namely, where electrolysis occurs): typically, < 1 V for aqueous or < 3 V for organic electrolyte [12]. ECs have some common features with electronic capacitors. Both of them store electric energy by charge separation, but the media between two electrodes normally contains no movable charged particles for electronic capacitors. When plotting the power density of dif- ferent basic energy storage devices as a function of energy density (see Figure 1.1), EC is at the middle between electronic capacitors and batteries and can fairly operate in a wide range of power density and energy density. Electronic capacitors have very high power and very fast response time and almost unlimited cycle life while their energy density is as low as 0.1 Wh/kg or even lower [13]. Hence electronic capacitors are often used in the situations that need fast response and to detect tiny signals but not require significant energy storage. Compared to an electronic capacitor, an electrochemical capacitor has a higher energy den- sity. ECs can have specific capacitances exceeding 200 F/g, cycle lives above 1 × 106 cycles. 5 Figure 1.1 Power density as a function of energy density for various energy storage devices(figure taken from [13]) If including pseudocapacitances, the total capacitance will be much higher as pseudocapaci- tances can contribute more than 100 times double-layer capacitance [14]. Pseudocapacitance is a process of battery-like energy storage mechanism. This means it involves charge transfer processes, but behaves like a electrochemical capacitance. Since the bolder between pseudo- capacitances and battery behavior is still subtle, Simon and Gogotsi [15] gave their principles to distinguish pseudocapacitances and battery behaviors. Using fast and highly reversible redox reactions at the surface of active material, is defined as the pseudo-capacitive behavior. Transition metal oxides such as RuO2 , MnO2 as well as electronically conducting polymers are usually used as electrodes for pseudo-capacitors. However, this thesis does not consider pseudocapacitances in the modeling and simulations. ECs have already been found to have a wide range of applications in hybrid energy storage 6 Figure 1.2 Schematic diagram of the electrochemical desalination device (figure taken from [16]) system [17], desalination [18] and biochemical system design [19]. An illustration of using ECs in water desalination is shown in Fig. 1.2. The capacitance of a modern EC can range from µF to kF. ECs are used in the places where rapidly varying load power profiles are needed such as hybrid electric vehicles. Influenced by battery technologies, porous carbon electrodes and propylene carbonate electrolyte that are commonly used in Li-ion batteries are directly used in ECs. Nevertheless, little attention has been paid on modeling ECs. For example, there have been numerous modeling studies applying to batteries [20, 21, 22, 23, 24] but only few was specific for ECs. Since the diffuse double layer, the ECs’ energy storage mechanism, also plays an important role in battery’s EIS behavior, developing a new continuum modeling and simulation capability of double layers is essential to this thesis work. This requires an understanding of the fundamental principle of the physical and chemical processes that occur 7 in ECs. 1.2.2 Physical Models of Double Layer Structure Helmholtz [25] was the first to realize that a charged electrode would form an ionic interface on its surface and called such a surface structure as the electric double layer in 1853. In Helmholtz’s model, all counter-ions are assumed to directly attach to the electrode surface. This structure is analogous to that of conventional dielectric capacitors with two planar electrodes separated by a distance H. See Fig. 1.3(a) for an illustration. Figure 1.3 Schematics of the electric double layer organization revealing the distribution of solvated anions and cations around the electrode/electrolyte interface. (a) Helmholtz model, (b) Gouy-Chapman model, and (c) Gouy-Chapman-Stern model. (figure taken from [26]) Gouy and Chapman (1910) developed an electric double layer model accounting for the fact that the ions are mobile in the electrolyte solutions and are dirven by the coupled influence of diffusion and electrostatic forces. This results in the so-called diffuse layer shown in Fig. 1.3(b). In this model, the ions are treated as point-charges and the equilibrium 8 concentration ci of ion species ‘i’ is given by the Boltzmann distribution as zi eϕ ci = ci∞ exp(− ) (1.1) kB T where zi and ci∞ are the valency and bulk molar concentration of ion species ‘i’, respectively. The absolute temperature is denoted by T , e is the elementary charge, while kB is the Boltzmann constant (kB = 1.381 × 10−23 m2 kgK−1 s−2 ). In the Gouy-Chapman model, the local electric potential ϕ in the diffuse layer is determined by the Poisson’s equation assuming a constant electrolyte permittivity. Stern (1924) combined the Helmholtz model and the Gouy-Chapman model and described the electric double layer as two layers seen as Fig. 1.3(c). The Stern layer referring to the compact layer of immobile ions strongly adsorbed to the electrode surface and the diffuse layer is where ions are mobile and the Gouy-Chapman model applies. Note that there are no free charges within the Stern layer. The total electric double layer capacitance consists of the Stern layer and diffuse layer capacitance in series. Today, the Stern model is considered the standard model of double layer structure. For most of electrolytes in practical applications, the predicted thickness of the diffuse double layer is less than 50 nm. The models above are only valid for dilute solutions and low electric potential because they didn’t count finite ion size. In actuality, ions have finite ion size and cannot be treated as point charges. Bikerman (1942) developed the first equilibrium modified Poisson-Boltzmann model accounting for finite ion size [27]. The ion concentration is limited to the concentra- tion of close packed structure which is given by cmax = 1/(NA a3 ) if we assume it’s simple cubic packing. Here NA is the Avogadro constant and a is the diameter of ions. There- fore, the equilibrium ion concentration given by Equation (1.1) cannot exceed cmax . The corresponding maximum surface potential can be calculated by kB T ϕmax = ln(NA a3 c∞ ). (1.2) zi e The magnitude of local electric potential |ϕ| in the diffuse layer should not exceed ϕmax referenced to the potential at infinity (in the deep bulk of the electrolyte). For instance, a 9 = 0.3 nm, c∞ = 10−5 M and z = 1 in solvated bulk K+ –Cl− interactions, then ϕmax = 0.33 V [28]. Figure 1.4 Two flat plates at constant surface potentials of ψa and ψb (figure taken from [29]) Analytical solution exists for the Stern model on one-dimensional case and for some simple two-dimensional geometries. For example, the Poisson equation for the case of parallel plates is given by d2 ϕ ϵ 2 = −ρ (1.3) dx where ϵ is the dielectric constant and ρ is the charge density. After applying Debye-Huckel approximation [29], it becomes d2 ϕ = κ2 ϕ (1.4) dx2 p where κ = 2e2 z 2 c∞ /(ϵkB T ). The solution of Eq. (1.4) gives ϕ = A1 cosh κx + B1 sinh κx (1.5) 10 Making use of the boundary conditions of ϕ at x = ±h/2 as shown in Fig. 1.4, where (ϕa − ϕ(h/2))/λ = ϕ′ (h/2) and (ϕb − ϕ(−h/2))/λ = −ϕ′ (−h/2), the coefficients are determined as ϕa + ϕb ϕa − ϕb A1 = and B1 = (1.6) 2 cosh κh/2 + 2λκ sinh κh/2 2λκ cosh κh/2 + 2 sinh κh/2 where λ represents the Stern layer thickness which is often several layers of ions thick. 1.2.3 Simulations of Double Layer Formation Due to the importance of double in many applications, simulations of double layer formation have been developed and performed using several methods, such as equivalent-circuit-based, atomic simulation-based, statistics-based, macrohomogeneous-electrochemistry-based tech- niques, etc. The equivalent RC circuit models and more complex transmission line models are used to investigate the performance of supercapacitors and batteries, in which the dou- ble layer is represented by a capacitor component. However, the equivalent circuit models (ECM) require prior knowledge of macroscopic parameters such as the resistance and ca- pacitance of the devices, which are usually determined by experiments or other empirical methods. In fact, the ECM is typically used to fit the experimental data to extract an ef- fective capacitance of the device [30, 31, 32] rather than to predict or explain the formation of the double layer. In statistical mechanical theories, the electric double layer is often simulated by the primitive model where the ions are treated as charged hard spheres. The solvent and the electrode are represented as dielectric continuum with the same dielectric constant because Poisson-Boltzman theory assumes that the dielectric constant has no influence on the ionic profile. Henderson et al. [33] used a classical density functional theory (DFT) to study the electric double layer formed by charged hard spheres near a planar charged surface. Their DFT predictions are found to be in good agreement with the results from Monte Carlo simu- lations. The primitive model can also applied to model electrodes with a dielectric constant different from that of electrolyte, e.g., a value approach to infinity for metal electrodes. It 11 was shown by Torrie et al. [34] that the assumption of Poisson-Boltzmann theory is not valid, especially when the electrolyte contains multivalent ions. Nagy et al. [35] modified the primitive model by adding an inner layer (or called Stern layer in this thesis). Their Monte Carlo simulation results reported that the inner layer with lower dielectric constant than that of the bulk electrolyte, bringing the results into agreement with experiments. They also showed the ionic profile of diffuse layer depends on the dielectric constants of both inner layer and diffuse layer. Bhuiyan et al. [36] reported a Monte Carlo simulation result based on a new density functional theory, in which each cation is represented by two touching hard spheres: one positively charged and the other neutral, and each anion is still represented by a negatively charged hard sphere. However, due to the physical scale of those statistical me- chanical simulations, they were all equivalent to one-dimensional cases and cannot account for the 3D irregular morphology/geometry of electrodes. Homogeneous models were developed to investigate the dynamics of double layer for- mation and charge/discharge processes of supercapacitors. Those models treat the hetero- geneous microstructure of electrodes as a homogeneous medium with effective macroscopic properties such as porosity, tortuosity and capacitance per area or per volume. Newman pioneered the theory of porous electrodes to treat the complex porous solid as a homoge- neous phase [37]. He and coworkers successfully fitted theoretical current-time curves with experimental data and determined the capacitance of double layer. Bazant et al. further implemented the porous electrode theory by nonequilibrium thermodynamics methods [38]. Additionally, they treated porosity accounting for interparticle void and intraparticle micro- pores differently by a modified-Donnan approach when dealing with electrodes composed of primary particles that are porous themselves [39]. Despite the simplicity of using macroho- mogeneous properties to simulate macroscale electrochemical processes, they cannot consider the detailed three-dimensional electrode morphology and microstructure. Simulations based on Poisson-Boltzmann theory assume that the profile of ions in the diffuse layer obeys Boltzmann distribution when the system reaches steady-state. Such a 12 treatment makes predicting double-layer structure at equilibrium convenient. For example, the diffuse layer around a nm-sized micropore in porous electrodes has been studied by solving Poisson-Boltzmann equation [40] although micropores in that work were still viewed as cylin- ders. Other simulations accounted for the finite ion size [28, 41] and electrode morphology have been reported. These studies largely extended the applicability of Poisson-Boltzmann theory into cases with high applied voltage and concentrated solutions. However, one of the limitations is that they cannot reveal the dynamic process of the double layer formation. And many of these simulations are still confined to one-dimension or two-dimension with regular geometries. The classical Nernst-Planck-Poisson (NPP) model is physics-based on the ionic trans- port driven by diffusion and electromigration, in which the electrostatic potential field is determined by the charge density distribution. Thus, it can describe the ionic concentration evolution in the diffuse double layer and the transient electric potential evolution. Details of the model will be introduced in the next section. While the Nernst-Planck-Poisson model is suitable to track the dynamics of double layer formation, most of the related simulation work was limited to one-dimensional problems; for example, the 1D dynamic simulation of double layer formation between two planar metal plates reported by Bazant et al. [42]. The main difficulty of implementing multidimensional NPP simulations is to resolve the extremely thin diffuse double layer (typically less than 50 nm) while keeping track of the ionic concentration and electrostatic potential evolution spanning several order of magnitudes larger than the diffuse layer (e.g., few or tens of µm space between particles and between electrodes). Very fine mesh for the diffuse layer region will be required for numerical methods to solve the governing equations. In Chapter 2, we will introduce an approach that allows a straightfor- ward reformulation of governing equations that can be solved on non-body conforming mesh on complex electrode microstructures. The reformulated equations are solved on adaptively refined grid systems that have fine grids near the diffuse layer regions and coarse grids in the particle or electrolyte bulk phases. Thus, it allows us to simulate the dynamics of double 13 layer formation on complex 2D and 3D electrode particles which is difficult for the conven- tional sharp-interface numerical methods that require generating mesh systems conformal with the complex microstructure. 14 1.3 Electrochemical Impedance Spectroscopy and Electrochemi- cal Modeling for Batteries There are electrochemical characterization techniques widely employed to study battery electrodes, such as potentiostatic intermittent titration (PITT) [43], galvanostatic intermit- tent titration (GITT) [44], cyclic voltammetry (CV) [45], and electrochemical impedance spectroscopy (EIS) [6]. In this work, we place our focus on electrochemical impedance spectroscopy (EIS) because it involves an electrode’s charge-transfer reaction, double-layer formation, and diffusion/migration process simultaneously in one measurement. A detailed modeling will reveal those coupled multiphysics mechanisms occurring in complex electrode microstructures. This knowledge is crucial for optimally designing battery electrodes. As mentioned earlier, EIS is widely used to measure electrodes’ (or batteries) resistance, capac- itance, and Warburg impedance, which correspond to the charge transfer reaction, polariza- tion, and mass transport, respectively, during the measurement process. Since each of those physical responses has a different time scale, as the source frequency sweeps from low to high values (or vice versa), the dominant mechanism will be separated at its corresponding frequency. Thus, values of different properties can be obtained from the EIS curve. Electrochemical processes can be controlled by voltage (potentiostat) or current (gal- vanostat) loadings, which maintains a potential difference between the working electrode and a reference electrode or maintains a current through the two electrodes, respectively. A current will be measured as the response for the potentiostat mode and a voltage will be the response for the galvanostat mode. EIS measurements can be conducted by either mode. In potentiostat experiments, a voltage adder is used to combine the direct-current potential corresponding to the polarization point (equilibrium potential) and the alternating-current potential generated by the frequency function analyzer [6]. Currents are measured as the output signal. Both the input and output signals are oscillatory in time domain, and they are related by a transfer function of frequencies. If the transfer function takes the form of a ratio of potential over current, the transfer function is called an impedance. The amplitude 15 of the source signal is usually small to ensure that the system is not perturbed to be far away from the equilibrium, thus returning a linear response [46]. However, nonlinear EIS is sometimes also utilized to investigate more complicated situations [47]. Figure 1.5 Bode plot and Nyquist plot for expressing EIS of an RC circuit [48]. The value of impedance depends on the loading frequency and is a complex quantity that can be expressed by an amplitude (absolute value) and a phase shift between the input and output signals. Bode plot and Nyquist plot are the two most common ways to express EIS results. Bode plot contains two pieces of information in one figure. The abscissa is frequency in logarithmic scale and one ordinate is the amplitude of the impedance Z in logarithm scale while the second ordinate is the phase shift ϕ. The lower left panel in Fig. 1.5 illustrates a Bode plot of an RC circuit (a capacitor in parallel with a resistor). The advantage of Bode plot is that all information is displayed together. At high frequencies, the phase shift is approaching to 90 degree while the impedance amplitude is down to nearly zero. On a ′′ Nyquist plot, negative of the impedance’s imaginary part (−Z ) is plotted versus the real ′ part of the impedance (Z ). In this case, the distance between an impedance point and the origin is the impedance’s amplitude, and the angle between the that connecting line and the abscissa is the phase shift. The lower right panel in Fig. 1.5 illustrates a Nyquist plot. At high frequencies, the impedance is close to the origin, also showing a 90-degree phase 16 delay and zero impedance as in the Bode plot. The low frequency impedance is away from the origin. At zero frequency, the phase delay is zero. While it is difficult to directly read impedance values from a Nyquist plot, differences between different EIS curves can be easily noticed on a Nyquist plot. The shapes of EIS curves on a Nyquist plot can also be quickly recognized to be related to relevant equivalent circuits. Thus, it is widely used in presenting EIS results. Additionally, some parameters can still be read from the plot. For instance, the resistance of RC circuit can be obtained by the diameter of semicircle as on lower right panel of Fig. 1.5. One of the applications of EIS measurements is to estimate battery’s degradation [49]. The internal resistance of electrode can increase either due to the growth of solid electrolyte interphase [50, 51], crack of electrode particles [52], or change of crystal structure of the host crystals [53, 54]. Such increases in resistance over cycling can be easily observed on EIS curves. However, only a total value of the whole system can be obtained, which may be convoluted with changes in both the intrinsic material properties and the microstructures. However, any further information to separate the contributions from material properties and microstructural features is absent in the EIS results. As mentioned, although ECM can be used to extract batteries’ overall properties, it does not provide physical insights of electrochemical processes. Such insights can only be obtained by physics-based electrochemical simulations [8, 9, 10] that explicitly consider the relevant electrochemical processes. However, due to the complexity of electrode microstruc- tures, porous electrode theory (PET) [55, 56], also known as pseudo-2D (P2D) modeling, is currently the most prevalent simulation tool for studying electrochemical energy storage or conversion devices. The PET models treat porous electrodes as homogeneous media that have uniform microstructural features, such as tortuosity of interparticle space, porosity, and reactive surface area density, throughout the entire electrodes. As a result, Li salt diffusion and electron migration through an electrode are calculated using only uniform effective salt diffusivity and effective electric conductivity of the electrode. The Li intercala- 17 tion/deintercalation is modeled by diffusion in/outward of spherical particles with uniform radii but uniformly located at different places in the electrode. Thus, 3D electrode simula- tions are decoupled to three 1D simulations (salt concentration, electrostatic potentials in the solid and electrolyte) along the electrode thickness direction and many 1D simulations (Li concentration) along the particle radial direction [57]. This types of models are also called P2D models because there are only two coordinate directions: one along the electrode thick- ness, x, and the other along the radial direction, r. See Fig. 1.6 for an illustration. Although electrochemical simulations are greatly simplified in the PET models [58, 59, 60], all mi- crostructural phenomena are absent. For example, small and large particles will (de)lithiate differently and salt diffuses non-uniformly in regions with different porosities and tortuosities in real electrode microstructures. Those microstructural phenomena are difficult to be con- sidered in PET models. To interpret the significance of EIS measurements, electrochemical processes in the electrode microstructures must be resolved. Figure 1.6 Schematic of the PET model for an LiC6 /LiCoO2 cell Nowadays image processing techniques have well advanced such that three-dimensional (3D) microstructure reconstructions become mature experimental methods to characterize electrode materials. Both focus-ion-beam scanning-electron-microscopy (FIB-SEM) [8, 61, 62] and nano/micro X-ray computing tomography (µ-CT) [9, 63, 64] are widely utilized in 18 the battery research community for 3D electrode microstructure reconstruction. As opposed to the destructive FIB-SEM technique, nano/micro X-ray CT is a non-destructive method, which has been employed to in situ study electrode’s microstructural evolution, porosity, and tortuosity variation [64]. The abundant microstructure data open a window for directly simulating the electrochemical processes in complex electrode microstructures. To investigate the relationship between macroscopic electrode performance and micro- scopic processes, complex microstructures must be explicitly considered in electrochemical simulations. In the conventional sharp-interface modeling, governing equations are numeri- cally solved on mesh systems that conform with the domain where the equations are solved; for instance, as in typical finite element method (FEM) [61, 65] or finite volume method (FVM) [8, 63] simulations. However, since electrodes are made of packed particles with very tortuous interparticle space, generating mesh conformal with those complex microstructures are very time-consuming. We acknowledge that mesh generation for electrode microstruc- tures are possible, as being demonstrated using commercial software packages such as COM- SOL [65], Ansys [8], GeoDict [62], etc. Multiphysics phenomena coupled with electrochemical processes and EIS response can be simulated using conventional methods [66]. However, in many cases, manually fixing thousands of broken elements is required, which is extremely tedious [8, 67]. As such, mesh generation is still usually the most time-consuming task in microstructure simulations. Recently, Shodiev et al. [68] demonstrated a framework of using coarse-grained molecular dynamics method to generate a composite porous electrode, mim- icking the slurry-drying process of fabricating NMC (Lix Ni1/3 Mn1/3 Co1/3 O2 ) cathodes, and split voxels that represent the microstructure into elements for FEM simulations [69]. They performed a set of electrochemical simulations with explicit considerations of three phases (NMC active material, carbon-binder domain, and electrolyte in pores) in an composite elec- trode and extracted EIS curves from the simulations. However, this method may also face the challenge as in other voxelization-mesh-generation processes, such as a large number of voxels. 19 As will be presented in Chapter 2, the diffuse interface approach developed for simulating double layer formation is also employed to circumvent the requirement of using conformal mesh in simulating the electrochemical processes in complex electrode microstructures. The capacitance due to double-layer formation is calculated based on the Nernst-Planck-Poisson model. The obtained capacitance is incorporated into the electrochemical simulations to ex- tract the total current responding to oscillating cell voltage loadings. Thus, this simulation tool allows a very efficient implementation of physics-based microstructure-level electrochem- ical impedance simulations. The simulated EIS curves on the Nyquist plot well agree with physical interpretation and demonstrate the capability of physics-based microstructure sim- ulations for EIS behavior. 20 1.4 Dissertation Outline In this thesis work, a model of electrochemical double layer formation is developed. The calculated capacitance is input to microstructure-level electrochemical simulations with os- cillating loadings. This simulation framework is applied to investigate EIS behavior of dif- ferent electrodes under various states. The main content of this thesis is briefly summarized below. In Chapter 2, the methodologies of Smoothed Boundary Method (SBM), Finite Difference Method (FDM), and Adaptive Mesh Refinement (AMR), in the simulations are introduced. The SBM is used to reformulate governing equations sucha the new equations can be solved on non-body conforming meshes. FDM is employed to approximate differential operators. AMR is utilized to create fine grids near the diffuse interface regions These methods are developed for solving governing differential equations in irregular and complex domains. In Chapter 3, electrochemical double layer formation is simulated using the methods mentioned above. The Nernst-Planck-Poisson equations are reformulated with SBM and are solved on tree-structured AMR grid systems. One dimensional simulations are conducted to verify the correctness and accuracy of the methodologies. Two dimensional simulations show how complex geometry surfaces affect the formation of double layers. Furthermore, a three dimensional simulation demonstrates the capability of the presented framework. In multi-dimensions, the simulation results exhibit a dynamic process of double layer formation. The charge separation starts on high curve surfaces and then spread throughout the entire particle surfaces. In Chapter 4, we describe the multi-physics electrochemical model for a half-cell config- uration, which couples the Li diffusion in the cathode particles, current continuity in the particles, ionic diffusion in the electrolyte, current continuity in the electrolyte, and the Butler-Volmer reaction on the particle surfaces. The double-layer capacitance is calculated using the methods in Chapter 2. This simulation framework is employed to investigate the EIS behavior of NMC cathode electrodes. We examine the effects of initial Li fraction 21 (state of charge) in the cathode, average salt concentration of electrolyte, and microstruc- tures of electrodes on the resulting EIS behavior. The results well reflect the connection between intrinsic material properties, microstructure, and the EIS behavior. For instance, the charge-transfer resistance is inversely proportional to the exchange current density at respective state of charge. The salt concentration in the electrolyte dominates the the ex- change current density. At the same state of charge, the resistance is inversely proportional to the active particle surface areas. In Chapter 5, we simulate the EIS curves of complex graphite electrode. Graphite is a multi-phase material and is widely used as the anode in Li-ion batteries. The Cahn-Hilliard phase field equation is used to model the phase transformation processes in graphite particles upon lithiation/delithiation. In single-phase stages, the obtained charge-transfer resistances reflect the total active surface areas as in typical solid-solution materials. In contrast, in two-phase coexisting graphite electrodes, when phase boundaries are present on the particle surfaces, the simulations exhibit an inductive loop on the EIS curve. Yet, in core-shell phase-morphology cases, the EIS results reflect only the properties of the shells. In Chapter 6, we combine a NMC cathode and a graphite anode to a full-cell configuration to simulate the EIS of a Li-ion battery at the microstructure level. The EIS curves of the full cell at different states of charge are extracted and analyzed. The results demonstrate that our simulations correctly reflect the input material property parameters. However, due to the uncertainty in material properties, the simulated EIS curves do not well agree with reported data. As result, the presented tool is utilized to back estimate the input material properties, which demonstrate the capability of our model in calibrating material properties. Chapter 7 summarizes this thesis work presented in Chapter 2–6. The key findings and challenging issues in the investigations of physical-based electrochemical simulation of lithium-ion batteries are described. Lastly, we provide an outlook for future research in the relevant field. 22 CHAPTER 2 NUMERICAL METHOD AND TECHNIQUES FOR SOLVING DIFFERENTIAL EQUATIONS IN SIMULATIONS 23 2.1 Introduction Numerical methods and techniques used in this dissertation work in simulating the electro- chemical processes in battery electrode microstructures are presented in this chapter. This project aims to develop a framework with a broad application involving complex microstruc- tures in the materials science field. In section 2.2, the Smoothed Boundary Method (SBM) is introduced to reformulate differential equations, such that the new equations can be solved on mesh/grid systems non-conformal to complex geometries. This technique is particularly advantageous in simulations involving real electrode particles with irregular shapes. Finite difference method (FDM) is used to implement the differential operators in solving the SBM- forumlated equations. Basic stencils for a uniform grid are described in section 2.3. With a diffuse interface of non-zero thickness in the SBM, errors between the SBM and the con- ventional sharp-interface results are inevitably introduced. Thus, tree-structured Adaptive Mesh Refinement (AMR) is employed to allow the use of thinner diffuse interfaces, which significantly reduces the modeling error of SBM simulations. Using AMR grids also greatly saves computational cost, compared to simulations with uniform fine-grid systems. The details of tree-refinement methods and associated FDM stencils are provided in section 2.4. 24 2.2 Smoothed Boundary Method Smoothed boundary method has been demonstrated to be a convenient tool for solving partial differential equations with boundary conditions imposed on irregular boundaries [70]. In the SBM, the complex geometries are embedded in the computational box, and a regular grid system is generated by discretizing that computational box. The immersed boundary is obtained by smearing the original sharp interface to a diffuse interface with four-to-six grid points across. Since the interface is no longer ‘sharp’, structural mesh conforming with geometries are not required. This circumvents the tedious processes of generating conformal meshes on complex geometries. In the conventional sharp-interface modeling, which are usually implemented by finite element method (FEM) or finite volume method (FVM), it is required to discretize the do- main with mesh (elements) conformal to the geometry of the domain. The equations are solved within the meshed domain. For instance, the FEM subdivides the domain into smaller and simpler parts, such as triangles (in 2D) and tetrahedrons (in 3D). These subdivisions are called the elements, and the total union of the elements is referred to mesh. To accu- rately represent a complex geometry, a high quality conformal mesh is required, especially in the regions with high curvatures. The meshing process is very time-consuming and te- dious. In many cases, manually creating or modifying elements are necessary. Furthermore, an individual mesh system must be generated for each complex geometry. To investigate microstructure effects using a series of different microstructures, the total efforts for gener- ating many conformal mesh systems will be extremely high. It is acknowledged that there have been many complex microstructure simulations performed on successfully meshed ge- ometries. However, the efforts for generating conformal meshes is still a main challenge in microstructure simulations. Figure 2.1 shows a typical finite element mesh in lithium-ion battery simulations [71]. From left to right, it consists of a lithium metal anode (gray), a solid electrolyte (green), and the cathode active material particles (anthracite). The dimension of composite cathode 25 Figure 2.1 Finite element mesh of a simplified solid state lithium-ion battery electrode. The active particles have a spherical shape. is about ten micrometers. The spatial discretization is based on tetrahedral elements using Coreform Cubit, and it consists of 2,430 nodes for the mesh as a reduced cluster for a realistic geometry. A representative electrode geometry will result in finite element mesh consisting of over 1.3 million tetrahedral elements. Figure 2.2 shows finite volume mesh of a three- Figure 2.2 A segment of reconstructed cathode with cut-cell mesh structure on the top face, taken from [8] dimensional cathode microstructure geometry generated by ANSYS TGRID software [8]. 26 Cartesian finite volume meshes are produced within the interiors of electrode particles. It requires approximately 1.5 million volume cells and 5 million cell faces to render this electrode volume. Despite highly automated algorithms for generating cells for finite volume method, there are sometimes pathologically skewed cells that need to be manually repaired. Such a manual repairing procedure is commonly faced in electrode microstructure simulations [63]. While the meshing process itself might be computationally modest, the repairing process is usually very tedious, which poises one of the main challenge in electrode microstructure level simulations. On the contrary, SBM uses a continuous domain parameter, ψ, to describe the domain of interest (denoted as Ω), where the differential equations are solved. Within the domain Ω, ψ is uniformly one. ψ(x, y, z) = 1 ∈ Ω. (2.1) It continuously transitions through the domain boundary and its value becomes uniformly zero outside of the domain Ω. ψ(x, y, z) = 0 ∈ / Ω. (2.2) Since the domain parameter is continuous, the domain boundary becomes a diffuse interface with a finite thickness. Thus, the region where 0 < ψ(x, y, z) < 1 ∈ ∂Ω (2.3) implicitly defines the location of the interface. The original boundary now becomes an embedded interface within the computational box. There are various choices of functions that can be used as the domain parameter ψ, as long as they are Sigmoid functions. They need to be continuous, differentiable, and monotonically transitioning through the diffuse interface. In this work, a hyperbolic tangent function is selected to be used as the domain parameter:    1 d ψ = 1 + tanh (2.4) 2 ζ 27 where d is a signed distance function, which describes the the shortest distance from a point to the original boundary. A positive distance indicates grid points inside the domain. If a grid point is outside the domain, the distance is defined as negative. ζ is a parameter that can control the thickness of the diffuse interface. Figure 2.3 (a) Domain parameter distribution in the computational domain, (b) modified hyperbolic tangent function profile ψ versus distance to the boundary. The original differential equations will be reformulated to include the their boundary conditions. For example, the mass conservation equation is expressed as ∂C = −∇ · ⃗j + S (2.5) ∂t where C is the concentration, t is time, ⃗j is the flux vector, and S is the source term. Multiplying both sides by ψ, and using the identity ψ∇ · ⃗j = ∇ · (ψ⃗j) − ∇ψ · ⃗j, we obtain the SBM-formulated equation as ∂C ψ = −[∇ · (ψ⃗j) − ∇ψ · ⃗j] + ψS. (2.6) ∂t Using diffusion as an example for the formulation, the mass flux is defined by Fick’s law as ⃗j = −D∇C, where D is the diffusion coefficient. Equation (2.6) becomes ∂C ψ = [∇ · (ψD∇C) − ∇ψ · D∇C] + ψS. (2.7) ∂t 28 The flux across the domain boundary is referred to as the Neumann boundary condition, which specifies the gradient normal to the boundary. ∂C jB = ⃗n · ⃗j = −D , (2.8) ∂n where ⃗n in the inward unit normal vector to the boundary. The inward normal vector is defined as ⃗n = ∇ψ/|∇ψ| in the diffuse interface description. Thus, we can obtain the relation: ∇ψ ∂C ⃗n · (−D∇c) = · (−D∇c) = −D =⇒ ∇ψ · D∇c = −|∇ψ|jB . (2.9) |∇ψ| ∂n Using Eq. (2.9), the SBM-formulated diffusion equation is ∂C   ψ = ∇ · (ψD∇C) + |∇ψ|jB + ψS. (2.10) ∂t The second term on the right side serves as the internal boundary condition imposed on the diffuse interface. The computational box is discretized with a regular grid system that does require grid points to be on the embedded boundary. The SBM-reformulated equation is solved in the entire computational box with the embedded boundary conditions. The SBM formulation is applicable for other types of mass transport equations by replacing different definitions of mass flux ⃗j in Eq. (2.6), e.g., using the Cahn-Hilliard equation of phase transformation. For a closed system of the complex geometry (no-flux boundary condition on the diffuse interface), the second term in Eq. (2.6) vanishes. The domain parameter can be employed to couple the governing equations in connect domains. For instance, in the electrochemical simulations, ψp = 1 is used to define the particle regions, and ψe = 1 − ψp is used to defined the electrolyte regions. Further details of SBM-formulated equations for the simulations are provided in the respective chapters later. 29 2.3 Finite Difference Method on Uniform Grid Systems SBM formulated differential equations can be solved on regular grid systems with the complex geometries embedded in the computational box. For a uniform, regular grid, finite difference method (FDM) is the most efficient numerical technique to implement. Although FDM on uniform grids is baraly used through this work, the basics of FDM is still provided here as a necessary background knowledge for the derivation of FDM stencils on adaptively refined grid systems later. The central difference scheme approximates a second order 1D derivative by 1 u′′ (x) ≈ [u(x + h) − 2u(x) + u(x − h)], (2.11) h2 where h is the grid spacing. Using Taylor series 1 1 u(x + h) = u(x) + hu′ (x) + h2 u′′ (x) + h3 u′′′ (x) + O(h4 ) (2.12) 2 6 1 1 u(x − h) = u(x) − hu′ (x) + h2 u′′ (x) − h3 u′′′ (x) + O(h4 ), (2.13) 2 6 The truncation error can be shown to be proportional to h2 . The first order 1D derivative is approximated by u(x + h) − u(x − h) u′ (x) ≈ (2.14) 2h with a truncation error in the order of h2 . In 2D, the Laplace operator is expressed as u(x + h, y) + u(x − h, y) + u(x, y − h) + u(x, y + h) − 4u(x, y) ∇2 u ≈ (2.15) h2 and the gradients are ∂u u(x + h, y) − u(x − h, y) ≈ (2.16) ∂x 2h ∂u u(x, y + h) − u(x, y − h) ≈ . (2.17) ∂y 2h The truncation error of the gradient calculated using central different method is also pro- portional to h2 . The FDM stencils for 3D cases can be derived similarly. 30 2.4 Adaptive Mesh Refinement In this thesis, quadtree and octree decomposition are utilized to create adaptive meshes for the 2D and 3D simulations, respectively. The adaptive mesh refinement is used to decrease the computational burden and to increase the model accuracy. As introduced in Chapter 1, the diffuse double layer spans only tens of nanometers near the electrode surfaces. In contrast, the dimensions of electrode particles and electrolyte space span over several (or tens of) micrometers. Using a uniformly fine grid that is able to resolve the ionic concentration and electropotential gradients within the double layer will lead to an enormous number of grid points for the entire system. Thus, we must use fine mesh in the region of double layer and use coarse mesh in the bulk of particles and electrolyte. Moreover, the SBM introduces an interface with a non-zero thickness. This treatment inevitably leads to deviations [70] in the results obtained using SBM and the conventional sharp-interface methods. This modeling error quickly diminishes if the interface is sufficiently thin. Generally four to six grid spacings are required to span over the diffuse interface to ensure the gradients can be properly calculated within the interface. Fine grid near the diffuse interface region allows the use of domain parameter that results in a thin interface, which can significantly reduce the modeling error. Note that even though AMR grid is used in the simulations, the grid points are not required to locate on the original sharp interface. In the following, we briefly describe the quadtree refinement and the data structure associated with the refinement. Two nouns are defined here: cell and node. In 2D cases, a cell is a square formed by four equal edges. For example, cell-a and cell-c in Fig. 2.4(a). A node is the intersect of different edges and it is also referred as a grid point. For instance, node-3 and node-7 in Fig. 2.4(a). The 2D domain is initially discretized into a regular Cartesian grid system, resulting in all equal-sized square cells. The resulting discretization is called the Root Level of the tree. A Cell-List is created, in which each row contains the labels of the four corner nodes in clockwise order. See Fig. 2.4(a): for example, cell-a with node-1 through node-4 are stored in the first row of the cell list. If the distance from a 31 Figure 2.4 An example of cell dividing and node labeling (a) cell and node label before dividing (b) cell and node label after cell-b was divided cell’s center to the embedded boundary is less than a threshold value, that cell is split to four equal-sized squares, and the original cell is eliminated. Consequently, four new cells are created and new nodes are inserted into the Cell-List. See Fig. 2.4(b) for instance, cell-b is eliminated and cell-e through cell-h are created which are referred to as the Level-1 cells. In the meantime, node-10 through 14 are inserted as shown in Fig. 2.4(b). The Level-1 cells can be further quadtree-decomposed into Level-2 cells, so on and so forth to higher levels of refinement. Note here special care must be taken to set the threshold distances such that there will be only one level of refinement between adjacent cells. A cross-level refinement will require a more complicated FDM stencils. Thus, for the simplicity of FDM 32 implementation, we constrain the refinements to ensure that all cells are adjacent to cells that are only one-level higher or lower in the refinements. Figure 2.5 An example Node-Neighbor-List of Fig. 2.4(b), regular node 3 and 11 have four neighbors; T-junction node 12 and 14 have five neighbors, two west neighbors and two south neighbors respectively. During the refinement process, two additional files are generated: (1) Node-List that stores the position coordinate of each node, and (2) Node-Neighbor-List that stores the labels of the neighboring nodes of each node along the coordinate axial directions. For instance, see Fig. 2.5, in which node-3 has four direct neighboring nodes in the four axial directions. However, due to the refinement, some nodes do not have direct neighbors along axial directions. For example, node-12 has two indirect neighbors in the west direction. As shown in Fig. 2.6(a), a typical node (e.g., u0 ) will have four neighbor nodes on its north, south, west and east directions, respectively. When a node inserted on the boundary of two different levels of refinement, it can miss some of its direct neighbor nodes, e.g., see u0 in Fig. 2.6(b). The u0 does not have a direct west neighbor node but it has two indirect neighbor nodes: u5 and u6 . The indirect neighbor nodes will also be stored in the Node-Neighbor-List (mentioned above) for calculating the derivatives using FDM. Instead of storing the values of u in a 2D array as in a typical FDM simulation on a uniform grid system, the values of u in this work are stored in a 1D vector, in which the indices are the node labels. Derivatives can be calculated by using the values at the neighboring nodes ui 33 according to the node labels stored in the Node-Neighbor-List. The FDM stencils for AMR grid systems will be introduced in the next paragraph. Figure 2.6 (a) Regular cross node Node u0 has four direct neighbor nodes, u1 , u2 , u3 , and u4 in the west, east, south and north directions respectively. (b) West facing T-junction, node u0 has three direct neighbor nodes, u2 , u3 , and u4 ; and two indirect neighbor nodes, u5 and u6 . As mentioned earlier, the finite difference method is used to solve governing differential equations. It is a node-based sampling method, meaning values are calculated at nodes (or grid points), instead of at the cell centers. In 2D cases, there are two types of nodes. One with four direct neighbor nodes as shown in Fig. 2.6(a). This type of nodes is referred to the regular nodes. The other type of nodes has only three or less direct neighbors. For instance, u0 on Fig. 2.6(b) has u2 , u3 , and u4 as the east, south, and north neighbors, respectively. There is no direct neighbor node on its west (u1 is a virtual node) but two indirect neighbor nodes u5 and u6 instead. This type of nodes is referred to the T-junction nodes, specifically west facing T-junction node in this case. Clearly, the difference in neighboring relations, compared to the regular node, will lead to a different FDM stencil for the T-junction node. Below, an example derivation of the FDM stencil for AMR grid system is presented. 34 We will begin with the Laplace operator ∇2 u. As in Fig. 2.6(a), for regular node u0 with four neighbor nodes, the Taylor series gives s21 u1 = u0 − s1 ux + uxx + ... (2.18) 2 s22 u2 = u0 + s2 ux + uxx + ... (2.19) 2 s2 u3 = u0 − s3 uy + 3 uyy + ... (2.20) 2 s2 u4 = u0 + s4 uy + 4 uyy + ... (2.21) 2 where the subscripts x and y denote differentiation. Eq. (2.18)×s2 + Eq. (2.19)×s1 leads to s1 s2 (s1 + s2 ) s1 (u2 − u0 ) + s2 (u1 − u0 ) = uxx (2.22) 2 which can be reorganized to   2 u2 − u0 u0 − u1 uxx = − . (2.23) s1 + s2 s2 s1 Eq. (2.20)×s4 + Eq. (2.21)×s3 gives   2 u4 − u0 u0 − u3 uyy = − . (2.24) s3 + s4 s4 s3 Finally, Eq. (2.23) + Eq. (2.24) results in the stencil for Laplace operator     2 2 u2 − u0 u0 − u1 2 u4 − u0 u0 − u3 ∇ u = uxx + uyy = − + − . (2.25) s1 + s2 s2 s1 s3 + s4 s4 s3 With the domain parameter function ψ as in the SBM equations, we have   2 ψ2 + ψ0 u2 − u0 ψ0 + ψ1 u0 − u1 ∇ · ψ∇u = · − · + s1 + s2 2 s2 2 s1   2 ψ4 + ψ0 u4 − u0 ψ0 + ψ3 u0 − u3 · − · s3 + s4 2 s4 2 s3     2 u2 − u0 u0 − u1 2 u4 − u0 u0 − u3 = ξ2 · − ξ1 · + ξ4 · − ξ3 · s1 + s2 s2 s1 s3 + s4 s4 s3 (2.26) where ξi = (ψi + ψ0 )/2 is the average of ψ between the center node and its individual neighbors, and i = 1, 2, 3 and 4. 35 Here, we use a west facing T-junction node as an example for deriving the FDM stencil for Laplace operator. For the T-junction node u0 on Fig. 2.6(b), the Taylor series of its indirect neighbor nodes u5 and u6 can be written as s21 s25 u5 = u0 − s1 ux − s5 uy + uxx + uyy + s1 s5 uxy + ... (2.27) 2 2 s2 s26 u6 = u0 − s1 ux + s6 uy + 1 uxx + uyy − s1 s6 uxy + ... (2.28) 2 2 Neglecting high order terms (after second order derivatives), Eq. (2.27)×s6 + Eq. (2.28)×s5 yields  2  s1 s5 s6 s6 u5 + s5 u6 = (s5 + s6 )(u0 − s1 ux ) + (s5 + s6 ) uxx + uyy 2 2 1  s2 s5 s6 s5 (u0 − u6 ) + s6 (u0 − u5 ) = s1 ux − 1 uxx −  =⇒ uyy s5 + s6 2 2 (2.29) s5 u6 + s6 u5 s21 s5 s6 =⇒ u0 − = s1 ux − uxx − uyy s5 + s6 2 2 s2 s5 s6 =⇒ u0 − u1 = s1 ux − 1 uxx − uyy 2 2 where u1 represents the value on a virtual node between u5 and u6 , and its value is an weighted average of u5 and u6 . Eq. (2.19)×s1 − Eq. (2.29)×s2 eliminates ux , yielding s1 s2 (s1 + s2 ) s2 s5 s6 s1 (u2 − u0 ) − s2 (u0 − u1 ) = uxx + uyy 2 2 (2.30) 2 u2 − u0 u0 − u1 s5 s6 =⇒ uxx = ( − )− uyy s1 + s2 s2 s1 s1 (s1 + s2 ) On the y direction, node u0 has two direct neighbors. Thus, uyy has the same form as the regular node as expressed in Eq. (2.24). Summing Eq. (2.30) and Eq. (2.24) gives uxx + uyy =       (2.31) 2 u2 − u0 u0 − u1 2 u4 − u0 u0 − u3 s5 s6 − + − · 1− . s1 + s2 s2 s1 s3 + s4 s4 s3 s1 (s1 + s2 ) This stencil is similar to that of a regular node, except for the correction factor β = s5 s6 /s1 (s1 + s2 ) that arises due to the lack of a direct neighbor. With the domain parameter function ψ, we start with second derivative in the y direction:   2 ψ4 + ψ0 u4 − u0 ψ0 + ψ3 u0 − u3 (ψuy )y = · − · s3 + s4 2 s4 2 s3   (2.32) 2 u4 − u0 u0 − u3 = ξ4 · − ξ3 · . s3 + s4 s4 s3 36 On the x direction, we use the two indirect neighbor nodes to represent the virtual west neighbor node, which yields  2 ψ2 + ψ0 u2 − u0 (ψux )x = · − s1 + s2 2 s2   s5 ψ0 + ψ6 u0 − u6 s6 ψ0 + ψ5 u0 − u5 s5 s5 · · + · · − (ψuy )y s5 + s6 2 s1 s5 + s6 2 s1 s1 (s1 + s2 )    2 u2 − u0 s5 u0 − u6 s6 u0 − u5 = ξ2 − · ξ6 · + · ξ5 · − β(ψuy )y . s1 + s2 s2 s5 + s6 s1 s5 + s6 s1 (2.33) Again, due to the lack of a direct neighbor, the correction factor β appears in the stencil. We can define   s5 ξ6 + s6 ξ5 1 ψ0 + ψ6 ψ0 + ψ5 ξ1 = = s5 · + s6 · s5 + s6 s5 + s6 2 2   (2.34) ψ0 1 s5 ψ6 s6 ψ5 1  = + + = ψ0 + ψ1 2 2 s5 + s6 s5 + s6 2 where ψ1 = (s5 ψ6 + s6 ψ5 )/(s5 + s6 ) is the weighted average value of ψ on the virtual neighboring node. Substituting ψ1 into Eq. (2.33), we obtain     2 u2 − u0 1 s5 ξ6 s6 ξ5 (ψux )x = ξ2 · − ξ1 · u0 − · · u6 + · · u5 − β(ψuy )y s1 + s2 s2 s1 s5 + s6 ξ1 s5 + s6 ξ1   2 u2 − u0 u0 − u1  = ξ2 · − ξ1 · − β ψuy y s1 + s2 s2 s1 (2.35) where the value at the virtual neighbor node is obtained by a weighted average: s5 ξ6 s6 ξ5 u1 = · · u6 + · · u5 . (2.36) s5 + s6 ξ1 s5 + s6 ξ1 Note that if ψ is a uniform value, u1 equals that in Eq. (2.29). Combining Eq. (2.32) and Eq. (2.35), we obtain     2 u2 − u0 u0 − u1 2 u4 − u0 u0 − u3  ∇ · ψ∇u = ξ2 · − ξ1 · + ξ4 · − ξ3 · · 1−β . s1 + s2 s2 s1 s3 + s4 s4 s3 (2.37) 37 Stencils for T-junction nodes in other directions can be derived in a similar manner. A general form of the Laplace operator with domain parameter on the AMR grid system is   2 u2 − u0 u0 − u1 (∇ · ψ∇u)T = ξ2 · − ξ1 · · (1 − β1 )+ s1 + s2 s2 s1   (2.38) 2 u4 − u0 u0 − u3 ξ4 · − ξ3 · · (1 − β2 ) s3 + s4 s4 s3 where β1 and β2 are the correction factors for nodes missing direct neighbor nodes in the x and y directions, respectively. For a west facing T-junction, β2 = s5 s6 /(s1 (s1 + s2 )), β1 = 0, u1 = (s6 u5 + s5 u6 )/(s5 + s6 ) which represents the value at a virtual neighbor on the west side. ξ1 = (ψ1 + ξ0 )/2, where ψ1 = (s6 ψ5 + s5 ψ6 )/(s5 + s6 ) is the average of ψ at the two indirect neighbors. It is noted that in a given direction, a T-junction node may have immediate indirect neighbors located at different distances away from the node. Consider node u0 in Fig. 2.6(b) as an example, where the immediate northwest neighbor (u7 ) is closer than the immediate southwest neighbor (u5 ) in that direction. In a such case, instead of using the closest neighbor, the next neighbor in that direction is chosen in order to ensure that distance from a T-junction to all its neighbors in a given direction is the same. This means that for a node like u0 on Fig. 2.6(b), u6 will be chosen as the northwest neighbor rather than u7 . For an east-facing T-junction, β2 = s5 s6 /(s2 (s1 + s2 )), β1 = 0, and u2 and ψ2 are the respective average values of u and ψ from the two indirect neighbors. For south or north- facing T-junctions, β2 = 0 and β1 is nonzero to account for the fact that a direct neighbor in the south or north direction is missing. The ui and ψi values at their virtual neighbors are obtained by respective average from the two indirect neighbors. The first order derivatives (gradient) of regular nodes and T-junction nodes will be derived as the following. For a regular node, u2 − u0 s1 u0 − u1 s2 ux = · + · (2.39) s2 s1 + s2 s1 s1 + s2 u4 − u0 s3 u0 − u3 s4 uy = · + · . (2.40) s4 s3 + s4 s3 s3 + s4 38 Note that since the distances between u0 and its neighbor are not always the same, the accuracy of approximation for those points is no long second order. For the T-junction node in Fig. 2.6(b), we have s6 u 5 + s5 u 6 u1 = = s5 + s6 s21   1 s5 s6 (s6 + s5 )u0 − s1 (s6 + s5 )ux + (s6 + s5 ) uxx + (s6 + s5 ) uyy (2.41) s5 + s6 2 2 s2 s5 s6 = u0 + s1 ux + 1 uxx + uyy 2 2 which leads to s21 s5 s6 u1 = u0 − s1 ux + uxx + uyy . (2.42) 2 2 On the other hand, for node u2 , its Taylor series is given as Eq. (2.19). Eq. (2.42)×s22 + Eq. (2.19)×s21 results in s21 s5 s6 s21 (u2 − u0 ) + s22 (u0 − u1 ) = (s22 s1 + s21 s2 )ux − uyy (2.43) 2 Reorganizing the equation, it is obtained that u2 − u0 s1 u0 − u1 s2 s1 s5 s6 ux = · + · + uyy . (2.44) s2 s1 + s2 s1 s1 + s2 2s1 (s1 + s2 ) The last term represents a correction factor due to the lack of a direct neighbor. Generalized expressions of FDM stencils for the gradient operator ∇u can be written as ∂u u2 − u0 s1 u0 − u1 s2 = · + · + α2 uyy (2.45) ∂x s2 s1 + s2 s1 s1 + s2 ∂u u4 − u0 s3 u0 − u3 s4 = · + · + α1 uxx . (2.46) ∂y s4 s3 + s4 s3 s3 + s4 The correction factors, α1 and α2 , account for the effect of missing direct neighbor nodes in the y and x direction, respectively. For a regular node, α1 and α2 are zero. For a west-facing T-junction as in Fig. 2.6(b), α1 = 0 and α2 = s5 s6 s2 /(2s1 (s1 + s2 )). For an east- facing T-junction, α1 = 0 and α2 = −s5 s6 s1 /(2s2 (s1 + s2 )). Similarly, the correction factor α1 = s5 s6 s4 /(2s3 (s3 + s4 )) for a south-facing T-junction and α1 = −s5 s6 s3 /(2s4 (s3 + s4 )) 39 node type α1 α2 β1 β2 ss sn se ss sn Tw 0 2sw (sw +se ) 0 sw (sw +se ) Te 0 − 2sess(sswn s+s w e) 0 ss sw se (sw +sw ) sw se sn ss se Ts 2ss (ss +sn ) 0 ss (ss +sn ) 0 ss se Tn − 2snsw(ssse+s ss n) 0 sn (ss +sn ) 0 Table 2.1 Correction factors used in the finite difference stencils for 2D quadtree AMR. α1 and α2 are the coefficients that appear in Eqs. (2.46) and (2.45). β1 and β2 are the coefficients that appear in Eq. (2.38). The subscripts w, e, s, and n in T and s denote west, east, south, and north, respectively. For example, sw , se , ss , and sn correspond to the notations s1 , s2 , s3 , and s4 , respectively, for the west-facing T-junction (Tw ) at node 0 in Fig. 2.6(b). For a regular node, all correction factors are zero. for a north-facing T-junction, while α2 = 0. Table 2.1 lists the correction factors of gradient operator and Laplace operator for different facing T-junctions in 2D cases. In 3D cases, Octree refinement is performed, in which cells are cubes instead of squares in 2D. Cubic cells will be split into eight equal-sized cubes. The Cell-List contains the labels of eight corner nodes. Similar procedures as in the quadtree refinement are conducted to generate nodes and identify the neighboring nodes. A cell list is created to store the cell labels and their eight vertex nodes in a consistent order. A node list is also created to store the labels and positions of the nodes. If the distance from the center of a cell to the embedded boundary is less than a threshold value, that cell will be divided into eight new equal-sized cubes. The original cell is eliminated from the cell list and the newly generated eight cells are inserted into the cell list. The new cells are referred to as the Level-1 of the tree. New nodes are added to the node list. This process can be continued to multiple levels of refinement. Once the refinement is completed, the neighboring nodes of each node can be determined according to the node positions and stored the labels to a node neighbor list. The FDM stencils for 3D differential operators can be derived via similar procedures as in 40 the 2D cases presented earlier, with a general form as   2 u2 − u0 u0 − u1  ∇ · ψ∇u = ξ2 · − ξ1 · · 1 − β21 − β31 + s1 + s2 s2 s1   2 u4 − u0 u0 − u3  ξ4 · − ξ3 · · 1 − β12 − β32 + (2.47) s4 + s3 s4 s3   2 u6 − u0 u0 − u5  ξ5 · − ξ5 · · 1 − β13 − β23 . s6 + s5 s6 s5 Note here, in the 3D case, u5 and u6 are used to denote the neighbors (direct or virtual (a) (b) (c) Figure 2.7 Examples of configurations of T-junctions: (a) west-facing T on the x-z plane, (b) a node that is simultaneously a west-facing T on the x-y plane and a bottom-facing T on the y-x plane, and (c) west-facing face-centered T, which has four indirect neighbors in the west direction. The black, magenta, and cyan dots indicate the center node, direct neighbors, and indirect neighbors. The green circles indicate virtual neighbors. neighbors) in the down (−z) and up (+z) directions, respectively. ξi = (ψi + ψ0 )/2 is the average value of ψ between the center node and its neighbor. For regular nodes having 6 41 direct neighbors in 3D, all correction factors β vanish and Eq. (2.47) reduces back to a typical 7-point 3D FDM stencil. For a T-junction missing a direct neighbor in the x direction, β12 and β13 are the correction factors calculated using the indirect neighbors on the x-y plane and x-z plane, respectively. Their values are listed in the Table 2.2. Figure 2.7(a) shows an example of a west-facing T-junction on the x-z plane. In this case, the center node has two indirect neighbors in the bottom-west and top-west directions. The correction factor of β13 is calculated according to β13 = sB sT /(sW (sW + sE )) and zero for all other βs. The subscript ‘13’ indicates this T-junction is along the x-axis and on the x-z plane. The value of ψW at the virtual west neighbor node can be calculated by averaging those at the two indirect neighbors in the west direction as ψW = (ψBW + ψT W )/2, where ψBW and ψT W are the values of ψ in the two indirect neighbors on the bottom-west and top-west directions, respectively. The value of uW (on the virtual west neighbor indicated by the green circle in Fig. 2.7(a)) is calculated by a weighted average of those at the two indirect neighbors as (ψBW + ψC )uBW + (ψT W + ψC )uT W uW = . (2.48) (ψBW + ψC ) + (ψT W + ψC ) It is possible that a node can be a T-junction on two orthogonal planes in 3D, as shown in Fig. 2.7(b). In this case, the values of u at the two virtual neighbors can be calculated separately using a similar method as in Eq. (2.48). If a node is simultaneously a T-junction on the x-y and y-z planes, β12 and β23 are nonzero. Furthermore, if a node has four indirect neighbors, as shown in Fig. 2.7(c) for an example of west-facing face-centered T-junction, the values at the virtual neighbor is obtained by a weighted average according to uW = P P i (ψi + ψC )ui / i (ψi + ψC ). The formulae of βs are provided in Table 2.2. Similar to the 2D case, the generalized FDM stencils for 3D gradient operator can also be derived, but they are not repeated here. Details can be found in Ref. [72, 73]. For the time derivatives, the fully forward Euler explicit method and fully backward Euler implicitly are used for different transport equations in this thesis. For a uniform grid system, parallel computing can be straightforwardly implemented by the technique of domain decomposition. The entire computational box is divided into several 42 sections, and each section is allocated to one CPU for the calculations within the subdivision. Figure 2.8 illustrates a domain decomposition along the x-axis. Communication between dif- ferent CPUs is conducted using MPI (Message Passing Interface), as calculating derivatives for the boundary points of one section will need the information from the neighboring sec- tions. On AMR grid systems, decomposition can be performed similarly. As the grid points are not aligned as in the uniform grid case, an additionally procedure is conducted to detect the labels of the neighbor nodes across the dividing boundary. These labels are stored such that the node values on those points can be sent to and received by the neighboring sections by MPI send and MPI receive functions during the computation . Figure 2.8 Illustration of domain dividing and section assigning to different ranks 43 2.5 Conclusions In this chapter, SBM is introduced to handle complex geometries such that conformal mesh is not required in relevant simulations. A general SBM formulation for imposing Neumann boundary conditions on the diffuse interface is demonstrated. Finite difference method is used in this work to approximate derivatives in solving differential equations. To reduce the computational burden and increase the modeling accuracy tree-structure adaptive mesh refinement is employed such that fine mesh is used near the diffuse interface but coarse mesh is used in the bulk phase. The relevant FDM stencils for AMR is also presented. While those FDM stencils have been previously derived for general non-uniform grid systems [74, 75], for graded tree-structure AMR systems, we organize those 2D and 3D stencils into a generalized form, which is very similar to the standard uniform grid FDM stencils. This greatly simplifies the development of simulation code. SBM introduces a diffuse interface with a non-zero thickness, which deivates the conventional sharp-interface modeling. Such a modeling error in double layer formation is examined in the next chapter. 44 Table 2.2 Correction factors for the finite difference stencils used in this work. Adapted from Ref. [72]. Note that αij (T∗ ) and βij (T∗ ) refer to the αij (T∗ ) and βij (T∗ ) value from T-junction type of T∗ . node type α12 α13 β12 β13 ss sn se ss sn Tw−xy 2sw (sw +se ) 0 sw (sw +se ) 0 sb st se sb st Tw−xz 0 2sw (sw +se ) 0 sw (sw +se ) Tw−xy−xz 0.5α12 (Tw−xy ) 0.5α13 (Tw−xz ) 0.5β12 (Tw−xy ) 0.5β13 (Tw−xz ) Tw−4i α12 (Tw−xy ) α13 (Tw−xz ) β12 (Tw−xy ) β13 (Tw−xz ) ss sn sw ss sn Te−xy - 2se (sw +se ) 0 se (sw +se ) 0 sb st sw sb st Te−xz 0 - 2se (sw +se ) 0 se (sw +se ) Te−xy−xz 0.5α12 (Te−xy ) 0.5α13 (Te−xz ) 0.5β12 (Te−xy ) 0.5β13 (Te−xz ) Te−4i α12 (Te−xy ) α13 (Te−xz ) β12 (Te−xy ) β13 (Te−xz ) node type α21 α23 β21 β23 sw se sn sw se Ts−xy 2ss (ss +sn ) 0 ss (ss +sn ) 0 sb st sn sb st Ts−yz 0 2ss (ss +sn ) 0 ss (ss +sn ) Ts−xy−yz 0.5α21 (Ts−xy ) 0.5α23 (Ts−yz ) 0.5β21 (Ts−xy ) 0.5β23 (Ts−yz ) Ts−4i α21 (Ts−xy ) α23 (Ts−yz ) β21 (Ts−xy ) β23 (Ts−yz ) sw se ss sw se Tn−xy - 2sn (ss +sn ) 0 sn (ss +sn ) 0 sb st ss sb st Tn−yz 0 - 2sn (ss +sn ) 0 sn (ss +sn ) Tn−xy−yz 0.5α21 (Tn−xy ) 0.5α23 (Tn−yz ) 0.5β21 (Tn−xy ) 0.5β23 (Tn−yz ) Tn−4i α21 (Tn−xy ) α23 (Tn−yz ) β21 (Tn−xy ) β23 (Tn−yz ) node type α31 α32 β31 β32 ss sn st ss sn Tb−yz 0 2sb (sb +st ) 0 sb (sb +st ) sw se st sw se Tb−xz 2sb (sb +st ) 0 sb (sb +st ) 0 Tb−yz−xz 0.5α31 (Tb−xz ) 0.5α32 (Tb−yz ) 0.5β31 (Tb−xz ) 0.5β32 (Tb−yz ) Tb−4i α31 (Tb−xz ) α32 (Tb−yz ) β31 (Tb−xz ) β32 (Tb−yz ) ss sn sb ss sn Tt−yz 0 - 2st (sb +st ) 0 st (sb +st ) Tt−xz - 2sstw(ssbe+s sb t) 0 sw se st (sb +st ) 0 Tt−yz−xz 0.5α31 (Tt−xz ) 0.5α32 (Tt−yz ) 0.5β31 (Tt−xz ) 0.5β32 (Tt−yz ) Tt−4i α31 (Tt−xz ) α32 (Tt−yz ) β31 (Tt−xz ) β32 (Tt−yz ) 45 CHAPTER 3 DOUBLE LAYER FORMATION ON ELECTRODES WITH COMPLEX MICROSTRUCTURE 46 3.1 Introduction The formation of double layers is the key mechanism for energy storage of electrochemical capacitors (ECs). The double layer capacitance also play a crucial role in batteries’ electro- chemical impedance spectroscopy measurements. Upon applying external voltage to the two electrodes of an EC, cations and anions in the electrolyte move to the surfaces of the cathode and anode, respectively, in order to counter the charges on the electrodes. This leads to a thin layer of charge separation in the proximity of the liquid-solid interfaces. The illustration of formation and disintegration is shown on Fig. 3.1. The same mechanism can take place in all other electrochemical systems that contain mobile ions to counterbalance polarization at electrodes; e.g., the electrochemical processes in liquid-cell batteries [76], flow batteries [77], electrochemical desalination [16], as well as the ion diffusion in biochemical systems. Figure 3.1 (a) ions move directionally under the electric field. (b) ions are attached to the electrodes accordingly. (c) electrolyte returns neutral without external electric field [78] The double-layer region usually spans 20 to 50 nanometers, whereas other significant length scales, e.g., particle size or inter-particle space in electrodes, are in tens to hundreds of microns. Thus, a direct numerical simulation in the continuum scale must resolve the ionic concentration and potential gradients in spatial scales across 3∼4 orders of magnitude 47 difference and is highly challenging. As mentioned in Chapter 1, double layer models and Nernst-Planck-Poisson equations have been developed to describe the ionic concentration and electrostatic potential evolution due to charge separation in the double layer. However, because of the challenges in handling a huge discrepancy in spatial dimensions, only 1D analytical solution or numerical solutions on extremely simple geometries were reported. In this chapter, we use the smoothed boundary method (SBM) introduced in Chapter 2 to formulate the NPP governing equations such that conformal mesh is not required for simulating the formation of electrochemical double layers on complex electrode particle ge- ometries. Specifically, the SBM formulated equations are solved on grid systems generated by AMR as described in the previous chapter. Thus, the ion concentration and electric potential variations are resolved in both the very thin diffuse layer and in the interparticle space that are several orders of magnitude larger. This approach allows the simulations of the dynamics of double layer formation in complex 2D and 3D cases, which is difficult for the conventional numerical methods. One-dimensional simulations are performed to verify the accuracy of the presented method. Two- and three-dimensional simulations are conducted to demonstrate the importance of considering the geometric effects on multidimensional cases. The results show that, depending on the charges, ions first rapidly adsorb onto or are repelled from the particle surfaces, followed by long-distance diffusion to alleviate the concentration gradient until the system reaches the steady state. The concentration and potential evo- lutions highly depend on the particle geometries. These geometric effects have not been investigated computationally before this work. 48 3.2 Model Construction and Computational Methods 3.2.1 Governing Equations Ionic diffusion and migration in electrolyte, referred to as electrodiffusion, are governed by the ionic concentration gradient and the electric field gradient, as described by the Nernst- Planck-Poisson equations [79]:   ∂c+ = −∇ · J+ = ∇ · D+ ∇c+ + z+ u+ F c+ ∇ϕ , (3.1) ∂t   ∂c− = −∇ · J− = ∇ · D− ∇c− + z− u− F c− ∇ϕ , (3.2) ∂t F ∇2 ϕ = − (c+ − c− ), (3.3) ε where ci , Di , zi , and ui denote concentration, diffusivity, charge number, and transport mobility of the ions, respectively, and Ji = −(Di ∇ci + zi ui F ci ∇ϕ) is the flux vector. The subscripts, + and −, indicate cations and anions, respectively. F is the Faraday constant, ϕ is the electrostatic potential, and ε is the dielectric constant (or relative permittivity) of the electrolyte. Here, ε is assumed to be a constant throughout the electrolyte. In this chapter, we focus only on binary electrolyte (i.e., the electrolyte contains only one species of cation and one species of anion) with monovalence ions (i.e., z+ = 1 and z− = −1). In a closed system, the boundary conditions for the two electrodiffusion equations are no- flux, such that n · Ji = 0, where n is the inward unit normal vector of the domain boundary. The Poisson’s equation is subject to a Robin boundary condition (combination of specified value and gradient) [42]: ϕ − ϕd = n · ∇ϕ, (3.4) λs where ϕd is the electropotential on the electrode surface and λs is the thickness of the Stern layer. λs is illustrated as d2, the thickness of Helmholtz layer, in Fig. 3.2. Within the outer Helmholtz plane, the gradient of electrostatic potential is assumed to be linear. Equation (3.4) ensures that the gradient of potential across the Stern layer (the left hand side) matches 49 that in the electrolyte phase (the right hand side) at the boundary (outer Helmholtz plane) of diffuse layer. Figure 3.2 Boundary condition for electrostatic potential near the double layer region [80]. As discussed earlier, the NPP equations for 1D problems have been solved by researchers. Several studies have implemented and modified the NPP model to fit different conditions, such as concentrated solution and high applied potential [40, 28, 41]. However, there are only a few 2D simulations available in the literature [41, 81, 40] and they are for regular-shaped geometries. In the conventional numerical methods with sharp interfaces, solving differential equations requires discretizing the domain of interest with conformal mesh or grid systems. Such meshing processes for complex geometries are usually time-consuming and difficult. 50 Here, we employ the SBM [70] to circumvent the tedious mesh generation process required for sharp-interface methods. A continuous domain parameter, ψ, is introduced to describe the domain of interest (i.e., the electrolyte) where the differential equations are solved. Within the electrolyte phase, ψ is uniformly one. It continuously transitions through the electrolyte-electrode interface and becomes uniformly zero in the electrode phase. There are multiple choices of mathematical of Sigmoid functions, such as hyperbolic tangent function, error function, etc., that can be used as the domain parameter ψ. In this work, we choose a hyperbolic tangent function as the ψ function (as in Ref. [70]), see Fig. 3.3. Figure 3.3 Hyperbolic tangent function used as domain parameter Assuming the electrolyte is a dilute solution, the mobility and diffusivity are related by the Einstein relation: ui = Di /RT , where R is the ideal gas constant. Multiplying Eqs. (3.1) 51 and (3.2) with ψ, and using the product rule of differentiation, we obtain     ∂c+ ψ =∇ · ψ D+ ∇c+ + z+ u+ F c+ ∇ϕ − ∇ψ · D+ ∇c+ + z+ u+ F c+ ∇ϕ ∂t   (3.5) F  =∇ · ψ D+ ∇c+ + D+ c+ ∇ϕ + ∇ψ n · J+ , RT     ∂c− ψ =∇ · ψ D− ∇c− + z− u− F c− ∇ϕ − ∇ψ · D− ∇c− + z− u− F c− ∇ϕ ∂t   (3.6) F  =∇ · ψ D− ∇c− − D− c− ∇ϕ + ∇ψ n · J− , RT where the last terms in Eqs. (3.5) and (3.6) are obtained by n · Ji = (∇ψ/|∇ψ|) · [Di ∇ci + zi ui F ci ∇ϕ] according to n = ∇ψ/|∇ψ| for the inward unit normal vector of the electrolyte- electrode interface. |∇ψ| is nonzero only at electrolyte-electrode interface and thus it indi- cates the location of the interface. The last terms in Eqs. (3.5) and (3.6) represent electro- chemical reactions occurring at the electrolyte-electrode interface and can be replaced with |∇ψ|ri , where ri is the insertion rate of cation or anions through the electrolyte-electrode interface. For a closed system (no reaction at electrolyte-electrode interface), ri = 0 and thus the last terms vanish. Similar procedures are performed to reformulate Poisson’s equation, Eq. (3.3), to the SBM form: ϕd − ϕ F ∇ · ψ∇ϕ + |∇ψ| = −ψ (c+ − c− ), (3.7) λs ε in which the second term on the left-hand side represents the Robin boundary condition for the electropotential across the electrolyte-electrode interface. Equations (3.5), (3.6), and (3.7) can be nondimensionalized with a reference diffusivity (D0 ), a reference concentration (c0 ), an electropotential unit (U = 1V), a time scale τ = l2 /D0 , and a length scale l, such that D± = D̂± D0 , c± = ĉ± c0 , ϕ = ϕ̂U , t = t̂τ , and x = x̂l, where ‘∧’ denotes dimensionless quantities. Substituting these variables into Eqs. (3.5) through (3.7), and organizing terms, we obtain   ∂ĉ+ ˆ · ψ D̂+ ∇ĉ ˆ + + A1 D̂+ ĉ+ ∇ ˆ ϕ̂ , ψ =∇ (3.8) ∂ t̂ 52   ∂ĉ− ˆ ˆ ˆ ψ = ∇ · ψ D̂− ∇ĉ− − A1 ĉ− D̂− ∇ϕ̂ , (3.9) ∂ t̂ ∇ˆ · ψ∇ˆ ϕ̂ + |∇ψ|A ˆ ˆ 2 (ϕd − ϕ̂) = −A3 ψ(ĉ+ − ĉ− ), (3.10) where A1 = F U/RT , A2 = 1/λˆs , and A3 = F c0 l2 /(εU ) are three dimensionless parameters, and ∇ˆ = l∇. Note that we have assumed a closed system for these equations. Hereafter, we will drop the ‘∧’ notation in the equations for the readability. The SBM formulated NPP equations are the governing equations for simulating the electrodiffusion through com- plex geometries in which ψ defines the regions where evolutions of ionic concentration and electropotential field occur. While the dynamics is of importance in understanding the double layer formation, in some cases, only the steady-state ionic concentration and electropotential distributions are sought, for which the Poisson-Boltzmann model is used to obtain the equilibrium concentration and electropotential distributions (as described in Chapter 1). The SBM can also be used to solve the Poisson-Boltzmann equations. At equilibrium, the concentration of cation and anion no longer change with time. Thus, the mass transport equations (3.8) and (3.9) become ∇ · ψ(D+ ∇c+ + A1 D+ c+ ∇ϕ) = 0 =⇒ ∇ · ψ(∇c+ + A1 c+ ∇ϕ) = 0, (3.11) ∇ · ψ(D− ∇c− − A1 D− c− ∇ϕ) = 0 =⇒ ∇ · ψ(∇c− − A1 c− ∇ϕ) = 0, (3.12) in which the diffusivities in the equations on the left are assumed to be constant in a dilute electrolyte solution, and thus cancelled. The Poisson’s equation remains the same as Eq. (3.10). Equations (3.11), (3.12), and (3.10) are coupled differential equations, which are extremely unstable if solved by typical iterative methods (i.e., solving c+ and c− based on ∇ϕ, then solving ϕ based on the obtained c+ and c− , and iterating until numerical equilibrium). Thus, to solve the static equations, we introduce the Slotboom variables: η+ = exp(A1 ϕ + ln c+ ) = c+ exp(A1 ϕ) and η− = exp(−A1 ϕ + ln c− ) = c− exp(−A1 ϕ), as in Slotboom’s work [82]. Substituting them into Eqs. (3.11), (3.12), and (3.10) and rearranging terms, we obtain ∇ · (ψ exp(−A1 ϕ)∇η+ ) = 0 and ∇ · (ψ exp(A1 ϕ)∇η− ) = 0, (3.13) 53 ∇ · ψ∇ϕ + |∇ψ|A2 (ϕd − ϕ) = −A3 ψ[η+ exp(−A1 ϕ) − η− exp(A1 ϕ)]. (3.14) Furthermore, Slotboom [82] introduced a linearization procedure with ϕ = ϕ0 + δϕ and used Taylor expansion to approximate exp(±A1 δϕ) ≈ 1 ± A1 δϕ to obtain  ∇ · ψ∇δϕ − |∇ψ|A2 δϕ − ψA3 A1 η+ exp(−A1 ϕ0 ) + η− exp(A1 ϕ0 ) δϕ =  (3.15) − ∇ · ψ∇ϕ0 − |∇ψ|A2 (ϕd − ϕ0 ) − A3 ψ η+ exp(−A1 (ϕ0 ) − η− exp(A1 (ϕ0 ) . An iterative scheme is utilized where ϕ0 is the value from the former iteration and δϕ is a perturbation on ϕ0 . δϕ is solved to obtain a new ϕ = ϕ0 + δϕ, which is then substituted into Eq. (3.13) for solving η+ and η− . The process is repeated until η+ , η− , and ϕ all reach numerical equilibrium. This iterative method above was proved to be convergent by Varga [83] since the matrix is symmetric and positive definite if written in the matrix form. 3.2.2 Simulation Parameters The parameters used in the simulations are discussed here. Chemical diffusivity of ions can  be expressed as Dichem = Di 1 + ∂ ln γi /∂ ln ai , where Di is the diffusivity at the dilute limit, ai = γi ci is the chemical activity, and γi is an activity coefficient [84]. The second term in the parenthesis represents the thermodynamic factor. An assumption is made in this work that the temperature is constant at 298 K, and the electrolyte is a dilute solution, such that the thermodynamic factor is zero and thus ionic diffusivities are constant. The cation diffusivity for Li+ is D+ = 1.25 × 10−6 cm2 /s, and the anion diffusivity PF− 6 is D− = 4.00 × 10 −5 cm2 /s [85, 86, 87]. The dielectric constant ε in Poisson’s equation, Eq. (3.10), is set to be the value of dielectric constant of ethylene carbonate at 298 K [88]. Based on the prediction that the dielectric constant increases as salt concentration increases [89], we extrapolate the dielectric constant to a projected value at 1 M of LiPF6 dissolved in ethylene carbonate which is 8.85 ×10−8 F/m. Debye length is used to approximate the characteristic thickness of double layer (the Stern layer thickness λs ) [90]: s εRT λD = P 2 (3.16) F2 i zi ci0 54 where ci0 is the initial average concentration of species i in the electrolyte. Thus, we obtain λs = 0.5 nm. This value is kept as a constant in the following simulations by neglecting the variation of Stern layer thickness with loading potential and electrode surface properties. In this work, the width of the diffuse double layer is approximately 20 nm. Roughly 10 grid cells are set to span the diffuse double layer to ensure a sufficient resolution. In the 2D simulations, where the root level cell width is 256 nm, this requires 7 levels of refinement. It should be also noted that since the Euler explicit time scheme is used for solving Eqs. (3.8) and (3.9), the stable time step size must be smaller than the square of the length of the smallest cell size (∆t < min(∆x)2 ). As a result, the time step size used in the simulations needs to scale 1/4 as one more level of AMR is implemented. 55 3.3 Results and Discussion 3.3.1 One-dimensional (Pseudo-1D) Simulations To verify the correctness and accuracy of the presented method, a series of pseudo-1D sim- ulations were performed. (Simulations are performed on 2D domains but the results are effectively 1D. Hereafter, they are referred to as 1D simulations.) In these 1D simulations, the left and right boundaries represent two metal plates and the region between is filled with a uniform, stable electrolyte. This setup is the same as in Bazant’s work [90]. Figure 3.4(b) illustrates the domain with a 3-level refinement. It can be clearly seen that the grid becomes denser as being close to the boundary in the magnified view. The centers of dif- fuse interfaces (ψ = 0.5) in the SBM represent the left and right electrode surfaces and are located at x = 13 and 113 nm, respectively. Thus, the effective domain for electrodiffusion and electropotential calculations is 100 nm (with ∆x = 1 nm for the Root-Level grid). A constant voltage of 0.05 volt is imposed by setting the electropotential, ϕd in Eq. (3.10), to be 0.05 volt on the right boundary and 0 volt on the left boundary. Thus, the left and right boundaries serve as the cathode and anode surfaces, respectively. Snapshots of ion concentrations and electropotential profiles taken during the simulation with a 3-level refinement are provided in Fig. 3.5. Because the external potential is lowest at the cathode surface, the cations are quickly adsorbed and aggregated while the anions are repelled there, resulting in a peak of cation concentration and a depletion of anions on the cathode surface (see Fig. 3.5(a)). As more cations aggregate on the cathode surface, a small depletion of cations from the bulk value can be observed at the location slightly away from the cathode surface; see the red curve near x = 22 nm corresponding to t = 50 s in Fig. 3.5(b). t is a non-dimensional number used in 1D simulation indicating the output frame. The anion concentration profile corresponding to the same time is still smooth at the same location. This cation concentration depletion is an indication that the cation diffusivity is low such that it requires a longer time for cation supply from the bulk of electrolyte to the cathode 56 Figure 3.4 (a) Adaptive mesh refinement applied to electrode boundaries with 7 levels of quadtree refinement. The grid system is used in the 2D simulation. (b) Domain and refinement used in the 1-D simulations. There are three levels of quadtree refinement for the grid system. surface. The concentration profiles of cations and anions gradually reduce back to the bulk value in the electrolyte as the distance from the cathode surface increases. In the vicinity of the cathode surface, cation concentration is much higher than that of anion, resulting in a charge separation region with a thickness of approximately 8 nm (14 < x < 22 nm). This region is the diffuse double layer. On the anode surface, charge separation occurs in an opposite manner where anion concentration is much higher than cation concentration. Similar to the cation case, depletion of anion concentration immediately outside the diffuse double-layer region also occurs but in a much lesser degree due to the anions’ high diffusivity. In the early stage, electropotential drop across the electrode surface is relatively small and the gradient of electropotential in the bulk of electrolyte is large. See the red curve in Fig. 57 Figure 3.5 (a) 1D cation and anion concentration profile evolution obtained by SBM with 3-level quadtree AMR. (b) Magnified view of (a) from 15 nm to 35 nm (c) Evolution of electropotential profile corresponding to the same times in (a). (d) Magnified view of (c) from 13 nm to 30 nm. A discontinuity at x = 13 nm represents the electropotential difference across the Stern layer. 3.5(c). This electropotential drop represents the potential drop across the Stern layer. As the double layer is developed, the drop of electropotential across the cathode surface becomes large while the gradient of electropotential in the bulk of electrolyte decreases. Once reaching the steady state, the electropotential profile becomes completely flat in the bulk electrolyte region (see the green curve in Fig. 3.5(c)), indicating that the charge separation in the double layer (including Stern layer and diffuse layer) has completely counterbalanced the polarization induced by the imposed voltage. To investigate the effect of interfacial thickness in the SBM simulation of double layer formation, five sets of simulation results, using Level-0 through Level-3 refinement and a 58 sharp interface method, for the same domain are presented. The cation concentration profiles at the steady state are shown in Fig. 3.6(a) & (b). These results are obtained from solving Eqs. (3.8) through (3.10) until ∂c+ /∂t, ∂c− /∂t, and the variation of ϕ between time steps are all negligible. Level-0 is the original uniform grid system. As the refinement level increases, finer grids are used in the region of SBM diffuse interface. Since the width of the SBM diffuse interface is kept to be approximately 4 grid spacings, the width of SBM diffuse interface for Level-3 is only 1/8 of that in Level-0 case. The sharp-interface method is implemented by using a plain 1D finite difference method to solve the NPP governing equation, as in Ref. [42], with the domain boundaries at x = 13 and 113 nm. In the sharp- interface simulation, the grid spacing is set to be uniform and equal to the smallest grid spacing in Level-3 refinement, which is 1/8 of ∆x in Level-0 case. The analytical steady- state cation concentration distribution, calculated using the electropotential that is obtained by solving the Poisson-Boltzmann equation [29], is also provided in the same figure. The corresponding electropotential profiles in the electrolyte are provided in Fig. 3.6(c) & (d). Note that since the Poisson-Boltzmann model has assumed the electrolyte is an infinite reservoir, the ion depletion in the bulk of electrolyte (observed in Fig. 3.5(b)) due to charge separation does not appear. The shapes of all the simulated cation concentration profiles in Fig. 3.6(a) are similar, but they are slightly different near the boundaries. Simulations with more levels of refinement at the boundary regions give results closer to the analytical ones. The increase in accuracy can be clearly observed in the magnified view in Fig. 3.6(b). More accurate SBM results with higher levels of refinement are also observed in the simulated electropotential profiles; see Fig. 3.6(c) & (d). Because a higher level of refinement effectively reduces the diffuse interface thickness in the SBM, the curves for Level-3 refinement closely approach those of the sharp-interface result. Both the Level-3 and sharp-interface results almost overlap the analytical ones. The relative error between the steady-state results of the Level-3 AMR SBM and the sharp-interface is less than 0.05%. The errors between SBM and sharp-interface 59 Figure 3.6 (a) 1D cation distribution at the steady state obtained by the sharp interface simulation and diffuse interface SBM simulations with four different levels of refinement. (b) Magnified view of (a) from 14.6 nm to 16.1 nm. (c) Potential field corresponding to (a). (d) Magnified view of (c) from 14.75 nm to 17.40 nm. The x scales in (b) and (d) are selected to highlight the differences between the curves. methods during transient state are observed to be in the same order of magnitude. Thus, the comparisons here demonstrate that the SBM with a sufficiently thin diffuse interface can reach approximately the same result as that from conventional sharp-interface method. In a pure 1D domain (grid system), a Level-3 refined case will lead to a total of 166 grid points. On the other hand, a uniform grid system, with the grid spacing equal to the smallest grid spacing in the Level-3 refined case, will result in a total of 1040 grid points. Since the computational cost is approximately proportional to the number of grid points 60 in the simulations, the AMR reduces the computational effort by roughly 6 folds, but still maintains the numerical accuracy comparable to the sharp-interface results. 3.3.2 Two-dimensional Simulation with Complex Geometry With the AMR-SBM approach verified, we extended the simulation to 2D case. Figure 3.4(a) shows the geometries of two arbitrarily shaped particles used in the 2D simulation. The average diameters of the two particles are approximately 10 and 15 µm, and the center- to-center distance between the two particles is roughly 35 µm, which is almost 4375 times the diffuse layer thickness of the double layer. Quadtree refinement of 7 levels is taken to generate the grid system. The smallest grid spacing is 2 nm, and the largest grid spacing is 256 nm. The particles on the left and right serve as the cathode and anode particles, respectively, in a virtual electrochemical capacitor. The rest of the domain is occupied by electrolyte. Initially, the cation and anion concentrations are uniformly 1 M throughout the electrolyte region. Electropotential on the cathode and anode particles was set to be 0 V and 0.05 V, respectively, to induce a voltage across the two electrodes. No-flux boundary conditions are imposed on the four sides of the computational domain. Figures 3.7(a) through (d) shows the snapshots of anion distribution corresponding to physical times from t = 0.133 to t = 2.13 ms. Here, anion concentration is presented because of more prominent evolution during the simulation. The anion concentration rises surrounding the anode particle but drops around the cathode particle. This behavior is similar to that in the 1D cases. However, in two dimensions, the concentration can be observed to be nonuniformly distributed surrounding the particles. Two distinct time scales for the concentration evolution process are observed. First, ions are adsorbed on the particle surface very rapidly without long-distance diffusion. The ion concentration changes more rapidly in response to the imposed voltage when the local dis- tance between cathode and anode surfaces is shorter. For example, the anion concentration increases more rapidly on anode surfaces that are located closer to the cathode than on anode 61 Figure 3.7 Simulated anion concentration in the electrolyte taken in side-view (2D projection in the y-direction) at t = (a) 0.133 ms (b) 0.266 ms (c) 0.532 ms and (d) 2.13 ms.(e)Anion concentration taken in a different view angle to highlight the difference in the concentration at the convex and concave region. The image is taken at the same time as (a) surfaces that are located farther away from the cathode; see Fig. 3.7(a). Additionally, it can be observed that the anion concentration at the convex surface facing the cathode particle is higher than that at the nearby concave surface during the evolution because the convex regions are effectively closer to the cathode particle, as indicated in Fig. 3.7(e). Similarly, the anion concentration decreases more rapidly on cathode surfaces that are located closer to the anode than on cathode surfaces that are located farther away from the anode, as can be seen in Fig. 3.7(a). Later, the effect of long-distance diffusion plays a role. For exam- 62 ple, as shown in Fig. 3.7(b) and (c), the anions start to diffuse around the anode particle from the side facing cathode particle to the back of the anode particle, smoothing the anion concentration variation around the particle. As the system reaches the steady state, the anion concentration eventually reaches a uniform distribution over the cathode surface as well as the anode surface; see Fig. 3.7(d). Similar dynamics, but in an opposite manner, are observed for cation concentration evolution. The electropotential distributions corresponding to Fig. 3.7 are provided in Fig. 3.8. At the early stage, a significant gradient of the potential over the entire domain is observed. A large potential drop across the double layer is seen on the anode surface that is located closer to the cathode particle while the potential drop on the back side of the anode particle is small; see Fig. 3.8(a). The high electropotential drops on the particle surfaces near the regions with the narrowest inter-particle-surface distance causes the rapid ion absorption and repelling there during the first stage mentioned above. As the simulation progresses, the electropotential field within the electrolyte begins to flatten. The potential drop across the double layer becomes more uniform over the whole particle surface. See Fig. 3.8(b) and (c). At the steady state the electropotential is completely flat over the entire electrolyte domain and the potential drop becomes uniform on the whole particle surfaces, as shown in Fig. 3.8(d). It is important to point out that the simulation shows that this system requires at least 2 ms to reach its steady state. This suggests that, in an electrochemical impedance spec- troscopy (EIS) measurement for a system with a similar physical dimension, the steady state cannot be reached if the oscillation frequency of the load is higher than 500 Hz. Therefore, the measured double layer capacitance from an EIS study is far away from the equilibrium capacitance when frequencies on the order of kHz or MHz are used. Furthermore, as discussed in Section 3.2.1, the steady-state solution of cation concen- tration, anion concentration, and electropotential distribution can be directly solved with Slotboom variables. This is convenient when the transient behavior is not of interest and 63 Figure 3.8 Potential field in the electrolyte at (a) 0.133 ms (b) 0.266 ms (c) 0.532 ms and (d) 2.13 ms. Each subfigure corresponds to that in Fig. 3.7. only the steady-state solution is desired. Figure 3.9 shows the steady-state results from solving Eq. (3.13) and (3.15) for the presented 2D system. These results are the same as the steady-state results from the transient SBM simulation. This demonstrates the versatility of SBM formulation for solving NPP types of coupled equations. 64 Figure 3.9 Steady state concentration and potential distribution in the electrolyte obtained using Slotboom method: (a) cation concentration (b) anion concentration (c) potential distribution 3.3.3 Three-dimensional Simulation The presented method is also extended to 3D simulation. Two electrodes are placed within a 3D domain initially of 128 × 64 × 56 (x × y × z) root-level grid points; see Fig. 3.10(a). A planar plate serves as the cathode on one end, and a dendrite-like particle is used as the anode particle on the other side. The cathode-electrolyte interface is set at x = 96 nm. The shape of anode particle is inspired by phase field modeling of dendrite growth [91, 92]. The Root-Level grid spacing is equal to 8 nm, and the effective radius of the anode particle is roughly 160 nm. As a demonstration simulation, only 2 levels of refinement are performed. As a result, the smallest grid spacing is 2 nm. The potential on the cathode and anode is set to be 0 V and 0.005 V, respectively. Similar to the 2D simulation, no-flux boundary conditions are imposed on the 6 sides of the computational domain. The dynamic behavior of ionic concentration and electropotential evolution is similar to the 2D case. Therefore, we do not repetitively describe it. Figure 3.10 shows one result taken before the system reaches steady state. A pronounced diffuse double layer has formed on the cathode and anode surfaces; see the bright yellow region near the cathode surface and the dark blue region near the particle surface in Fig. 3.10(b). In the meantime, the electropotential field still shows a noticeable gradient throughout the electrolyte phase; see Fig. 3.10(c). This simulation clearly demonstrates the capability of the presented SBM-AMR 65 Figure 3.10 (a) Domain used in 3D simulation. The dendrite-like particle served as anode particle and the plate on the other side served as the cathode. (b) Simulation results for cation concentration. The bright yellow region in the cathode and the dark blue region on the particle indicate double layer regions. (c) Simulation results for potential distribution in electrolyte. approach in simulating double layer formation in 3D complex geometries. 66 3.4 Conclusion In this chapter, the SBM with AMR approach is applied for simulating the double layer formation. It can resolve the thin double layer while handling the complex geometry of electrode particles. One-dimensional simulations are performed to verify the correctness and accuracy of this method. Two-dimensional simulations with irregular shaped electrodes are performed to demonstrate the importance of considering the geometric effect that cannot be included in 1D simulations. Both dynamic process during transient states and the steady- state results obtained from the direct solver are studied. A clear two-stage ion concentration evolution has been observed during the double-layer formation. Three-dimensional simula- tion is also provided to highlight the capability of presented method. With this method, further investigation into the effects of material properties, such as diffusivity, dielectric con- stant, as well as more complicated geometries, on electrochemical double layer formation can be performed. 67 CHAPTER 4 PHYSICAL-BASED SIMULATION OF ELECTROCHEMICAL IMPEDANCE SPECTRUM ON CATHODE MATERIAL 68 4.1 Introduction A battery’s macroscopic performance is determined by many microscopic processes simul- taneously occurring in the electrodes. For example, deintercalation of the cathode involves Li outward diffusion in the particles, electron migration through the particle percolation network, electrochemical reaction on the particle surfaces, and Li-ion diffusion in the elec- trolyte, all of which take place in a complex electrode microstructure. As can be envisioned, the microstructural features, such as porosity [1], tortuosity [2] of interparticle space, re- active surface area [3], and particle size [1, 4], all affect the battery’s performance. These microstructural features sometimes vary through the electrode, which can result in substan- tially non-uniform reaction current density and might deteriorate cell performance or even cause safety concerns [5]. Without a comprehensive understanding of the relationship be- tween electrode microstructure and battery performance, it is impossible to improve electrode designs or optimize battery performance. Electrochemical impedance spectroscopy (EIS) is a widely used technique to measure properties of electrochemical devices, such as batteries and fuel cells. The EIS measures the response current (or voltage) to an oscillating voltage (or current) loading [6]. The device’s resistance, capacitance, and Warburg impedance are evaluated by fitting the obtained EIS curve to an equivalent circuit model (ECM) comprised of resistors, capacitors, and constant- phase elements [7]. These resistance, capacitance, and Warburg impedance correspond to the charge transfer reaction, polarization, and diffusion (mass transport) involved in the whole electrochemical process. In this chapter, we employed SBM and AMR as described in Chapter 2 to handle the complex electrode microstructures so that the implementation of simulations can be greatly accelerated. This work emphasizes the ease of using non-conformal mesh in simulating coupled electrochemical phenomena in complex microstructures. The SBM-reformulated governing equations can be solved using other numerical methods (e.g., FEM or FVM), not limited to those presented (FDM) in this work. Further improvements using more advanced 69 numerical methods can greatly accelerate the complex microstructure simulations and make them accessible tools for the research community. In order to extract EIS curves from physics-based microstructure electrochemical simula- tions, capacitance due to double-layer formation was calculated based on the Nernst-Planck- Poisson model separately [73], which is presented in Chapter 3. The obtained capacitance is incorporated into the electrochemical simulations in calculating the total current respond- ing to the oscillating cell voltage loadings. Similar to Refs. [68, 69, 93], synthetic NMC microstructures were generated using discrete element method. We examine the effects of initial Li fractions (state of charge) in the cathode, average salt concentration in the elec- trolyte, and microstructures on the EIS curves. The obtained charge-transfer resistance is inversely proportional to exchange current density, which is a function of state of charge of the cathode particles. While the salt concentration in the electrolyte can simultaneously affects double layer capacitance, ionic diffusion, and exchange current density, the simulation results indicate that the change in the exchange current density dominates the variation of semicircle diameter of the EIS curve. With different cathode microstructures, the change in active surface area determines the variation of total charge-transfer resistance of the elec- trode. The simulated EIS curves on the Nyquist plot well agree with physical interpretation and demonstrate the capability of physics-based microstructure simulations for EIS behavior. 70 4.2 Model Construction and Simulation Details 4.2.1 Electrochemical Model During an EIS measurement, a current arises responding to an oscillating voltage loading. Generally, the amplitude of loading is sufficiently small such that the amplitude of response (concentration and electrostatic potential variations) are also small. Thus, the variations ap- proximately linearly respond to the loading. While this macroscopic behavior appears intu- itive, there are multiple coupled electrochemical processes occurring in the complex electrode microstructures. In a full-cell battery, those processes take place in five regions: cathode particles, cathode-electrolyte interface, electrolyte, electrolyte-anode interface, and anode particles. Even in a simplified half-cell case, which consists of only a cathode, electrolyte, and metallic Li anode, multiple electrochemical processes still need be considered, which in- clude (1) Li-ion transport and electric-current flow in the solid cathode particles, (2) solute ion transport and ionic-current flow in the liquid electrolyte, and (3) the electrochemical intercalation reaction at the cathode-electrolyte interface. Those processes are described by the classical electrochemical governing equations [79] within their respective domains. They are briefly presented below. More detailed descriptions can be found in Ref. [93]. In cathode particles, Li diffuses through interstices in the host crystal, as described by ∂Cp  ∂X  = ∇ · Dp ∇Cp ∈ Ωp =⇒ = ∇ · Dp ∇X ∈ Ωp (4.1) ∂t ∂t where Cp , X, and Dp are the Li concentration, Li site occupancy fraction, and diffusivity in the particles, respectively. Cp = ρX, where ρ is Li site density of the cathode crystal. t and Ωp denote time and the domain of the particle. Here, a simple Fickian diffusion is assumed for the Li transport. For more complicated Li transport behavior, the phase field method can be utilized to describe the concentration evolution [94, 95, 96, 97]. Li insertion/extraction occurs on particle surfaces via electrochemical intercalation reaction, which is described as a Neumann boundary condition (specifying flux or gradient): rxn /ρ = np · jp ∈ ∂Ωp for Eq. (4.1), where rxn is the surface reaction rate, np is the inward surface normal vector, and jp is 71 the intercalation flux. The electron current continuity in the particle regions is described as  ∇ · κp ∇ϕp = 0 ∈ Ωp , (4.2) where κp and ϕp are the electrical conductivity and electro-potential in the particles. This  equation is subject to the boundary condition, np · ip = np · − κp ∇ϕp = z+ F rxn ∈ ∂Ωp , where ip is the electrical current, z+ is the charge number, and F is the Faraday constant. Here, we have ignored the additive phases for the clarity of presenting equations, which can be included by introducing another domain that allows electron conduction but no Li diffusion. Here, for simplicity, the electrolyte is assumed to be a monovalent binary electrolyte, which contains only one species for cation and one species for anion with +1 and −1 charge number, respectively. The diffusion and migration processes of ions can be described as ∂Ce  ie · ∇t+ = ∇ · De ∇Ce − ∈ Ωe , (4.3) ∂t z+ ν+ F where Ce , De , and ie are the salt concentration, the ambipolar diffusivity of the salt, and the ionic current in the electrolyte, respectively. Ωe indicates the domain of electrolyte. νi and ti are the dissolution number and transference number, respectively, where the subscript + denotes cation. The salt concentration is related to ion concentrations by Ce = ν+ C+ = ν− C− . The ambipolar diffusivity is related to the ionic mobilities and diffusivities by De = (z+ m+ D− − z− m− D+ )/(z+ m+ − z− m− ) where mi is the transport mobility and Di is the diffusivity of the ions. The transference number of cation is t+ = z+ m+ /(z+ m+ − z− m− ) = 1−t− . In this work, z+ = 1 and ν+ = 1. Equation (4.3) is subject to the boundary condition: rxn = ν+ (ne · je ) ∈ ∂Ωe , where je = −De ∇Ce + t+ ie /(z+ ν+ F ) is a virtual salt flux. Note that ne is the inward normal to the electrolyte and ne = −np . If t+ is constant, the second term in Eq. (4.3) vanishes. In the simulation work presented later, the Einstein relationship (mi = Di /RT ) is assumed for simplicity in calculating the values of ambipolar diffusivity. The double layer thickness (a few tens of nm) is negligible in the microstructure scale (µm). Thus, current continuity is still assumed in the electrolyte at the microstructure-scale 72 equations, i.e., ∇·ie = 0 ∈ Ωe , where ie = −F z+ ν+ [(z+ m+ −z− m− )F Ce ∇ϕe +(D+ −D− )∇Ce ]. This leads to the governing equation of electrostatic potential in the electrolyte as ∇ · [(z+ m+ − z− m− ) F Ce ∇ϕe + (D+ − D− ) ∇Ce ] = 0 ∈ Ωe . (4.4) The coefficients in the first term of Eq. (4.4) are related to the electrolyte’s electric con- ductivity by κe = (z+ m+ − z− m− )F 2 Ce . The boundary condition for Eq. (4.4) is ne · ie = z+ F rxn ∈ ∂Ωe . The intercalation rate through particle-electrolyte interfaces is determined by the Li concentrations and electrostatic potentials on the different sides of the interfaces and is expressed as the Butler-Volmer equation [79]:     −αz+ F p (1 − α) z+ F p rxn = kf Ce exp [ϕ]e − kb Cp exp [ϕ]e (4.5) RT RT where kf and kb are the forward and backward rate constants, respectively, α is the symmetry factor, R is the ideal gas constant, T is the absolute temperature, and [ϕ]pe = ϕp − ϕe ∈ ∂Ωp is the electro-potential drop across the particle-electrode interface. This equation provides the necessary boundary condition for solving Eqs. (4.1) through (4.4). 4.2.2 Smoothed Boundary Method Solving partial differential equations in the conventional sharp-interface methods requires discretization of the domain with a mesh system conformal to that domain (e.g., using FEM). The mesh generation processes for complex geometries are very time-consuming, which hinders complex microstructure simulations. In this chapter, we continue employ- ing the SBM [70, 73] to circumvent the tedious mesh generation process that are required in the sharp-interface methods. A brief derivation of the SBM-formulated electrochemical governing equations is provided below. Here, ψp is used to define the regions of electrode particles: ψp = 1 inside the particles and ψp = 0 outside. Multiplying ψp on both sides of Eq. (4.1), we obtain ∂Xp ψp = ψp ∇ · (Dp ∇Xp ). (4.6) ∂t 73 Using the product rule of differentiation on the right-hand side, this equation can be further written as ∂Xp ψp = ∇ · (ψp Dp ∇Xp ) − ∇ψp · (Dp ∇Xp ). (4.7) ∂t The second term on the right-hand side serves as an ‘internal’ boundary condition embedded inside the computational domain. In the diffuse-interface description, np = ∇ψp /|∇ψp |. The Neumann boundary condition on the particle surface (rxn /ρ = np ·jp ∈ ∂Ωp ) can be expressed as rxn ∇ψp = np · jp = · (−Dp ∇Xp ). (4.8) ρ |∇ψp | Thus, we obtain ∂Xp 1  |∇ψp | rxn = ∇ · ψp Dp ∇Xp + (4.9) ∂t ψp ψp ρ as the SBM equation for the Li fraction evolution in the particles. Following a similar procedure, the SBM formulation for Eq. (4.2) can be obtained as ∇ · (ψp κp ∇ϕp ) − |∇ψp |z+ F rxn = 0. (4.10) Similar to the derivation of Eq. (4.9), we multiply Eq. (4.3) with ψe , which is the domain parameter of electrolyte and ψe = 1 − ψp . The obtained equation is written to ∂Ce 1 1 ie · ∇t+ = ∇ · (ψe De ∇Ce ) − ∇ψe · (De ∇Ce ) − . (4.11) ∂t ψe ψe z+ ν+ F Recall that De ∇Ce = −je + t+ ie /(z+ ν+ F ), rxn /ν+ = ne · je , z+ F rxn = ne · ie , and ne = ∇ψe /|∇ψe |, Eq. (4.11) is reorganized to ∂Ce 1 |∇ψe | rxn t− ie · ∇t+ = ∇ · (ψe De ∇Ce ) + − . (4.12) ∂t ψe ψe ν+ z+ ν+ F In the cases where the transference numbers are constant, the last term vanishes. Lastly,   using the relation z+ F rxn = ne · ie = (∇ψe /|∇ψe |) · − z+ ν+ F (z+ m+ − z− m− )F Ce ∇ϕe +  (D+ − D− )∇Ce , the SBM equation for current continuity in the electrolyte is obtained as rxn ∇ · [ψe (z+ m+ − z− m− ) F Ce ∇ϕe ] + |∇ψe | = ∇ · [ψe (D− − D+ ) ∇Ce ] . (4.13) ν+ 74 In summary, Eqs. (4.9), (4.10), (4.12), and (4.13) are the SBM-formulated equations based on the classical Eqs. (4.1), (4.2), (4.3), and (4.4), respectively. These equations can be solved on grid systems (mesh) nonconformal to the complex electrode microstructures, while imposing the reaction flux (Butler-Volmer equation) at the diffuse interfaces. At each time step, Xp and Ce are updated based on rxn by Eqs. (4.9) and (4.12), respectively. Since De is typically more than five orders of magnitude larger than Dp , the stable time step size for Eq. (4.9) is too large for Eq. (4.12). Thus, a fully implicit time scheme is used for Eq. (4.12) in order to use the same ∆t of Eq. (4.9). Within each time step, the static Eq. (4.10) is solved for ϕp with rxn from Eq. (4.5) as a flux boundary condition on the particle-electrolyte interface and ϕp |c as a Dirichlet boundary condition (specifying value) on the cathode current collector (computational domain boundary). Similarly, static Eq. (4.13) is solved for ϕe with rxn as a flux boundary condition on the particle-electrolyte interface and ϕe |a as a value boundary condition on the metallic anode surface (on the opposite side of computational domain boundary). The obtained ϕp and ϕe are substituted back to Eq. (4.5) to calculate a new rxn . This process is repeated until ϕp , ϕe , and rxn all reach numerical equilibrium. Then, the next time step starts. In this work, the implicit time evolution Eq. (4.12), static Eqs. (4.10) and (4.13) are solved using a standard Jacobi relaxation method. More aggressive relaxation methods can be employed to accelerate simulations if needed, which however is beyond the scope of the current work. 4.2.3 Double Layer Capacitance Electrochemical double layers form near the charged electrode surfaces to balance the ex- ternally imposed electrostatic potential field. Charge separation between cations and anions in the diffuse double layer generates the capacitance observed in EIS measurements. The ion concentration and electrostatic potential evolution during double layer formation can be described by the Nernst-Planck-Poisson (NPP) equations as shown in Chapter 3. For 75 readers’ convenience, they are provided again here:   ∂C+ = −∇ · J+ = ∇ · D+ ∇C+ + z+ m+ F C+ ∇ϕe , (3.1) ∂t   ∂C− = −∇ · J− = ∇ · D− ∇C− + z− m− F C− ∇ϕe , (3.2) ∂t F ∇2 ϕe = − (C+ − C− ), (3.3) ε where Ji = −(Di ∇Ci + zi mi F Ci ∇ϕe ) is the flux vector and ε is the dielectric constant (or relative permittivity) of the electrolyte. The SBM has been utilized to solve the NPP equa- tions for complex geometries previously. The relevant SBM equations and numerical details can be found in Chapter 3 and Ref. [73]. Based on the simulations, the charge separation reaches the steady state to form double layers in about 2 ms for the electrolyte considered in this work and the double layer thickness is approximately 20 nm. Since the thickness of double layer is several orders of magnitude smaller than other characteristic scales, it is difficult to solve the NPP equations for explicit ion concentrations and electrostatic potential within the double layer regions while obtaining the concentrations and potential distribu- tions spanning throughout the electrode scales. Although it is possible to accommodate both double layer formation and microstructure electrochemical processes together in one SBM simulation, extremely high levels of mesh refinement is required to resolve the huge spatial discrepancy, which would demand an enormous computational resource. Thus, we did not pursue such routes. There are two methods to calculate the resulting capacitance from the charge separation in the double layers [98]. In the time domain, the capacitance can be computed from the ratio between the total separated charges within the double layer and the potential difference across the double layer: Cdl = ∆q/∆ϕ. This type of calculations are performed using the simulation results from the steady states. In the frequency domain, the capacitance can be extracted from the impedance as Cdl = −1/(f ZIm ), where ZIm is the imaginary part of the total impedance and f is the ordinary frequency in the unit of Hz. The angular frequency is ω = 2πf in the unit of radian per second. In the relevant simulations, an oscillating 76 potential (ϕe = Φ0 sin(ωt) ∈ ∂Ωe ) are imposed on the electrodes as the boundary condition, which is equivalent to the expression of Φ0 exp(jωt) ∈ ∂Ωe , where j is the imaginary unit. The response current is calculated by dq/dt = i = I0 exp(j(ωt + η)) ∈ ∂Ωe , where η is the phase shift. The impedance is then obtained by Z = Φ0 exp(jωt)/{I0 exp[j(ωt + η)]}. As discussed later, our simulations show that the double layer capacitance becomes constant when the frequencies are below ∼500 Hz. The difference between the results from those two methods mentioned above are negligible. Since the double layer formation (∼20 nm near the electrode particle surface) reaches the steady state within milliseconds, it is reasonable to assume that the double layer capacitance is a constant value in the EIS measurements that the loading frequencies are under 100 Hz. Therefore, in this chapter, the SBM with AMR techniques were used only on the electrochemical governing equations (Eqs. (4.9), (4.10), (4.12), and (4.13)) for simulating the diffusion and surface reaction phenomena. The double layer capacitance was calculated separately as a constant. 77 4.3 Results and Discussion 4.3.1 Specific Capacitance The NPP equations in Section 4.2.3 were solved to simulate the double-layer evolution under sinusoidal voltage boundary conditions in pseudo-1D cases. The oscillating voltage loading was imposed by setting ϕe = 0 on the left electrode surface and ϕe = V0 sin(ωt) on the right electrode surface, where V0 = 50 mV. The values of D+ and D− in Eqs. (3.1) and (3.2) were set to be 1.25 × 10−6 and 4.0 × 10−6 cm2 /s [99, 73], respectively. The dielectric constant in Eq. (3.3) was set to be 8.85 × 10−10 F/cm [73]. The details of SBM simulations of the NPP model can be found in Ref. [73]. Figures 4.1(a) and (b) show snapshots of C+ and ϕe profiles between the two electrode plates in the pseudo-1D simulation taken at t = 2.5 × 10−5 and 3.0×10−5 s at a loading frequency f = 10000 Hz. The two electrode plates were separated apart in a distance of 1.5 µm. At t = 2.5 × 10−5 s, cell voltage was highest, the C+ rose on the left plate and dropped on the right plate to balance the ϕe variations near the plate surfaces; see the red dashed line and green solid line in Fig. 4.1(a) and (b), respectively. At t = 3.0 × 10−5 s, the boundary value of ϕe on the right decreased to 45 mV and ϕe changed accordingly; see the green dashed line in Fig. 4.1(b). Responding to the voltage loading, C+ on the left decreased and that on the right increased, see the blue line in Fig. 4.1(a). C− profiles have similar shape to C+ but with opposite signs. For clarity of view, C− profiles are not presented in the figure. The amount of charge separation in the double layer can be R calculated by F (C+ − C− )dx. Here, the response current was calculated according to  I  dq d  I= = F C+ − C− ds , (4.14) dt dt where the circular integration is over one half of the electrolyte domain. This is because only the charge separation on one of the two electrode surface is enough for calculating the resulting current. The positive and negative charge separations will neutralize each other if intergated over the whole domain. The time derivative is taken using the first-order forward difference. The response current was also a sinusoidal function in time with a phase shift to 78 the voltage loading. See Fig. 4.1(c), in which the the red curve is the response current to the voltage loading (the blue curve). We used MATLAB curve-fit function to extract the phase shift and amplitude of the simulated current for constructing the Nyquist plot. Note that Figure 4.1 Pseudo 1D simulation of double layer evolution under oscillating loads: (a) the cation profile at t = 2.5 × 10−5 and 3.0 × 10−5 s, and (b) the electrostatic potential profiles corresponding to (a). (c) The voltage loading and current response in the simulation of f = 10000 Hz. (d) The Nyquist plot for electric double layer capacitor with different inter-plate distances extracted from the simulations. The impedance curve for (a)-(c) is not shown in (d). since the sinusoidal form of the response current were not completely developed in the first wave, the curve-fitting started from the third period. Figure 4.1(d) shows impedance curves extracted from pseudo-1D simulations with loading frequencies from 10 to 10k Hz. The 79 four curves are for inter-plate distances of 100, 200, 500, and 1000 nm. At low frequencies, the impedance curves are nearly vertical lines, indicating a constant capacitance. See the black curve for an example. Here, when f < 500 Hz, the evolution of the double layer is no longer limited by diffusion. The capacitance current (change rate of the amount of separated charges) exhibits a constant 90-degree phase delay to the voltage oscillation, thus leading to a constant capacitance. In contrast, a Warburg impedance behavior is observed in the high frequency regime: the curve bends toward the lower-left direction; see the black curve below f = 500 Hz. This can be attributed to the fact that the ionic diffusion is limiting the double layer evolution. Such an effect is more pronounced when the two plates are separated further. Based on the simulation results presented in Fig. 4.1(d), ZIm is nearly independent of the distance between the two electrodes, which can be attributed to the fact that the steady- state ion concentrations in the double layer is independent of the inter-plate distance due to the assumption in Eq. (3.3) that ϕe profile can be established regardless how far the two plates are separated. The capacitance is almost constant when the frequency is lower than ∼500 Hz because most of ionic diffusion occurs only within the double layer. As in Section 4.2.3, the impedance was calculated according to Z = V /I, from which  the capacitance in the frequency domain was obtained by Cdl = −1/ f ZIm . Note here Cdl is the total double-layer capacitance on one electrode surface. The specific capacitance is the capacitance per electrode surface area: Csp = Cdl /A. The values of specific capacitance obtained from the pseudo-1D simulations is 140 µF/cm2 , which is approximately four times an experimentally measured value, 29.4 µF/cm2 [100]. This deviation may be due to the extrapolated, high dielectric constant used in the NPP simulations [73]. Calculations with concentration-dependent ion diffusivities in Eqs. (3.1) and (3.2) were also performed to ex- tract the specific capacitance. The concentration-dependent ionic diffusivities are shown in the next section later. The obtained Csp was of negligible difference to that from constant- diffusivity simulations since the voltage loading amplitude (50 mV) was small. Thus, it is 80 valid to use the Csp values obtained from constant salt diffusivity simulations in further EIS simulations. Here, we also neglect possible variation of Csp due to change in loading potential. Moreover, specific capacitance was also calculated using two circular electrodes in 2D configurations. See Fig. 4.2 for simulated cation, anion, and electrostatic potential distributions taken at t = 1.505 × 10−4 s (under a loading frequency of f = 10000 Hz), and the Nyquist plot constructed based on the calculated impedance at several frequencies. The Figure 4.2 Simulation of two-circular-electrode case: (a) the cation, (b) the anion concentration, (c) the electrostatic potential distributions over the electrolyte and electrodes, and (d) Nyquist plot generated from impedance calculations under different frequencies impedance curve in Fig. 4.2(d) at the lowest simulated frequency bends slightly to the left due to the numerical error in curve fitting. The obtained Csp (∼ 136 µF/cm2 ) are close to those from the pseudo-1D simulations and is also nearly independent of frequencies when the frequency is below 500 Hz. Therefore, we will use the value of Csp from the two-plate simula- tions for the electrochemical simulations presented in the following sections. The capacitance 81 current is calculated according to ∂[ϕ]pe Z   Ic = Csp dA, (4.15) ∂t where [ϕ]pe is the voltage across particle-electrolyte interface, and dA indicates all active electrode particle surfaces. In the diffuse interface SBM method, we use ∂[ϕ]pe Z   Ic = Csp |∇ψ| dΩ (4.16) ∂t instead for the calculation, since the interface has a nonzero thickness. The time derivative is calculated explicitly between two successive time steps in the simulations. 4.3.2 Electrochemical Simulation In this chapter, we use NMC-111 (Lix Ni1/3 Mn1/3 Co1/3 O2 ) as the model cathode material, for which the secondary particles have a nearly spherical shape. It is convenient to use the discrete element method (DEM) [101, 102] to computationally generate synthetic mi- crostructures. 166 spheres following a truncated log-normal distribution of radius (with a mean radius µ = 5.8 µm, and variance σ 2 = 0.25) were created and initially randomly placed in the 3D domain. The particle radius distribution is shown in Fig. 4.3(a). Those spheres were relaxed with a downward body force such that they eventually descended to an agglom- erate at the bottom when equilibrium was reached; see the green spheres in Fig. 4.3(b). The region around the particle agglomerate was selected as the computational domain (height to 100 and 70×90 in the horizontal plane). The domain was truncated in the dimension of the horizontal plane such that the particle domain in the simulations will be similar to those cropped from a large electrode. This computational domain was discretized into a uniform grid system. The distances from each grid point to the particle surfaces were calculated using a level-set distancing method [103], and the domain parameter function ψ was obtained by   substituting the distance function to a hyperbolic tangent as ψ = 1 + tanh(d/ζ) /2, such that ψ = 1 in the particles and ψ = 0 outside, where d is the shortest distance to the particle surfaces and ζ is used to control the interfacial thickness. 82 Figure 4.3 (a) particle radius generated based on log normal distribution (b) NMC cathode accumulated by small spherical particles (c) half cell model with grey particles being NMC cathode, brown plate being the current collector and green plate being lithium metal (d) side view of AMR with 2 layers of refinement The virtual cell used in the simulations is shown in Fig. 4.3(c), in which the NMC agglomerate was truncated to fit the rectangular-prism computational domain as mentioned earlier. The root-level grid system contains 70 × 100 × 90 point in the x, y, and z directions, respectively, with the grid spacing ∆x = 1 µm. Note that this virtual cell has been rotated 90◦ in the z-x direction from that in Fig. 4.3(b). Two levels of octree refinement were implemented, which resulted in a total of 10,551,898 grid points. The value of ζ was set to 83 be ∆x/(22 ), such that approximately four Level-2 grids spanned over the diffuse interface. The time step size was ∆t = 1 × 10−4 s. Figure 4.3(d) shows the AMR mesh on the domain boundary at y = 0. The semitransparent cyan slab in Fig. 4.3(c) represents a metal Li anode and the dark brown slab represents the current collector on the cathode side. The transparent region is the space between particles and is filled by electrolyte. The empty space in the region (0 < x < 35 µm) serves as the separator. Based on our experiences [93], two levels of refinement are sufficient for SBM electrochemical simulations to reach results close to those obtained from traditional sharp interface simulations. The value of ϕe |a was set to be zero on the metallic anode as the boundary condition for solving ϕe and the ϕp |c = 0.03 sin(ωt) + ϕOCV (V) was imposed on the current collector to create an oscillating voltage loading, where ϕOCV is the open circuit voltage (OCV) at the respective state of charge. Li diffusivity in Eq. (4.6) and electric conductivity in Eq. (4.10) are functions of Li fraction in NMC crystals. The Li diffusivity in NMC particles used in this work is DLi = (0.0277 − 0.0840X + 0.1003X 2 ) × 10−8 cm2 /s (4.17) where X is the Li fraction. This function is fitted from experimental data taken from Ref. [104], in which Li diffusivity was measured from NMC disk pellets using EIS techniques. The electric conductivity of NMC as a function of Li fraction is given as κ = 0.0193 + 0.7045 tanh(2.399X) − 0.7238 tanh(2.412X) S/cm. (4.18) It is also fitted from experimental data in Ref. [104]. Since no data point is available in the low Li fraction region (X < 0.2), the curve was extrapolated to the low Li fraction region. At equilibrium, the net reaction described by Eq. (4.5) is zero, from which the reaction rate constants can be calculated according to     i0 αz+ F i0 (α − 1)z+ F kf = exp ϕeq and kb = exp ϕeq , (4.19) z+ F C+ RT z+ F Cp RT where i0 is the exchange current density and ϕeq is the equilibrium voltage drop across the particle-electrolyte interface. Here, ϕeq ≈ ϕOCV since the open circuit voltage is measured 84 Figure 4.4 (a) Li diffusivity in NMC, (b) electric conductivity in NMC, (c) diffusivity in the electrolyte, (d) OCV as a function of XLi , (e) exchange current density as a function of XLi at Ce =1 M, and (f) calculated forward and backward reaction constants from the exchange current density and OCV. The image is taken from Ref. [93]. at a near equilibrium condition. Both i0 and ϕOCV are functions of Li fraction in the NMC particles. Here, we assume α = 0.5. The open circuit voltage function of NMC cathode materials is ϕOCV = 1.095X 2 − 8.234 × 10−7 exp (14.32X) + 4.692 exp (−0.5389X) V. (4.20) This OCV function is fitted from the data taken from Ref. [105] using the formula in Ref. [106]. Since open circuit voltages are measured at quasi-equilibrium conditions in which the 85 Li concentration is nearly uniform throughout the whole cathode and the salt concentration is also nearly uniform at 1 M, the OCV can be viewed as the equilibrium potential. The exchange current density function is i0 = 10−0.2(X−0.37)−0.9376 tanh(8.961X−3.195)−1.559 mA/cm2 , (4.21) which is fitted from the data measured using EIS techniques on single NMC particles in Ref. [107]. Binary electrolyte was assumed in the simulations, in which the ionic diffusivities for Li+ and PF− 6 at 1 M of LiPF6 salt concentration in the electrolyte were 1.25 × 10 −6 cm2 /s and 4.0 × 10−6 cm2 /s, respectively, as in Ref. [99]. The ambipolar diffusivity of salt in the electrolyte is De = 0.004649 × exp(−7.02 − 0.83Ce + 0.05Ce2 ) cm2 /s. (4.22) The curve of this function has a shape following the measured data in Ref. [108], but the magnitude was scaled to match the values of D+ and D− (1.25 × 10−6 cm2 /s and 4.0 × 10−6 cm2 /s, respectively) at 1 M. In our test simulations, the difference between the obtained EIS curves from constant and variable ambipolar diffusivities was negligible because the amplitude of salt concentration variation was very small under the oscillating loading. Thus, ambipolar and ion diffusivities were later set to be constant values corresponding to the average salt concentration of 1 M in further simulations. As experimentally observed in Ref. [108], transference numbers were set to be constant. Thus, the last term in Eq. (4.12) vanished. The electrolyte conductivity can be calculated according to   D+ D− κe = z+ − z− F 2 Ce S/cm. (4.23) RT RT The details of parameterizing material properties can be found in Ref. [93]. The initial Xp was uniform 0.25 throughout all particles and the initial Ce was 1 M throughout the electrolyte. (Here, we use Xp and X as interchangeable notations for Li fraction in the electrode particles.) At Xp = 0.25, ϕOCV = 4.17 V. Figure 4.5 shows snapshots 86 of simulated Xp and ϕp , and Fig. 4.6 shows snapshots of simulated Ce and ϕe distributions under f = 0.1 Hz, taken at t = 1.51, 1.75, and 2.01 s. At lithiation phase, increased Xp can be observed on the particle surfaces in Fig. 4.5(a) because the Li diffusion time scale was comparable to the oscillation time scale. The Xp in the bulk of particles remains fairly uniform. As lithium was extracted from the particle surfaces, Li fraction on the particle surfaces decreases, as shown by Fig. 4.5(a), (b), and (c). Figure 4.5 Top row: Simulated Li fraction in the particles. Bottom row: Simulated electrostatic potential in the particles. The columns from left to right are taken at t = 1.51, 1.75, and 2.01 s, respectively, under a loading frequency of f = 0.1 Hz. The ϕp corresponding to Xp in Fig. 4.5(a) through (c) are shown in Fig. 4.5(d) through (f), respectively. The negative gradient of ϕp along the +x direction in Fig. 4.5(d) indicates an intercalation flux to the cathode, during which high Xp was seen on particle surfaces in in Fig. 4.5(a). At t = 1.75 s, the positive ϕp gradient along the +x direction in Fig. 4.5(e) indicates that Li was extracted from the particles. At this moment, Li near particle surfaces was extracted such that Xp in Fig. 4.5(b) was relatively uniform. Li extraction continued as 87 indicated by the positive ϕp gradient along the +x direction in Fig. 4.5(f), which led further Xp decrease on NMC particles, as observed in Fig. 4.5(c). Figures 4.6(a) through (c) are the simulated Ce taken at the times corresponding to Fig. 4.5(a) through (c), respectively. The negative Ce gradient near the x = 0 plane (anode- Figure 4.6 Top row: Simulated salt concentration in the electrolyte. Bottom row: Simulated electrostatic potential in the electrolyte. The columns from left to right are taken at t = 1.51, 1.75, and 2.01 s, respectively, corresponding to the three columns in Fig. 4.5. electrolyte interface) in Fig. 4.6(a) indicates that Li enters the computation domain at there to compensate those Li ions intercalated into NMC particles. At t = 1.75 s, as shown in Fig. 4.6(b), the gradient of Ce near the x = 0 plane was positive, which meant that Li ions were removed from that plane to maintain a conservation of ions during Li deintercalation. The Ce gradient remained positive in Fig. 4.6(c). Figures 4.6(d), (e), and (f) show the snapshots of simulated electrostatic potential field in the electrolyte, corresponding to Fig. 4.6(a), (b), and (c), respectively. Since the ϕe gradients are very small, individual color scales were used for each subplot to clearly show 88 the contrast of the gradients. The negative gradient near x = 0 in Fig. 4.6(d) indicates that Li ions enters the system through the anode-electrolyte interface, and the positive gradients near x = 0 in (e) and (f) indicate Li ion leave the system. The behavior of ϕe distributions are generally similar to that of Ce since their gradients are mainly determined by the Li ion flux in the system. However, some difference can still be discerned due to the different governing physics. For example, the undulation of Ce near x = 15 µm in Fig. 4.6(b) does not appear in (e). These results clearly demonstrate the versatility and capability of using SBM with AMR to properly simulate the coupled multiphysics electrochemical processes in the electrode microstructures. During the simulation, the total response current was recorded as shown in Fig. 4.7(a). The total current includes the reaction and capacitance currents and is calculated according to Z Itot = F |∇ψp |rxn dΩ + Ic , (4.24) where Ic is given in Eq. (4.16). MATLAB curve-fit function was used to extract the phase shift and amplitude of the response current, which allow the calculation of impedance using the method mentioned in Section 4.2.3. Similar simulations were performed for 20 other loading frequencies, ranging from 0.001 to 200 Hz. Figure 4.7(b) shows the EIS curve on a Nyquist plot constructed from the 21 sets of simulations. The obtained EIS curve has a well-developed semicircle in the intermediate-to-high fre- quency range: 0.25–200 Hz. In the results, 0.25 Hz is considered to be the minimal frequency as it is on the right-hand side end of the semicircle. The critical frequency is the point at the middle of the semicircle (the maximum height on the semicircle), which is approximately 5 to 6 Hz on the simulated EIS curve. On the right to the semicircle, the EIS curve enters the Warburg section, which usually exhibits a straight line with 45 degrees from the x-axis and eventually changes to a vertical line at very low frequencies. Here, however, due to the long computational time for further lower frequencies, our simulations stopped at 0.001 Hz, which does not yet reach the capacitance-dominated regime in the Warburg part. 89 Figure 4.7 EIS simulation for half cell with lithium fraction equal to 0.25 in cathode. (a) Applied voltage signal (blue), simulated response current (black circles) at frequency f = 1 Hz and fitted sine function (red). (b) EIS curve with frequencies from 200 to 0.001 Hz. As is observed above, when the amplitude of voltage loading is limited to a small value (here, less than 0.05 V), the entire system is not driven far away from the equilibrium condition. Thus, the current response can still be approximated to a linear response, which has the sinusoidal format similar to the loading voltage. However, when we increase the amplitude of loading voltage (e.g., to 0.5 V), the current no long linearly responds to the loading. Figure 4.8 shows an example of non-linear current response at f = 0.1 Hz for the electrode with average Xp = 0.5. Sharp peak and valley regions are observed on the response current curve. In this simulation, the surface Xp can vary between 0.17 and 0.91, which can lead to a significant material property variation (e.g., Li diffusivity, conductivity, and exchange current density) within the cyclic loading. The resulting current curve cannot be fitted with a single sine/cosine function. The analyses of non-linear EIS require much more simulation results and efforts, which is beyond the scope of our current work. Thus, we will not further discuss non-linear EIS behavior here. However, the presented simulation still demonstrates the potential of using SBM with AMR electrochemical simulations to investigate non-linear EIS behavior in future extensions. 90 Figure 4.8 Non-linear current response due to large loading amplitude: (blue) loading voltage with amplitude of 0.5 V, and (red) current response. The loading frequency is f = 0.1 Hz, which is the same as in Fig. 4.7(a), but the initial Li fraction in this set of simulations is Xp = 0.5. 4.3.3 Effect of Average Li Fraction The intrinsic material properties of cathode particles, such as open circuit voltage (Eq. (4.20)), diffusivity of lithium (Eq. (4.17)), and electric conductivity (Eq. (4.18)), are functions of Li fraction. Even with the same cell voltage imposed on the current collectors, voltage drop across the particle-electrolyte interface can be significantly different due to concentration- dependent conductivity. As in the Butler-Volmer equation, Eq. (4.5), different [ϕ]pe and Cp are expected to result in a substantial difference of the reaction rate on the particle-electrolyte interface. To investigate how electrode Li fractions affect EIS behavior, additional simulations with initial average Xp = 0.50, 0.75, and 0.90 were performed on the same cathode microstructure described in the previous section. Using Eq. (4.20), the values of ϕOCV corresponding to 91 Figure 4.9 (a) The OCV of NMC for different Li fractions used in the simulations. (b) Nyquist plot of simulated EIS curves with Li fractions equal to 0.25, 0.50, 0.75, and 0.90 in the NMC electrode Xp = 0.5, 0.75, and 0.90 are 3.85, 3.71, and 3.45 V, respectively, as shown in Fig. 4.9(a). Those are used to set up the electro-potential boundary conditions on the anode and current collector as the equilibrium cell voltage. The boundary conditions on the remaining four computational domain boundaries are kept the same as in the previous case (Xp = 0.25). Following the same simulation procedure described in the previous section, three additional EIS curves were obtained, which are shown in Fig. 4.9(b). As can be observed in the figure, the radii of the semi-circles increase as the Li fraction increases. This indicates the charge transfer resistance on the particle surface increases when surface Xp increases. The charge transfer resistance at the equilibrium condition can be written as [109] RT Rct = . (4.25) zF i0 Since the amplitudes of applied oscillating loading are sufficiently small, the lithium fraction on the particle surfaces will not significantly deviate from the initial (average) lithium frac- tion. The values of i0 calculated using Eq. (4.21) are 1.448 × 10−1 , 4.079 × 10−3 , 2.685 × 10−3 , and 2.495 × 10−3 mA/cm2 for Xp = 0.25, 0.50, 0.75, and 0.90, respectively. Those i0 values lead to a Rct ratio of approximately 1 : 35 : 53 : 58 for the four cases. This resulting ratio is consistent with that observed from the semi-circles on the Nyquist plot in Fig. 4.9(b), in 92 which the values of Rct can be extrapolated at where the extended semi-circles intersecting the horizontal axis or two times the radii of the semi-circles. Those values are approximately 0.26×106 , 8.5×106 , 12.8×106 , and 14.2×106 Ω. This result demonstrates that the presented electrochemical simulation and impedance calculation methods can properly reflect the input material parameters. Here, the low-frequency results of the three additional simulations do not well extend to the Warburg regime. Since the obtained impedance has already exhibited the effect of i0 variation at different initial Li fraction, we did not pursue further simulations at lower frequencies. 4.3.4 Effect of Electrolyte Salt Concentration As indicated in the Butler-Volmer equation, Eq. (4.5), the salt concentration (Ce ) affects the intercalation/deintercalation reaction rates. To investigate such an effect on the EIS behavior, two more simulations, one with average Ce = 0.5 M and the other with Ce = 2 M, were performed in addition to the case in Section 4.3.2. The ambipolar (De ) and ionic (D+ and D− ) diffusivities are functions of Ce as described in Eq. (4.22). Figure 4.11(a) shows those functions, in which the three diffusivities monotonically decrease as Ce increases. The calculated De values at 0.5, 1, and 2 M are 2.778×10−6 , 1.904×10−6 , and 9.650×10−7 cm2 /s, respectively. As in Eq. (4.23), the shape of κe curve is non-monotonic, with a maximum of 0.02066 S/cm near Ce = 1.2 M. The calculated κe at Ce = 0.5, 1, and 2 M are 0.01429, 0.01959, and 0.01985 S/cm, respectively. Thus, predicting the EIS behavior simply just by De and κe is difficult. Here, because we cannot acquire i0 data measured at Ce = 0.5 and 2 M in the literature, the functions of kf and kb obtained based on Ce = 1 M were still used in this set of simulations. The equilibrium potential, ϕeq , for the two additional simulations were calculated by setting the net reaction to be zero in Eq. (4.5) with Xp = 0.25 (as in Section 4.3.2). The obtained values are 4.151 and 4.187 V, which are utilized to set up the boundary conditions for the equilibrium cell voltage. The dielectric constant (ε) in Eq. (3.3) also varies with Ce . As in Ref. [73], we also simply assumed that ε linearly increases as the 93 Figure 4.10 Calculated κe based on Eq. 4.23 with a maximum point around 1.20 mol/L. salt concentration increases [110]. The calculated specific capacitances using the method in Section 4.3.1 are 107.8 and 190.0 µF/cm2 at Ce = 0.5 and 2 M, respectively. The simulated EIS curves, with that from Section 4.3.2, are plotted in Fig. 4.11(b). Again, 21 loading frequencies ranging from 200 to 0.001 Hz were sampled. Three main observations are summarized from the results. First, the Warburg impedances are almost the same: nearly horizontal translation toward the right on the Nyquist plot as Ce decreases. This indicates that the diffusional impedance is limited by Li diffusion inside the particles because De is six orders of magnitudes larger than Dp . In the low frequency regime, salt concentration can easily reach its equilibrium. The variation of salt diffusivity due to Ce change in the electrolyte has minimal impact on the Warburg part of EIS curve. Second, 94 Figure 4.11 (a) Salt ambipolar diffusivity and ionic diffusivity as functions of salt concentration used in the simulations. (b) Nyquist plot of simulated EIS curves for different salt concentrations in the electrolyte. the radii of the semi-circles significantly increase as the Ce decreases, indicating that the reaction rate sufficiently decreases. This effect is likely attributed to the fact that Ce in the first term of Eq. (4.5) dominates the total reaction rate. Thus, i0 decreases as Ce decreases, which leads to the increase in Rct . Lastly, the critical and minimal frequencies in the three simulations are all approximately at 5 and 0.25 Hz, respectively. There are only slight shifts of their positions on the semi-circles. This manifests that the variations of double-layer capacitance in the salt concentration range examined here has only a negligible impact on the EIS behavior. The specific capacitance varies between 107 and 190 µF/cm2 . The presented simulation framework exhibits the potential of computationally investigating the role of salt concentration on the EIS behavior. However, because variation of Ce simultaneously affects De , κe , ε, and also likely i0 , further study with more material properties is highly expected. 4.3.5 Effect of Microstructures Battery electrodes are agglomerates of active particles for Li storage and non-active additive particles for structural adhesion and electrical percolation. The electrode microstructures 95 are highly complex with tortuous inter-particle space. Many microstructure-related fea- tures, such as active surface areas, particle size distribution, electrical percolation network, tortuosity of inter-particle space, and channel width of electrolyte regions, can all affect the electrode performance. Here, we use the presented simulation tool to explore the effect of microstructures on the resulting EIS. Due to the fact that the time scale for salt diffusion is much smaller than that of loading frequencies examined in this work, as already discussed in Section 4.3.4, the effects of tortuosity of inter-particle space and the effects of channel width will be unable to be examined in the simulated EIS. Nevertheless, the effects of surface areas and diffusion in the particles should be able to be observed in the EIS curves. Another electrode microstructure was generated using the DEM mentioned earlier with 172 spheres, as shown in Fig. 4.12(a). Hereafter, this microstructure is referred to as Ge- ometry 2 (Geo-2) and that in the previous section is referred to as Geometry 1 (Geo-1). The particles in Geo-2 have the same mean radius (5.8 µm) as in Geo-1 on the log-normal distribution, but with a smaller variation (variance: σ 2 = 0.0625). This new electrode is thinner compared to that in Fig. 4.3(c) because the spheres packed more tightly. The parti- cle surface area in this electrode is calculated to be 49151 µm2 , which is approximately 28% less than that of Geo-1 (68061 µm2 ). The surface area is calculated by summing the areas of all triangular patches of the isosurface at ψ = 0.5 generated using MATLAB software. The total volumes of Geo-2 and Geo-1 are 1.733 × 105 and 3.135 × 105 µm3 , respectively. These numbers lead to surface-volume ratios of 0.28 and 0.22 µm−1 in Geo-2 and Geo-1, respectively. The two geometries have a similar porosity of 0.72, calculating in the range of 66 < x < 100 µm. Following the same simulation procedures in Section 4.3.2 with average Xp = 0.25, the EIS curve of Geo-2 is extracted from the simulation results at a series of loading frequencies, as the red curve in Fig. 4.12(b). The charge transfer resistance of Geo-2 is 3.45×105 Ω (twice the semi-circle radius), which is significantly larger that of Geo-1 (2.60×105 Ω). The ratio between Geo-2 and Geo-1 resistances, Rct,2 /Rct,1 = 1.33, is consistent with the inverse of the 96 Figure 4.12 (a) Configuration of the virtual cell for Geo-2. (b) Simulated EIS curves of Geo-2 and Geo-1 electrodes. (c) Configuration of single-particle virtual cell with particle radius of 10 µm. (d) Simulated EIS curves of R2 and R1 cases. ratio between the two electrodes’ surface areas, AGeo1 /AGeo2 = 68061/49151 = 1.38. Because i0 , as a material properties, is the same for the two electrodes, the total reaction fluxes are proportional to the active surface areas. Under the same voltage loading, the resistance is inversely proportional to the current. As a result, the sizes of the semi-circles reflect the active surface areas of the electrodes. Interestingly, the Warburg impedance also approximately follows a similar ratio of Rct,2 /Rct,1 in the frequency range probed. For instance, at f = 0.001 Hz, the magnitudes of the imagi- nary component of the Warburg impedance are 5.62×105 Ω and 3.94×105 Ω for Geo-2 and 97 Goe-1, respectively, which lead to a ratio of 1.43. This observation is counter intuitive to that Geo-2 has a larger surface-volume ratio, meaning a shorter effective diffusion length for inserted/extracted Li. This discrepancy can be attributed to the fact that the diffusional flux in the particles is still limited by the reaction flux. In the frequencies probed in the Warburg regime here, the time scales are still sufficiently close to that of surface reaction. If at an extremely low frequency, where the diffusion time scale completely dominates the EIS behavior, we expect the Warburg impedance would eventually reflect the inverse of the volume ratio (VGeo1 /VGeo2 = 3.135/1.733 = 1.91) between the two electrodes as studied in Ref. [111], instead of the surface ratio. To verify the discussion above, another set of simulations with a simple geometry were performed. The virtual battery cell contains only one single spherical NMC particle: one with radius of R1 = 5 µm and the other with radius of R2 = 10 µm. Figure 4.12(c) shows the configuration of virtual cell of the R2 case. Note that a small region of the particle is truncated at the contact to the current collector. The simulated EIS curves are plotted in Fig. 4.12(d). As inferred from the curves, the semi-circle radii are approximately in a 4 : 1 ratio between the R1 and R2 cases, with 6.1 × 107 and 1.5 × 107 Ω for the charge transfer resistance, respectively. The 4 : 1 ratio is the reverse of the particle surface area ratio between these two cases. Also, the imaginary magnitudes of the Warburg impedance at 0.001 Hz are 1.36×108 and 2.88×107 Ω, which still follow an approximately 4 : 1 ratio but already shows the effect of finite space Warburg behavior. Therefore, the discussions for the Geo-2 and Geo-1 comparison are confirmed by the simple single-particle simulations. 98 4.4 Conclusion In this chapter, an electrochemical simulation framework with explicit considerations of elec- trode microstructures is demonstrated capable of properly capture the concentration and electro-potential evolution under sinusoidal voltage loadings. Through the use of smoothed boundary method, the simulations do not require mesh conforming to complex electrode microstructures. EIS curves are extracted from the physics-based simulations. With double- layer capacitance, surface reaction, and Li bulk diffusion, the obtained impedance curves resemble the typical curves that have a semi-circle in the high-to-intermediate frequency regime and a finite-space Warburg curve in the low frequency regime. Simulations with var- ious initial electrode Li fractions demonstrate the effects of Li-fraction-dependent i0 on the EIS curves. The total charge-transfer resistance is inversely proportional to i0 for the same electrode. This tool is also utilized to investigate the role of initial salt concentration in the electrolyte. The results suggest that i0 variation due to salt concentration dominates the EIS behavior, rather than concentration-dependent salt diffusivity or dielectric constant. Due to the fact that salt diffusivity is several orders of magnitude larger than Li diffusivity in the particles, the Warburg impedance is only determined by the Li bulk diffusion. Lastly, simu- lations with different electrode microstructures demonstrate that the EIS can be employed to probe the active surface areas. For electrodes have the same i0 , the total charge-transfer resistance is inversely proportional to the active surface area. This physics-based simulation framework is expected to be widely used in studying the relationship between EIS behavior, intrinsic material properties, and electrode microstructures. 99 CHAPTER 5 PHYSICAL-BASED SIMULATION OF ELECTROCHEMICAL IMPEDANCE SPECTRUM ON GRAPHITE 100 5.1 Introduction In Chapter 4, we investigate the EIS behavior of cathode electrode. In this chapter, we focus on simulating EIS of anode microstructures. The electrode microstructure will significantly impact a battery’s electrochemical performance. Different from NMC cathode that shows a solid-solution behavior, the anode material studied here (graphite) exhibits multiple phase transformation upon lithiation/delithiation. The common anode material, graphite, has a theoretical gravitational capacity (372 mAh/g [112]) significantly larger than that of cath- odes. Instead of increasing the capacity, research interests of anode materials are placed on increasing rate capabilities [113] or preventing Li plating[114]. Graphite has a layered struc- ture consisting of stacked graphene sheets. In the intercalation processes, Li ions migrate in the space between the graphene sheets, and thus the Li migration shows strong anisotropic behavior: fast in the direction parallel to, but slow in the direction perpendicular to the graphene sheets. The insertion leads to some expansion in the off-plane direction. During lithiation, graphite undergoes three interlayer-ordering phase transition processes: Phase 1’ to 3, 3 to 2, and 2 to 1 [115]. At the fully lithiated state (Phase 1), one Li atom is accommodated by six carbon atoms as LiC6 . Other carbonaceous anode materials include soft carbon and hard carbon. Unlike graphite, soft carbon and hard carbon are amorphous solid. Soft carbon is moderately disordered and can be converted to graphite around 2300 ◦ C. Hard carbon is highly disordered and can hardly be converted to graphite [116]. Because of the amorphous structures, soft and hard carbons exhibit solid-solution behavior during Li insertion processes, leading to a lithiation potentials higher than that of graphite. These three carbon materials all have the advantages of high electronic conductivity and low cost. While graphite’s intrinsic properties are difficult to alter, improvement of electrochemical performance of graphite anode can still be achieved by microstructure designs [117, 118], e.g., introducing tunnels to facilitate Li salt diffusion throughout the electrode. As in the classical electrochemical modelings, the double-layer capacitance is calculated by solving Nernst-Planck-Poisson equations. Li salt transport in the electrolyte is mod- 101 eling by Fickian diffusion. Li insertion/extraction on the particle surfaces is governed by the Butler-Volmer equation. However, we employ the Cahn-Hilliard equation to model Li transport in the graphite particles, as proposed by Guo et al [115] to account for the phase separation in graphite. Gao et al also pointed out that describing Li transport in graphite with phase-separation or solid-solution models will lead to different results [119]. The use of phase field model in the graphite particles differs our work from those using Fickian diffusion for graphite (e.g., as in Ref. [66]). We simulated the EIS curves of an experimentally reconstructed graphite electrode at different stages of lithiation to highlight the effect of Li-fraction dependent properties. The effect of microstructure was examined with simulations on two different experimentally re- constructed graphite electrodes. We also investigated the effect of tunnels in electrodes on the resulting EIS curves. The EIS behavior of two-phase coexisting case was investigated. Interestingly, a unique inductive loop was observed on the simulated EIS curves. The simu- lations also demonstrate that EIS cannot distinguish the difference between core-shell phase coexistence morphology and single phase morphology. Drifting in the EIS simulation was also noticed in some of the simulations, in which the electrodes were not yet relaxed to near equilibrium conditions. Lastly, even though Fick’s diffusion is commonly mistakenly employed to simulate Li transport in graphite, our simulations show that, for an equilibrium system, the Cahn-Hilliard model and Fickian diffusion model produce the same EIS curve. These simulations exhibit a myriad of rich physical phenomena taking place in electrode microstructures, which will not be uncovered without detailed microstructure-level simula- tions. With the advantages of SBM electrochemical simulations, such as speed and ease of implementation, demonstrated in this paper, we expected the presented tool to be widely employed in simulating electrode dynamics to study a variety of electrochemical systems. 102 5.2 Model and Methods Graphite (LiX C6 ) exhibits multiple phase transformations upon lithiation (or delithiation). It undergoes four stages (or phases) from a fully delithiated to a fully lithiated state. At the first stage (1’ where X < 0.06), it can be considered a dilute Li solid solution in graphite with Li atoms sparsely distributed across all inter-graphene layers. These Li atoms are disordered and isolated (nearly noninterative) from each other. As more Li atoms are present in graphite, ordering occurs in the inter-graphene layers and simultaneously a sequence of one filled layer next to two consecutive empty layers appears. This stage is named Stage 3. Further insertion leads to a configuration that one filled layer in every two layers, i.e., Stage 2. At the last stage, every inter-graphene layers are filled with ordered Li atoms: Stage 1 (stoichiometrically denoted as LiC6 ). Different phases can coexist within a graphite particle. The classical Fick’s diffusion model, which is valid only for describing mass transport in solid solutions, cannot properly describe the processes of phase domain propagation/recession. Although the Li ordering within the phase boundaries is not as simple as a monotonic transition, as proposed by Guo et al [115], the phase transformation processes in graphite in the continuum scale can still be modeled by phase-field methods. Since no crystal structure change occurs and mass of Li is conserved in graphite. The Cahn-Hilliard equation will be adequate to model the second-order phase transformation process in graphite, which reads ∂Xp = ∇ · Mp ∇ µp (Xp ) − ε∇2 Xp  (5.1) ∂t where Xp is the Li occupancy fraction, Mp is the transport mobility, µp is the Li chemical potential in the bulk phase, ε is the gradient energy coefficient, and the subscript p denotes graphite particles. Different from typical phase field simulations, we do not use a deter- ministic double-well function for the thermodynamic free energy. Instead, µp in this work is parameterized from measured data in the literature. For a half cell, the Li chemical potential  in the electrode is related to the cell open-circuit voltage (OCV) by ϕOCV = µp − µ0Li /e, 103 where µ0Li is the chemical potential of metallic Li (a constant value), and e is the elemental electron charge. Here, we assume that no concentration gradient is present in the particles. Thus, the Li chemical potential in graphite is obtained by µp (Xp ) = −e · ϕOCV (eV). (5.2) The chemical potential is related to the Gibbs free energy by Figure 5.1 (a) Open circuit voltage (or open circuit potential) as a function of Li fraction in graphite electrodes. (b) Modified chemical potential of Li ions in graphite electrode. The blue dashed curve is the negative of OCV times electron charge. (c) Gibbs free energy of lithium in graphite electrode. The three black dashed lines are the common tangent lines, which indicate two-phase coexistence regions. (d) Exchange current density from kinetic Monte Carlo chronoamperometric simulation (gray dashed curve) and modified exchange current density (red curve) used in the simulations. 104 ∂G(Xp ) µp (Xp ) = , (5.3) ∂Xp where the Gibbs free energy function will have four energy wells (local minima) corresponding to the four stable phases of graphite. Note that we have assumed that metallic Li is the counter-electrode for our simulations and the potential of the reference electrode (Li metal) is zero. Figure 5.1(a) shows the OCV [120] used in this work, which shows three flat plateaus on the curve. Each plateau corresponds to a two-phase coexistence region in graphite. The dark green curve in Fig. 5.1(b) shows µp parameterized from the ϕOCV . Because ϕOCV in the two- phase coexistence regions are flat plateaus, we extrapolated µp to be non-monotonic in the three two-phase regions (i.e., 0.05 < Xp < 0.12, 0.23 < Xp < 0.51, and 0.56 < Xp < 0.97). Without those non-monotonic functions in the two-phase regions, no phase separation will occur. The extrapolated curve leads to a nucleation barrier about 15–40 e·mV. Here, we constructed the non-monotonic regions of the µp curve simply for numerical ease of phase- field simulations such that stable phase boundaries spanning about 4 grid spacings can form inside graphite particles. Some measured intrinsic voltage hysteresis on graphite OCV is reported to be approximately 10–40 e·mV [121, 122, 123]. The voltage hysteresis is equivalent to the barrier mentioned above. Since our values do not deviate much from the reported one, the simulation results using the parameterized µp should be acceptable. Integrating µp function, we obtained the Gibbs free energy as the light green curve in Fig. 5.1(c), in which there are four local energy wells corresponding to the four stbale phases and the three common tangent (dashed) lines indicate the three two-phase regions. Again, since we do not have access to the interfacial energy of the phase boundaries in graphite, the value of ε in Eq. (5.1) was chosen for numerical convenience, such that the the phase boundaries span approximately 4 root-level grid spacings.  The Li mobility in Eq. (5.4) is parameterized using the relation: Mp = Dp / ∂µp /∂Xp , where µp is shown as the dark green curve in Fig. 5.1(b). The composition-dependent mobility is plotted as the red curve shown in Fig. 5.2(a). As in Refs. [124, 125], we assume four constant values of Dp in the four respective single-phase regions (see Fig. 5.2(b)) and use the 105 Figure 5.2 (a) Li mobility in graphite used in the simulations. (b) Diffusivity values in the four single phases (gray solid lines) and linear interpolation in the two-phase regions (gray dashed lines). gradient of chemical potential function to estimate the mobility in the single phase regions. The mobility in the two-phase regions (the humps on the red curve) is extrapolated from the single-phase regions (the valleys on the red curve). Although Li transport in graphite is anisotropic, we used isotropic Mp in this work. This is justified based on Ref. [125], in which the simulations showed only a negligible difference between the cell performance results of anisotropic and isotropic models if the isotropic model uses an effective mobility averaged over the anisotropic mobility. We used a continuous domain parameter, ψp , to define graphite particles: ψp = 1 in the particles and ψp = 0 outside the particles. The particle-electrolyte interfaces are the regions where ψp varies from 1 to 0. Note that this diffuse interface defined by ψp is a mathematically smeared boundary, not a physical interface. The domain parameter was used to reformulate Eq. (5.1) to ∂Xp 1 1 |∇ψp | ∇ · ψp Mp ∇ µp (X) − ε∇2 Xp +  = rxn (5.4) ∂t ψp ψp ρ where ρ is Li site density in graphite (0.0312 mol/cm3 ), rxn is the (de)intercalation rate on the particle surfaces, and |∇ψp | indicates the particle surfaces. Because the particle regions 106 are defined by ψp function, Eq. (5.4) can be solved on numerical mesh non-conformal to the particles. Thus, the implementation of simulations can be performance much easier than the conventional sharp-interface methods. As shown in Chapter 4, the rest of SBM electrochemical governing equations in the simulations are provided here again for readers’ convenience, ∇ · (ψp κp ∇ϕp ) − |∇ψp |z+ F rxn = 0, (4.10) ∂Ce 1 |∇ψe | rxn t− ie · ∇t+ = ∇ · (ψe De ∇Ce ) + − , (4.12) ∂t ψe ψe ν+ z+ ν+ F rxn ∇ · [ψe (z+ m+ − z− m− ) F Ce ∇ϕe ] + |∇ψe | = ∇ · [ψe (D− − D+ ) ∇Ce ] , (4.13) ν+ where κp is the electric conductivity of graphite, ϕp is the electrostatic potential in graphite, Ce is the salt concentration in the electrolyte, De is the ambipolar diffusivity of the salt, ψe is the domain parameter for the electrolyte region, ti is the transference number, νi is the dissolution number, zi is the valence number, ie is the ionic current, F is the Faraday constant, mi is the ionic mobility, and Di is the ionic diffusivity. The subscripts + and − denote cations and anions, respectively. These equations have appeared in Chapter 4, but we provide them again here for readers’ convenience. Here, ψe = 1 − ψp . The material parameters in Eqs. (5.4) through (4.13) can be found in Refs.[93, 124, 125]. The surface (de)intercalation reaction (charge transfer process) is governed by the Butler- Volmer equation:     −αz+ F p (1 − α) z+ F p rxn = kf Ce exp [ϕ]e − kb Cp exp [ϕ]e , (4.5) RT RT where kf and kb are the forward and backward reaction rate constants, α is the symmetry factor, R is the ideal gas constant, T is the absolute temperature, Cp is the Li concentration (Cp = ρXp ), and [ϕ]pe is the electrostatic potential drop across particle-electrolyte interface. The rate constants can be calculated from the exchange current density as in Refs. [93, 124]. The exchange current density can be obtained by performing chronoamperometric simulations at different electrode potentials and letting the system approach towards the 107 steady state. At steady state, the net Li ion flow across the interface becomes zero, i.e., the oxidation current equals the reduction current. The gray dashed line in Fig. 5.1(d) shows an exchange current density curve obtained by Gavilan-Arriazu et al using kinetic Monte Carlo chronoampreometric simulations. This curve exhibits a nearly flat plateau in the region 0.2 < Xp < 0.8 and the value decreases towards the two ends of the horizontal axis. If this i0 curve were used in the EIS simulations, the charge transfer resistance would be almost the same within the entire plateau region. Thus, to highlight the effect of i0 on the EIS simulations, we modified the gray curve to the red curve in Fig. 5.1(d). As mentioned earlier, there are two stable single-phase regions of Li in graphite (Stages 3 and 2) within the plateau on the gray dashed curve. These single phase regions are thermodynamically stable. Thus, we hypothesize that the respective i0 values in the single phase regions are smaller than those in the two-phase regions. However, there is no available literature data for the parameterization. Here, we selected i0 values to be approximately 1.5 and 1 mA/cm2 for Stages 3 (Xp ∼ 0.19) and 2 (Xp ∼ 0.54), respectively. These values are moderately larger than those in Stages 1’ and 1, but enough to result in different EIS curves in the simulations. Here, we emphasize that the i0 values for Phases 3 and 2 are chosen only to contrast the differences in exchange current density in different phases. We acknowledge that the red curve in Fig. 5.1(d) was constructed for the ease of demonstrating the EIS simulations, which does not represent the real values of i0 on graphite surfaces. If available literature data exist, a more reliable i0 function could be parameterized for the simulations. We used experimentally reconstructed 3D graphite electrode microstructures in the sim- ulations. These microstructures were obtained using X-ray computed nano-tomography and are publicly accessible [126, 127]. Electrode II on the online repository consists spherical par- ticles and Electrode IV has flaky particles. We selected these two microstructures to examine the effect of microstructures on the EIS behavior. The voxel centers in the 3D microstructure data arrays were directly used as the root-level grid points in the SBM simulations. With the voxel size of 0.325 µm, the root-level grid spacing is ∆x = 0.325 µm. 108 To increase the simulation accuracy, octree adaptive mesh refinement (AMR) [93, 72] was performed to generate the mesh systems, in which fine mesh was used near the dif- fuse interfaces while coarse mesh was used in the particle and electrolyte bulk regions. Figure 5.3 Experimentally reconstructed 3D graphite electrode microstructures: (a) Electrode II with spherical particles, (b) Electrode IV with flaky particles, (c) Electrode II with a tunnel at the center, and (d) Electrode IV with a tunnel at the center. With one-level of refinement, the diffuse interface thickness was controlled to be approx- imate 2 root-level grid spacings. We solved the electrochemical governing equations, Eq. (5.4) through Eq. (4.5), to obtain the reaction current responding to the oscillating voltage 109 loadings. The specific capacitance stemming from double layer formation on the particle surfaces was calculated separately by solving the Nernst-Planck-Poisson (NPP) equations as in Refs. [73, 128]. The total response current consists of reactive and capacitive currents. The phase angle and amplitude of the impedance was fitted using MATLAB as in Ref. [128] and in Chapter 4. 110 5.3 Results and Discussion 5.3.1 Electrochemical Simulation of Graphite Lithiation Two graphite electrode microstructures [126] were used as the input geometries for the simulations. We cropped a region of 180×120×180 voxels out of the data, which corresponds to a region of 58.5×39×58.5 µm3 . Fig. 5.3(a) and (b) show the microstructures taken from Electrode II and Electrode IV, respectively, from the online data. Hereafter, these two microstructures used in the simulations are referred to as E II and E IV in this text. The total particle surface areas are approximately 3.54 × 10−4 and 3.67 × 10−4 cm2 for E II and E IV, respectively. The values were obtained by summing all triangular patches on the isosurface of ψp = 0.5 generated using MATLAB. An empty space with a thickness Figure 5.4 The virtual cell used in the simulations. The yellow color indicates graphite particles, the semitransparent cyan slab indicates the Li foil anode, and the brown slab indicates the current collector. The transparent region is filled with electrolyte. The space between Li foil and graphite electrode serves as the separator. of 50 root-level ∆x (16.25 µm) was included between the Li foil and the electrode as the separator layer. Figure 5.4 shows the virtual cell for E II simulations, for which Eqs. (5.4) 111 and (4.9) are solved in the particles (yellow region), Eqs. (4.12) and (4.13) are solved in the electrolyte (transparent region), and Eq. (4.5) is calculated on the particle surfaces. The semi-transparent cyan and solid brown slabs represent the Li foil and the current collector, respectively. We assume an EC/EMC electrolyte with 1 M of LiPF6 salt dissolved within. The ionic diffusivities were set to be 1.25×10−6 and 4.0×10−6 cm2 /s for D+ and D− , respec- tively, at 1 M as in Ref. [99]. The electrostatic potential on the Li foil was set to be constant (ϕe |a = 0 V). The electrostatic potential on the current collect (ϕp |c ) was constantly adjusted to maintain a constant current loading for lithiation (or delithiation). Figure 5.5 Simulated Li fraction in graphite electrode (a) at average XLi = 0.54 and (b) XLi = 0.85. Clear multi-phase coexistence in graphite particles can be seen in (a) and (b). Simulated (c) electrostatic potential in graphite particles, (d) salt concentration in electrolyte, and (e) electrostatic potential in electrolyte corresponding to (a). (f) Simulated cell voltage curve at 0.5C lithiation of E II. The circular and square markers correspond to (a) and (b), respectively. Figures 5.5(a) and (b) show simulated Li fraction in E II at two different times during a 0.5C (or C/2) lithiation. 1C (or C/1) means a rate that completely delithiates or lithiates the whole electrode in 1 hour. As Li is inserted into the particles, Li fraction increases near 112 the surface regions and the particle cores have a lower Li fraction. Clear phase coexistence in graphite particles can be observed. For example, Stage 2 (green) and Stage 3 (light blue) in Fig. 5.5(a). See the cyan dashed oval. In Fig. 5.5(b), three phase coexistence can be distinguished by the yellow, green, and blue colors. See the red dashed ovals. Other physical fields, i.e., electrostatic potential in the particles, electrostatic potential in the electrolyte, and salt concentration in the electrolyte are also solved in the simulations. They are plotted in Fig. 5.5(c), (d) and (e), respectively. The cell voltage curve during the lithiation process is plotted in Fig. 5.5(f). In a slow lithiation process, like 0.5C here, the cell voltage curve exhibits multiple steps as graphite’s OCV curve. These steps indicates lithiation via two- phase processes. The details of electrochemical simulations of charge/discharge cycling of graphite electrodes can be found in Ref. [125]. E IV has flaky graphite particles with a slightly larger total particle surface area than that of E II. Because the plate-like particles are perpendicular to the primary electrolyte diffusion direction (the electrode thickness direction), the tortuous channels for electrolyte will hinder Li ion migration through the electrolyte at high lithiation rates, thus leading to a lower cell performance of E IV. Figure 5.6(a) shows the cell voltage (CV) curves of E IV and E II at 6C lithiation. Compared to the CV curve of 0.5C lithiation, the 6C CV (the red curve in Fig. 5.6(a)) drops much quicker. It reaches 0 V near XLi ∼ 0.6, which means only 60% of graphite material is utilized up to the cutoff point. The details of cycling simulations can be found in Ref. [125]. The E IV CV curve is below the E II CV curve, indicating E IV has a poorer performance compared to E II. Here, because the electrode is only 39-µm thick, the hindrance of cell performance due to pore tortuosity is not prominent. In thicker electrodes, the performance of E IV will be much deteriorated. Tunnels have been introduced in cathodes [129] and anodes [130, 117, 118] to enhance the ion transports in electrolyte at high-rate operations. The effect of including tunnels has also been studied in simulations [66]. Here, we present the simulated electrochemical performance of the two electrodes, each of which has a tunnel through the middle. Figures 5.3(c) and (d) 113 Figure 5.6 Simulated cell voltage curve at 6C lithiation of (a) E II and E IV, and (b) E IV and E IV with tunnel. show the microstructures of E II and E IV, respectively, with a straight tunnel each along the thickness direction. The tunnel diameter is 13 µm. The two microstructures are referred to as E II T and E IV T. The simulated CV curve of E IV T at 6C lithiation is plotted as the black curve in Fig. 5.6(b), with E IV CV curve (red) in the same figure for comparison. The E IV T cell voltage curve is above the E IV one, indicating that E IV T has a better rate performance due to the enhancement of electrolyte transport through the tunnel. Again, the effect of enhancement will be much magnified in thicker electrodes. The case of E II T versus E II also shows an improvement of rate performance, but in a much smaller degree relative to the E IV T-to-E IV case. The tunnel did not significantly enhance the long range ion transport in E II T. This is because the pore tortuosity in E II (∼1.6) is significantly smaller than that in E IV (∼2.6). The E II T 6C cell voltage curve is not presented here to avoid repetitive descriptions. In a brief summary, an electrode with less tortuous pore channels has a better high rate performance (e.g., E I) versus E IV). An electrode with a tunnel that can facilitate salt ion transport throughout will also has a better high rate performance (e.g., E IV T versus E IV). 114 5.3.2 EIS of Single Phase Graphite The same electrochemical simulation code was employed to simulate EIS processes. While the electrostatic potential on the Li foil was still set to be constant (ϕe |a = 0 V), an oscillating voltage load was imposed by setting ϕp |c = V0 sin(ωt) on the current collector, where ω = 2πf was the angular frequency, f was the ordinary frequency, and the amplitude was V0 = 15 mV. The specific capacitance (Csp ), originating from double-layer formation, was 140 µF/cm2 [128]. The details of simulation implementations can be found in Refs. [93, 128] and in Chapter 4. The total response current under the oscillating voltage loading includes reactive and capacitive currents as itot = irxn + idl (5.5) where irxn represents reactive current and is obtained from the Butler-Volmer reaction, Eq. (4.5). The capacitive current idl originates from double layer formation and is calculated R  using idl = Csp ∆V /∆t dA [128] as in Chapter 4, where ∆V is the voltage drop across particle-electrolyte interfaces, ∆t is the time step size, and A is the particle surface area. These currents were integrated over the entire electrode region. The total current was fitted to a sinusoidal function using MATLAB curve fitting function. An example of voltage loading and total current response is shown on Fig. 5.7, in which the Li fraction is 0.54 and frequency is 32 Hz. Finally, the impedance was calculated as Z(f ) = ϕp |c /itot . We examine EIS behavior over a frequency range from 512 to 0.01 Hz. Again, the details of the implementations can be found in Ref. [128] and in Chapter 4. The first set of simulations were performed to examine graphite electrodes with stable, uniform single phases, i.e., Stage 1’, 3, 2, and 1. The Xp was set to be uniformly 0.02, 0.19, 0.54, and 1 throughout all particles in E II. Hereafter, we use XLi and Xp interchangeably since both of them are Li fraction in graphite particles. The corresponding equilibrium cell voltage are 0.694, 0.159, 0.109, and 0.0 V for these four simulations. These values were selected based on the OCV at the respective XLi . Figure 5.8(a) shows the simulated EIS 115 Figure 5.7 The blue curve is the voltage loading imposed on the half cell, the black circles are the current response at the sampling timings, the red line is the fitted current response curve. There is a clear phase shift between the loading and response. curves with magnified views in Fig. 5.8(b). All the curves exhibit a well-developed semicircle in the intermediate-to-high frequency range: 2.0–512 Hz. The four curves (black, green, blue, and red) all show recognizable Warburg impedance: a line 45 degrees from the horizontal axis. On the blue curve (XLi = 0.54), the critical frequency is between 32 and 64 Hz. This value agrees with that estimated from the total charge-transfer resistance and total capacitance. The critical frequency can be analytically expressed as fc = 1/(Rct Ctot )/(2π), which for E II is 43.97 Hz. The total capacitance for E II would be Ctot = Csp A = 4.96 × 10−8 F. At XLi = 0.54, i0 = 1 mA/cm2 , and the total charge-transfer resistance [109] for E II would be RT Rct = , (5.6) zF i0 A which leads to 7.29 × 104 Ω. The Rct calculated from i0 and A is close to the diameter (7.16×104 Ω) of the blue semicircle. Note that the semicircles do not intersect the horizontal axis because of the overlap with Warburg impedance. Hereafter, E II with a uniform Stage 2 116 Figure 5.8 EIS curves extracted from 3D graphite microstructure simulations. (a) EIS curves for single-phase Stage 1’, 3, 2, and 1 on E II electrode. (b) Magnified view of (a). (c) Simulated EIS curves for E II, E IV, E II with tunnel, and E IV with tunnel. (d) Magnified view of (c). phase (XLi = 0.54) is selected to be the baseline case in this chapter unless otherwise stated. The cross-sectional area of this electrode is 3.42 × 10−5 cm2 , which leads to a cross-sectional resistance of the electrode to be 2.45 = (7.16 × 104 ) × (3.42 × 10−5 ) Ω · cm2 . Such a unit may be more commonly used in the battery community. However, for the clarity of expressing results, we still use the total resistance in this chapter. For Stage 3 and Stage 1, the values of Rct calculated based on i0 and A are 4.86 × 104 and 1.82 × 105 Ω, respectively. The diameters of the green and red semicircles are 4.68×104 and 1.75×105 Ω, respectively, which agree with 117 those values estimated from i0 and total surface area. These semicircle radii were obtained from circle fitting using MATLAB. (Only the data points on the semicircle are used to fit into the circle function. See Fig. 5.9 for an example. The diameter of the circle represents the charge transfer resistance). The obtained values are slightly smaller than those directly predicted from i0 and A. The difference could originate from a slight overestimate of the total particle surface area in generating the isosurface for visualization in MATLAB. Figure 5.9 The red semicircle is fitted from the calculated impedance data points of the baseline case (with XLi = 0.54 in E II). The diameter is 7.16 × 104 Ω. For Stage 1’ (XLi = 0.02), the semicircle diameter (9.18 × 105 Ω) also agrees with the estimated Rct (9.35 × 105 Ω). However, a significant overlap with the Warburg part can be observed on the black curve, which is caused by the shift of Warburg part to the left. For a typical resistor-capacitor-Warburg circuits, increasing diffusivity will translate the Warburg part toward the origin [131]. The diffusivity used in parameterizing Li mobility at XLi = 0.02 is more than one order of magnitude larger than those in Stages 3 through 1. This explains the large overlap between the semicircle and Warburg part, as well as the extension of 118 Warburg part toward the upper right corner, on the black curve. These results demonstrate that the SBM electrochemical microstructure simulations reliably capture the EIS behavior of single-phase graphite electrodes qualitatively and quantitatively. Although E IV has a poorer performance, its simulated EIS curve (the red curve in Fig. 5.8(c)) shows that E IV has a smaller total resistance than E II (the blue curve in Fig. 5.8(c)). The red curve in Fig. 5.8(c) was obtained from simulations with XLi = 0.54 (Stage 2) uniformly throughout all particles in E IV. Figure 5.8(d) shows the magnified view of Fig. 5.8(c) near the region corresponding to the minimum frequencies. The Rct (semicircle diameter 6.68×104 Ω) estimated from simulated Nyquist plot agrees with the value (7.04×104 Ω) calculated based on the i0 and particle surface area. The smaller resistance of E IV stems from the fact that E IV has a larger active surface area than that of E II. Since EIS is conducted in near equilibrium conditions, it can only probe properties at near equilibrium conditions. High C rate performances, which in this case are limited by Li-ion transport in the electrolyte, cannot be reflected in the EIS results. As in the comparison between E IV and E II, EIS behavior cannot reflect the improve- ment of high-rate performance enabled by tunnels. The simulated EIS curve of E IV T is shown as the black curve in Fig. 5.8(c). Although E IV T has a better rate performance, its total resistance is larger than that of E IV. The semicircle diameter (7.06 × 104 Ω) of the black curve is slightly larger than that of the red curve. The diameter of the black semicircle agrees with the value (7.23 × 104 Ω) calculated using i0 and the active surface area (3.58 × 10−4 cm2 ). While E II T has a similar 6C rate performance to E II, its EIS curve has a larger diameter (7.34 × 104 Ω) than E II. See the green curves in Figs. 5.8(c) and (d). The tunnel in E II T results in a decrease of total particle surface area from 3.54 × 10−4 to 3.46 × 10−4 cm2 , which leads to the calculated Rct for E II T to be 7.47 × 104 Ω. Again, the resistance (semicircle diameter 7.34 × 104 Ω) of the simulated EIS curve agrees with that value. While the simulations above were all controlled to confine the particle XLi variations 119 Figure 5.10 Simulated EIS curve with a large loading amplitude (V0 = 30 mV). The Warburg region exhibits a high capacitive response. within their respective single-phase regions, one additional set of simulations were conducted to examine the scenario if XLi entered a two-phase region (miscibility gap between two stable phases). The simulations were performed with uniform XLi = 0.54 (Stage 2) throughout the graphite particles, but the voltage amplitude was much larger (V0 = 30 mV). Figure 5.10 shows the obtained EIS curve. The semicircle portion of the curve is very similar to that of the baseline E II case (the blue curve in Fig. 5.8(a)). However, the Warburg region exhibits a strong finite-space Warburg (FSW) capacitive impedance, as shown by a nearly vertical line at low frequencies, which indicates a low effective diffusivity. (This phenomenon under large loading amplitudes has been similarly observed in other continuum-scale simulation results [132] for mass-transport (Warburg) impedance.) The strong FSW behavior is very different from the baseline case (with V0 = 15 mV) even though the same material properties were used in the simulations. We attribute this phenomenon to the fact that the particle surface XLi exists the single-phase region (solubility limit) and enters the two-phase region 120 Figure 5.11 Snapshots of XLi in graphite particles. (a)–(c) Under a large loading amplitude (V0 = 30 mV) and (d)–(f) under a small loading amplitude (V0 = 15 mV). The surface XLi in (b) enters the miscibility gap between Stage-2 and Stage-1 phases. All surface XLi under a small loading amplitude remains within the solubility limit of Stage-2 phase. (a) through (f) use the same color bar. under the large loading amplitude. Diffusivity is related to mobility by D = M (∂µ/∂X). The thermodynamic factor ∂µ/∂X is usually small in a two-phase region. As a result, the effective diffusivity is smaller in the two-phase regions (miscibility gap) than in the single-phase regions (if modeled as a phase-separating material). This behavior of significant decrease of effective diffusivity in two-phase regions has been observed in many first-order phase-transformation materials [133, 134, 135]. In experimental observations, this would be understood as a sluggish phase boundary motion. Fig. 5.11(a) through (c) show snapshots of XLi in the graphite particles under the larger loading amplitude. At t = 4 s, the particle surface XLi (> 0.58) is above the solubility limit of Stage 2 (0.56), which supports our explanation of the strong FSW impedance behavior under large loading amplitude. However, we should note that there was still no phase boundary present on the particle surfaces. The 121 hindering of Li insertion/extraction was along the direction normal to the particle surfaces. Figures 5.11(d) through (f) show snapshots of XLi in the baseline case (V0 = 15 mV) for comparison, in which particle surface XLi remained within the solubility limit throughout the entire simulation. In these single-phase test cases, the SBM simulations well predict the EIS curves of the graphite electrodes. As EIS can only probe properties at near equilibrium conditions, the effect of transport in electrolyte, which is the limiting factor of the high-rate performance, cannot be reflected on the EIS behavior. Additionally, a low-tortuosity structure can usually has a smaller active surface area that results in a larger total Rct of the electrode. Thus, it could mislead to a counter-intuitive scenario that an electrode has a high Rct on the EIS but has a better high-rate performance. 5.3.3 EIS of Multi-phase Graphite It is common for graphite to stay in a two- (or possibly more) phase coexistence. Thus, the electrochemical simulations were used to examine the EIS behavior of graphite with a two- phase coexistence. To highlight the role of phase boundaries on the EIS behavior, first we examined a hypothetical phase morphology of spindoal decomposed phase distribution. We used the Cahn-Hilliard equation to generate an initial configuration that contained Stage-2 and Stage-1 phases in the particles in the E II microstructure. The composition was initially XLi = 0.75 with a small noise throughout the graphite particles while all other governing equations were switched off (i.e., no surface reactions involved). Spinodal decomposition occurred within the graphite particles, leading to phase separation and resulting in regions of XLi ∼0.56 and ∼0.98 as shown in Fig. 5.12(a). Note that the spinodal decomposition here is a hypothetical scenario that could occur when a lithiated graphite electrode is first heated to become a random solid solution of Li and then tempered down to the room temperature during which phase separation takes place. As discussed later, phase morphology resulting from typical lithiation/delithiation would exhibit a core-shell distribution, which is different 122 Figure 5.12 (a) Phase morphology in graphite particles generated using the Cahn-Hilliard equation. The equilibrium is established between Stage-2 and Stage-1 phases. (b) Simulated EIS curve for the configuration in (a). The curve exhibits an inductive loop. (c) The Bode plot corresponding to (b), which shows an decrease in the amplitude at the frequencies corresponding to the inductive loop. from the spinodal morphology. However, simulating the EIS curve of a spinodal morphology provides a potential way to identify phase separation behavior of graphite electrodes. The volumes occupied by the two phases were approximately equal. The simulated EIS curve is shown as the black curve in Fig. 5.12(b). Interestingly, this EIS curve exhibits a ‘loop’ near 123 the minimum frequency region. The curve bends inward into the semicircle when f < 4 Hz. As the frequency continued to decrease, the curve moved outward following a Warburg type curve when f < 0.1 Hz. This type of loops are commonly referred to as inductive loops, which have been observed in many EIS measurements involving phase transformations such as magnetization [136] and corrosion [137]. Recently, a chemical inductance mechanism has been introduced to explain the low-frequency inductive loop [138]. However, to the best of our knowledge, we are unaware of reported inductive loops on graphite electrode measurements. Inductive loops on battery EIS curves are usually in the high-frequency regime and are stemming from the induction of metal wires or metal casings. If inductive loops appear in a low-frequency regime, they are commonly attributed to the formation of solid electrolyte interphases (SEI). Nevertheless, we did not include SEI formation in this model. Here, we hypothesize that the low-frequency loop originates from the following mecha- nism. As can be seen in Fig. 5.12(a), the boundaries between Stage-2 and Stage-1 phases intersect the particle surfaces. Thus, the particle surfaces can be viewed as a parallel con- nection of Stage 2, Stage 1, and phase boundaries, in which the phase boundaries region (0.58 < XLi < 0.96) have a larger i0 value (∼2.77 mA/cm2 ) than the remaining two types of surfaces (∼1.9 and ∼0.72 mA/cm2 ). See the red i0 curve in Fig. 5.1(d). The semicircle radius reflects the average Rct of those three types of surfaces, thus being smaller (i.e., lower resistance) compared to those individually in the Stage-2 or Stage-1 cases. Near the mini- mum frequency regime, the phase boundaries dominated the response, which led to a smaller overall Rct . Figure 5.12(c) shows the Bode plot of the impedance response. A local minimum of impedance can be observed in the frequency range of 0.01 < f < 1 Hz. This frequency range corresponds to where the inductive loop appears on the Nyquist plot. As a result of the decreased impedance, the EIS curve bent inward. Further decrease in frequency led to the regime where the Li transport in the bulk of the particles dominated the impedance. Thus, the curve followed the typical Warburg line again. 124 Figure 5.13 Phase morphologies of Stage 3-2 coexistence generated using (a) the Cahn-Hilliard equation and (b) planar case. (a) and (b) use the same color bar. (c) Simulated EIS curves for the two Stage 3-2 cases, which also show inductive loops. To verify this hypothesis of inductive loop formation, another series of electrochemical simulations were performed on a configuration that contained coexisting Stage-3 and Stage-2 phases. The configuration was generated in a similar manner but with an average XLi around 0.375. The spinodal decomposition led to regions of XLi ∼0.26 (Stage 3) and ∼0.51 (Stage 2), as shown in Fig. 5.13(a). The volumes occupied by these two phases were similar. The simulated EIS curve is plotted as the blue curve in Fig. 5.13(c). Again, an inductive loop appears on the EIS curve. Yet, the loop is much smaller compared to that in the previous Stage 2-1 coexistence case. This may be due to the fact that the i0 value (∼2.75 mA/cm2 ) in the Stage 3-2 phase boundaries is less deviated from those in the Stage-3 (∼2.6 mA/cm2 ) or Stage-2 (∼1.6 mA/cm2 ) phases, compared to the previous Stage 2-1 case. Note that if i0 in the phase boundaries were the same as in those bulk phases, the loop would not appear and the EIS curve would exhibit a typical semicircle with Warburg part as in those single- phase cases. Here, the EIS curve of the Stage 2-1 case is also plotted as the black curve in Fig. 5.13(c) for comparison. The semicircle radius of the blue curve is smaller than that of the black curve. As discussed earlier, the radius reflects an average of Rct on the particle surfaces. The Stage 3-2 equilibrium here is established between XLi ∼0.26 and ∼0.51. (The Stage 2-1 equilibrium is between XLi ∼0.56 and ∼0.98.) The average i0 value in the Stage 3-2 case is higher than that in the Stage 2-1 case. Therefore, the semicircle radius of the 125 Stage 3-2 case is smaller. Another series of simulations were performed on an artificially created configuration, in which the region (x < 37 µm) close to the separator was assigned to be XLi = 0.25, while the rest of the electrode was assigned to be XLi = 0.5. The volumes of the two phases are approximately the same. See Fig. 5.13(b) for this configuration. The simulated EIS curve is plotted as the red curve in Fig. 5.13(c), which also shows an inductive loop. The red curve is overall similar to the blue curve but with a slightly smaller radius. This is expected because the average i0 value based on the XLi values here will be larger than that in the previous State 3-2 coexistence case. In this planar phase morphology case, the values of i0 are ∼2.48 and ∼2.5 mA/cm2 at XLi = 0.25 and 0.5, respectively, which are larger than those values in the previous Stage 3-2 case. We acknowledge that a recent Cahn-Hilliard phase-field simulation of a LiFePO4 single particle demonstrated strong inductive behavior on the EIS curve when phase boundaries are present on the particle surface [139]. The miscibility gap of Lix FePO4 is much wider than that of graphite. The inductance observed in Ref. [139] is much larger than those in our work. The authors attributed the inductance to the hysteresis arising from a high lithiation flux across the particle surface but a slow phase-boundary motion in the particle. In our simulations, the enhanced insertion rate is due to the increase i0 at the phase boundary present on the particle surfaces. Nevertheless, Ref. [139] and our work both demonstrate that, for phase- separating materials, inductance could occur in the EIS processes when phase-boundaries intersect with the active particle surfaces. This highlights the importance of physics-based 3D microstructure simulations because intersections between phase boundaries and particle surfaces can only be resolved in particle-level simulations. The conventional PET simulations are difficult to capture such detailed dynamics because they cannot consider concentration or phase distribution other than along the particle radial direction. Thus, phase boundaries cannot intersect spherical particle surfaces. Li is inserted/extracted via the particle surfaces during lithiation/delithiation processes. 126 Figure 5.14 Stage 3-2 core-shell phase morphologies of (a) thick shell and (b) thin shell. The shell layers have a XLi similar to that of Stage 2. (a) and (b) use the same color bar. (c) Simulated EIS curves for the two core-shell Stage 3-2 cases, which are similar to the single-phase Stage 2 case. The black curve in covered underneath the red curve. As shown in Section 5.2, XLi gradient will develop along the particle radial direction. A concentric core-shell phase coexistence is expected [115, 119, 125, 124] in this case. Since most EIS measurements are performed on electrodes that have been charged or discharged, we anticipate that a core-shell phase coexistence is a likely configuration of graphite electrode particles in typical EIS measurements. Thus, we employed the electrochemical simulations to examine the EIS behavior of core-shell phase coexistence. Figure 5.14(a) shows a manually created phase morphology configuration, in which the particle shell regions (within 3.25 µm to the surfaces) were assigned to be XLi = 0.54 (Stage-3 phase). The core regions were assigned to be XLi = 0.19 (Stage-2 phase). In this scenario, while the average Li fraction over the electrode was 0.507, all particle surfaces were covered by Stage-3 phase. The simulated EIS curve is plotted as the red curve in Fig. 5.14(c). The EIS curve of the baseline case (uniformly XLi = 0.54 throughout the electrode) is provided as the blue curve for comparison. The red and blue curves mostly overlap, except for the low-frequency region. Because the surface XLi in this core-shell configuration is the same as in the baseline case and that double layer capacitance and charge-transfer reaction occur only on the particle surfaces, the semicircles of the red and blue curves overlap. At a low frequency where Li transport in the particles dominated the response, the presence of phase boundaries in the 127 particles slightly increased the Warburg impedance, thus leading to the small difference between the red and blue curves. Note that no inductive loop appears on the EIS curve of the core-shell two-phase coexistence case because no phase boundaries intersect the particle surfaces. The configuration in Fig. 5.14(b) is another phase morphology created in the same man- ner but with a thinner shell (within 1.95 µm to the surfaces) of Stage-3 phase. The average Li fraction is 0.443. Again, all particle surfaces were covered with Stage-3 phase. The black curve in Fig. 5.14(c) is the simulated EIS curve. The black curve closely overlaps the red one and covered underneath the red curve. These concentric core-shell simulations demon- strate that even though the average Li fraction (equivalently the degree of discharge (DOD) or the state of charge (SOC)) of the electrode is substantially different, the resulting EIS curves are indistinguishable if the shell layers have the same XLi . This means that an EIS measured property will vary step-wisely versus the electrode’s DOD (or SOC). For example, the semicircle diameter remains a constant (due to the same surface XLi ) for a range of DOD, then followed by a rapid change (due to formation of a new phases of the shell layer). Therefore, if probed surface properties remain the same over a significant range of an elec- trode’s DOD (or SOC), that electrode could be undergoing a core-shell two-phase lithiation process. Conversely, if the semicircle diameter varies continuously versus the DOD, this is an indication that the material is undergoing a solid-solution lithiation. See a schematic il- lustration in Fig 5.15. Furthermore, as EIS measurements are usually performed on graphite particles with core-shell phase morphologies, the low-frequency inductive loops seen in the previous simulations are unlikely to be observed in experiments. The phase morphology gen- erated from spinodal decomposition (using the Cahn-Hilliard equation) probably only exists in the scenario of annealing partially lithiated graphite electrode. This may explain why low-frequency inductive loops have never been reported for graphite electrodes (in addition to the possibility that the manually constructed i0 function in this work could be overly de- viated from the physical one). Nevertheless, this work proposed a method to experimentally 128 Figure 5.15 A schematic illustration of step-wise variation of EIS measured properties versus electrode DOD. The red curve is continuously varying, which is an indication of a solid-solution material. verify whether spinodal decomposition will occur by EIS measurements. 5.3.4 EIS of Graphite Immediately after Lithiation Inheriting from charging/discharging processes, XLi gradients form along the particle radial direction. Moreover, degrees of lithiation/delithiation can vary heterogeneously in the elec- trode at high C-rate operations. Before a long-term rest, the composition distributions in the particles may not be in the equilibrium state. Thus, we here examine the EIS behavior of E II with Li distributions taken immediately after lithiation processes. Figures 5.16(a) and (b) show the simulated Li fraction distributions under 0.5C and 6C lithiation at average XLi = 0.54 (from initially XLi = 0.02). As expected, relatively uniform XLi distribution under a low-rate lithiation is seen in Fig. 5.16(a). The graphite particles near the separator exhibit a relatively higher concentric XLi gradient (with surface XLi ∼0.8, as indicated by the dark yellow color, and Stage-3 cores, as indicated by the light blue color) than in the 129 particles away from the separator. For example, see the particle in the magenta dashed circle. The particles near the current collector have a relatively small XLi gradient (with surface XLi ∼0.54 in Stage 2, as indicated by the light green color, and Stage-3 cores, as indicated by the light blue color). See the particle in the cyan dashed oval in Fig. 5.16(a). Under a high rate, a prominent core-shell phase morphology is observed in the graphite par- Figure 5.16 XLi distributed in the graphite particles immediately after (a) 0.5C and (b) 6C lithiation. Response currents (c) and (d) are for configurations (a) and (b), respectively, under a cell voltage of 0.109 V. (e) Response currents for (b) under a cell voltage of 0.0 V. ticles, in which the particle surface XLi is close to ∼1 (Stage 1 phase), as indicated by the bright yellow color. The particle cores are still in a low XLi (∼0.05, Stage 1’ phase). See the magenta arrow in Fig. 5.16(b). Since the average XLi is 0.54, we first used the cell OCV (0.109 V) corresponding to XLi = 0.54 as the equilibrium cell voltage for the oscillating voltage loading on the phase 130 configuration obtained after 0.5C lithiation. All other conditions remained the same as in the baseline case. The black markers in Fig. 5.16(c) show the response current in the 0.5C- lithiation case. The amplitude of the oscillation slightly decreased over time. Besides the oscillating response, the average current approached the horizontal axis during the simula- tion, as indicated by the gray dashed line. The oscillating response current curve is difficult to be fitted with a single sinusoidal function. Thus, we cannot obtain an EIS curve of this setup. The response current in the 6C-lithiation case is provided in Fig. 5.16(d). It can be observed that the amplitude of the response increased overtime but the average of the current slightly moved to the more negative direction (see the gray dashed line). These re- sults demonstrate a ‘drifting’ during the EIS processes. During the simulations, the systems did not reach their equilibrium states yet. In the 0.5C case, the inward Li transport, which occurred to equilibrate the Li distribution, decreased the surface XLi (∼0.8) towards the equilibrium Stage-3 value (∼0.54). Thus, i0 value on those particle surfaces decreased and the amplitude of the response current decreased accordingly. In the 6C case, the inward transport decreased the surface XLi from ∼1 towards ∼0.97 (equilibrium Stage 1 value). In this case, the i0 value on those particle surfaces increased, according to the i0 function in Fig. 5.1(d). As a result, the amplitude magnified. The overall variation of the current is also related to the slow Li transport across particles to equilibrate the whole electrode. For instance, the 0.5C case was under a slowly decelerated delithiation of the electrode at the designated cell voltage. The average current was still negative but approaching towards zero. The 6C case, with a slowly increasing negative average current, was under slow delithiation over time. (Note that the equilibrium XLi mentioned above depends on the phase fraction of the system and is not a constant value throughout all cases.) Moreover, we expect that the systems in the two simulations examined above will eventually evolve to a single-phase Stage-2 case because the average XLi corresponds to the value of Stage-2 phase. However, the drifting process will likely take a very long time. In the two cases above, the equilibrium cell voltages were away from those corresponding 131 to the particle surface XLi . Based on the experience in the core-shell phase morphology EIS simulations in the previous section, we conducted another test for the 6C-lithiation case. The equilibrium cell voltage was set to be 0.0 V, which was selected based on the OCV corresponding to the particle surface XLi . However, note that this cell voltage value was roughly estimated because XLi was not uniform through all particle surfaces. Figure 5.16(e) shows the response current, which again exhibits drifting behavior. The overall current decreased toward zero and the amplitude gradually decreased toward a constant value. The positive response current shows that the particles were under lithiation. The lithiation originated from the fact that Li migrated inward to equilibrate the Li distribution at the designated cell voltage. If the simulation were to continuing over a much longer time, the system could eventually reach its equilibrium state, such that a sinusoidal response current symmetric around the horizontal axis could be observed. In that case, the EIS curve would be able to be extracted from the simulations. This set of simulations demonstrate that EIS processes can only be conducted on an equilibrium system. If an electrochemical system is still evolving towards its equilibrium, the material properties can vary during the process, which will lead to a varying amplitude of the response current as observed in the simulations. In such cases, the obtained response current cannot be fitted to a simple sinusoidal function to extract the amplitude and phase angle. Perhaps nonlinear EIS studies would be necessary to analyze those scenarios, which is beyond the scope of this paper. 5.3.5 Comparison between Phase Separation and Solid Solution Models As mentioned earlier, Fick’s diffusion is sometimes utilized to simulate Li transport inside graphite particles. In that case, phase separation will not occur. Only a Li solid solution will be observed. Here, we use our simulations to examine whether such incorrect treatments will impact the EIS behavior. In this series of electrochemical simulations, we replaced the governing equation of Li transport in the electrode particles with Fick’s diffusion equation 132 as: ∂Xp 1  |∇ψp | rxn = ∇ · ψp Dp ∇Xp + , (5.7) ∂t ψp ψp ρ where the Li diffusivity Dp is shown as the gray solid-dashed curve in Fig 5.2(b). In the first test case, we examined the EIS behavior of a stable single-phase electrode. The value of Xp was set to be uniformly 0.54 (Stage 2) as in the baseline case. All other material properties and setups remained the same as in the baseline case. The red curve in Fig. 5.17(a) shows the simulated EIS curve using Fick’s diffusion model. The EIS curve of the between Figure 5.17 (a) Simulated EIS curves of graphite electrode in single-phase Stage 2 using Fick’s diffusion model and Cahn-Hilliard model. (b) Simulated EIS curves of two-phase coexisting graphite using Fick’s diffusion and Cahn-Hilliard models, as well as the EIS curve for an artificial single phase using Fick’s law. baseline case is plotted as the blue curve on the same figure for comparison. The semicircle region of the two curves almost completely overlap. This is because the same i0 function was used in the two simulations (multi-phase model and solid-solution model). The difference in modeling Li transport barely affects the Rct . On the Warburg part, a slight difference can be observed between the red and blue curves. We attribute that small deviation to the error of parameterizing Mp from Dp . Generally, the difference in the resulting EIS curves obtained using the Cahn-Hilliard model or Fickian diffusion model is very small when the 133 initial configuration is a single-phase electrode. In the second test, we used the configuration in Fig. 5.12(a), in which stable Stage-2 (XLi = 0.54) and Stage-1 (XLi = 0.98) phases coexisted. The simulated EIS curve obtained using Fick’s law on this configuration is plotted as the red curve in Fig. 5.17(b). The EIS curve obtained using the Cahn-Hilliard equation is provided as the black curve on the same figure for comparison. Similar to the single-phase test case, the semicircle region of two curves closely overlap. Both the two curves exhibit an inductive loop near the minimum frequency region. Again, a slight difference is observed in the Warburg region due to the same issue as in the single-phase test case. In the third test case, an artificial scenario was set with uniform XLi = 0.75 throughout the entire electrode. This configuration will not exist in reality because phase separation (spinodal decomposition) will occur in graphite at this composition. A Li solid-solution within the miscibility gaps can only exist if Fick’s law is employed to mistakenly model Li transport during charge/discharge cycles. The simulated EIS curve is plotted as the blue curve in Fig. 5.17(b), which exhibits a typical semicircle with a Warburg part as in those stable single-phase cases discussed earlier. No inductive loop appears on the blue curve. The semicircle radius is significantly smaller than that of the two-phase case (either red or black curves) because the i0 at XLi = 0.75 is greater than those in Stage-2 and Stage-1 phases. The Warburg impedance is close to that of the red curve (two-phase but using Fick’s law) because value of Dp at XLi = 0.75 is similar to those in Stages 2 and 1. Even though the scenario examined here is unphysical, our electrochemical simulations still properly reflected the impact of input material properties and phase morphologies on the resulting EIS curve. Based on the results of these simulations, we find that the EIS curves obtained using either the Cahn-Hilliard model or Fickian model on an equilibrium, stable configuration is very similar. Thus, since experimental EIS measurements are conducted at equilibrium con- ditions, it is difficult to distinguish whether Li transport proceeds via a two-phase coexisting or solid-solution process from the measured EIS data. As mentioned in the core-shell two- 134 phase case, perhaps EIS can detect a material’s phase transformation only if the measured properties (e.g., semicircle diameter) exhibit step-wise variation versus the electrode’s DOD. 135 5.4 Conclusion In this Chapter, we used the Cahn-Hilliard equation to model the Li transport in graphite particles such that two- (or multiple) phase coexistence in experimentally reconstructed graphite electrodes can be considered in physics-based simulations. We examined the EIS behavior of stable single phase cases on different electrode microstructures and found that EIS measurements could mislead to a counter-intuitive scenario that an electrode exhibiting a high resistance but posses an enhanced high-rate performance. The effective diffusivity, calculated based on the Warburg impedance, could be underestimated if the particles compo- sition enters the miscibility gap under a large loading amplitude. In two-phase coexistence cases, low-frequency inductive loops appear on the Nyquist plots if phase boundaries are present on particle surfaces. Furthermore, EIS measurements cannot distinguish the differ- ence between a core-shell phase coexistence and a single uniform phase, as EIS only probes the particle surface conditions. The simulations also demonstrate that drifting occurs if the system is away from equilibrium. Lastly, although Fick’s diffusion is widely mistakenly employed in simulating Li transport in phase-separating graphite particles, the EIS curves obtained using the Cahn-Hilliard equation and Fick’s diffusion equation is very similar if the system is in equilibrium. Determining if an electrode particle is phase-separating or solid- solution material should be based on whether the EIS semicircle diameter (or other measured properties) varies step-wisely versus the DOD. We expect this presented simulation tool to be widely used to investigate complex EIS behavior of many other electrochemical systems. 136 CHAPTER 6 FULL CELL SIMULATION 137 6.1 Introduction In the previous chapters, we have demonstrated that the presented simulation framework worked properly for half cell configurations with NMC in Chapter 4 and with graphite in Chapter 5. For both cases, lithium metal was used as the counter electrode. In those simulations, we solved the governing equations of electrochemical processes in the electrolyte, in the electrode particles, and on the particle-electrolyte interfaces. A full cell simulation will require to solve two sets of equations for the mass transport and current continuity conditions: one set for the cathode and the other set of the anode. Furthermore, there will be two sets of Butler-Volmer equations: one on the cathode-electrolyte interface and the other on the anode-electrolyte interface. While there is only one governing equation for salt concentration evolution in the electrolyte, This equation is subject to Neumann boundary conditions (specifying reaction fluxes) on both sides. Fig. 6.1 shows the full cell Figure 6.1 Full cell model used in simulation, the gray part represents graphite electrode and the yellow part represents cathode electrode; the rest pores and interparticle space are filled with electrolyte. 138 configuration used in the simulations, which consists of the NMC cathode and graphite anode used in Chapter 4 and Chapter 5, respectively. The yellow color indicates spherical NMC particles and the microstructure is generated using discrete element method. The gray color indicates the graphite particles reconstructed using X-ray computed nano-tomography [126, 127]. The transparent regions are interparticle space and filled with liquid electrolyte. The empty space between the cathode and anode serves as the separator to avoid direct contact between the two electrodes. In the conventional sharp-interface methods, the microstructures must be discretized to mesh conformal with the complex geometries. Electrochemical simulations of a symmetric cell [68] consisting of 3D electrode microstructures have been performed to investigate EIS processes and compared the resulting EIS with electric circuit model and macro-homogenized physical models (e.g., PET model). Full-cell 3D microstructure simulations were also per- formed to investigate the effect of structuring techniques (introducing tunnels in the elec- trodes) on the electrochemical performance and thermal behavior of lithium-ion batteries [66]. These simulations all involved tedious mesh generation processes. In this chapter, we apply our diffuse-interface method to simulate the EIS processes in a full cell configuration, such that no body-conforming mesh is required. The simulations reveal detailed electrochem- ical dynamics occurring in the cathode and anode, which allow us to connect the macroscopic EIS behavior to the intrinsic material parameters of the two electrode materials and to the electrode microstructures. As demonstrated in this work, the ease of implementing full-cell complex electrode microstructure simulations opens a door for investigating electrochemical processes using simulations on a large number of different microstructures. 139 6.2 Method 6.2.1 Material Properties and Microstructures In this chapter, NMC-111 and graphite are used as the cathode material and the anode material. Synthetic sphere agglomerate cathode and experimentally reconstructed graphite electrode II, which have been employed in half-cell simulations in Chapter 4 and 5, respec- tively, are used in the full-cell simulations. The material properties remain the same as in the previous chapters. The dimensions of computational box are 320×150×150 along the x, y, and z directions, which correspond to 104×48.75×48.75 µm in physical unit. The thick- nesses of the cathode and anode are similar to typical dimensions of real battery electrodes, as used in full-cell testing experiments. 6.2.2 Full Cell Solver Workflow As discussed in Chapter 4, NMC is treated as a solid solution of Li. Thus, Fick’s diffusion equation is used to describe the lithium concentration evolution in the NMC particles.  |∇ψc |ce rxn c ∂Xc 1 = ∇ · ψc Dc ∇Xc + . (6.1) ∂t ψc ψc ρc This equation is the same as Eq. (4.9). Note that the subscript c indicates cathode particles. The subscript ce indicates cathode particle-electrolyte interface. The (de)intercalation rate, c rxn , on cathode particle surface is described by the Bulter-Volmer equation:     c −αz+ F c (1 − α) z+ F c rxn = kfc Ce exp c [ϕ]e − kb Cc exp [ϕ]e , (6.2) RT RT where [ϕ]ae is the electrostatic potential drop across cathode particle-electrolyte interface. This equation also provides the flux boundary condition on the cathode particle-electrolyte interface for the salt concentration evolution in the electrolyte: c a ∂Ce 1 |∇ψe |ce rxn t− |∇ψe |ae rxn t− = ∇ · (ψe De ∇Ce ) + ++ , (6.3) ∂t ψe ψe ν+ ψe ν+ 140 where the second and third terms are for the reaction rate on the cathode particle surfaces and the anode particle surfaces, respectively. As discussed in Chapter 5, phase transformation occurs in graphite during (de)lithiation. The Cahn-Hilliard phase field equation is utilized to model the phase transformation pro- cesses:  |∇ψa |ae rxn a ∂Xa 1 = ∇ · ψa Ma ∇ µa (X) − ε∇2 Xa + , (6.4) ∂t ψa ψa ρa where the subscript a indicates anode particles, and the subscript ae indicates the anode particle-electrolyte interface. The reaction rate on the anode particle-electrolyte interface is also governed by the Butler-Volmer equation as     a −αz+ F a a (1 − α) z+ F a rxn = kfa Ce exp [ϕ]e − kb Ca exp [ϕ]e . (6.5) RT RT Note that there are three regions inside the entire computational domain: cathode, elec- trolyte, and anode. These three regions are defined by domain parameters: ψc , ψe , and ψa . Equations (6.1), (6.3), and (6.4) update the Li fractions and salt concentration in the respective regions. Forward Euler methods are employed for Eqs. (6.1) and (6.4) between time steps. Backward Euler method is used for Eq. (6.3). Within each time step, three Poisson’s equations are solved for the electrostatic potentials in the cathode, electrolyte, and anode regions: c ∇ · (ψc κc ∇ϕc ) − |∇ψc |ce z+ F rxn = 0, (6.6) c rxn ra ∇ · [ψe (z+ m+ − z− m− ) F Ce ∇ϕe ] + |∇ψe |ce +|∇ψe |ae xn ν+ ν+ (6.7) =∇ · [ψe (D− − D+ ) ∇Ce ] , a ∇ · (ψa κa ∇ϕa ) − |∇ψa |z+ F rxn = 0. (6.8) Equation (6.6) is subject to a Dirichlet boundary condition, ϕc |cc , on the computation box c boundary (x = 104 µm) and is subject to an internal boundary condition, rxn . The subscript ‘cc’ indicate cathode current collector. The second and third terms on the left hand side of Eq. (6.7) are the reaction rates on the cathode and anode particle surfaces, respectively. 141 Equation (6.8) is subject to a Dirichlet boundary condition, ϕa |ac , on the computation box c boundary (x = 0 µm) and is subject to an internal boundary condition, rxn . The subscript ‘ac’ indicates anode current collector. For the three electrostatic potential fields, Eqs. (6.6), (6.2), (6.7), (6.5), and (6.8) are solved in an iterative scheme. An oscillatory voltage loading is first imposed on the cathode current collector by setting ϕc |c = ϕc0 + V0 sin(ωt), where ϕc0 is the equilibrium potential, V0 is the amplitude, and ω is the angular frequency. Then, we adjust the value of ϕa |ac , and relax Eqs. (6.6), (6.2), (6.7), (6.5), and (6.8) until all these five equations reach their numerical equilibrium. Next, we calculate the total current on the cathode side which c c comprises of reaction and capacitance currents (Itot = Irxn + Icc , as mentioned in Chapter 4), and that on the anode side (Itot a a = Irxn + Ica ). These two total currents have opposite signs. a Positive/negative means insertion/extraction from the particle. If the magnitude of Itot is c different from that of Itot , we adjust ϕa |ac and relax the five equations again. This process c a is repeated until the difference between the magnitudes of Itot and Itot is less than some c threshold value. In this work, the threshold is 0.25% of the magnitude of Itot , below which c a the equilibrium ϕc , ϕe , ϕa , rxn , and rxn are accepted as the solutions and the cell voltage across the cathode and anode is obtained by ϕCV = ϕc |cc − ϕa |ac . (6.9) c a In the next time step, Xc , Ce , and Xa are updated with rxn and rxn . See the flow chart in Fig. 6.2. 6.2.3 EIS Curve Extraction The EIS curve is comprised of many data points, each of which is the impedance value at a loading frequency. The impedance value is calculated as Z = V /I. As mentioned in the previous section, the total cell voltage (ϕCV ) is different from the imposed sinusoidal voltage (ϕc |cc = V0 sin(ωt)) on the cathode current collector because ϕa |ac is adjusted to ensure the 142 Figure 6.2 Flowchart of the simulation scheme for solving the coupled equations in a full cell configuration. conservation of currents throughout the whole cell. For a close-to-fully-discharged cell, Li fraction on the cathode is set to be 0.92 and that on the anode is 0.04. Based on the OCV curve of NMC, the equilibrium potential for the cathode is ϕc0 = 3.350 V at Xc = 0.92. The equilibrium potential for the anode is ϕa0 = 0.246 V at Xa = 0.04. A slight sinusoidal potential perturbation is imported on the cathode current collector, with V0 = 10 mV and ordinary frequency f = 32 Hz, which is shown on Fig. 6.3(a). The electrostatic potential on the anode current collector is adjusted to maintain the conservation of current. The result ϕa |ac versus time is shown in Fig. 6.3(b), which also exhibits a sinusoidal function oscillating around the equilibrium potential ϕa0 = 0.246 V with an amplitude about 2.5 mV. Note that, 143 Figure 6.3 (a) Sinusoidal voltage loading imposed on the cathode. (b) The resulting electrostatic potential on the anode current collector in order to maintain a current balance. (c) The cell voltage across the cathode and the anode current collectors. due to the numerical scheme of adjusting ϕa |ac , the values in the first half of the period is very low. The sinusoidal form is well developed after two periods. Interestingly, there is a small phase shift between ϕc |cc in Fig. 6.3(a) and ϕa |ac in Fig. 6.3(b). Thus, the total voltage cannot be obtained by simply adding the two amplitudes. The cell voltage, Eq. (6.9), is plotted as the black circles in Fig. 6.3(c). The fitted sinusoidal curve is shown as the yellow curve in the same figure. The amplitude of the cell voltage is approximately 11.5 mV. Note that the response voltage on the anode side can change with frequency. Based on our simulations, at a higher frequency, the anode response voltage amplitude becomes larger. Furthermore, 144 because graphite is a conductive material, NMC at XLi = 0.92 is reasonably conductive, and ionic diffusivity in the electrolyte is high, ϕa , ϕc , and ϕe all have a small variation throughout their respective domains. Thus, ϕa |ac and ϕc |cc can be viewed as approximately the electrostatic potential drop across the anode particle-electrolyte interface and that across the cathode particle-electrolyte interface, respectively, to result in the same magnitude of total anode current and total cathode current, but with opposite signs. The simulation reveals the difference in the electrostatic potential drops on the cathode particle surfaces and on the anode particle surfaces, and reveals the phase shift between those sinusoidal voltage responses. These phenomena taking place inside a battery cell are difficult to be directly observed in experiments. Figure 6.4 Total current (a) across the cathode and electrolyte interface (b) across the anode and electrolyte interface As mentioned in the previous section, conservation of total current is maintained at each time step. The magnitudes of cathode current and anode current are the same, but the two currents have opposite signs. Figures 6.4(a) and (b) show the total current on the cathode and anode sides, respectively. These two curves can be fitted to almost the same sinusoidal function, but with different signs. With the obtained cell voltage and total current, the impedance is then calculated by Z(f ) = ϕCV /Itot . In the next section, the 145 full-cell simulations are performed over a frequency range from 512 Hz to 0.05 Hz to extract the impedance spectroscopy. The total current on each of the electrodes consists of reaction current and capacitance current. Figures 6.5(a) and (b) show the reaction current and capacitance currents across Figure 6.5 (a) Capacitive current and (b) reaction current across the cathode particle-electrolyte interfaces. (c) Capacitive current and (d) reaction current across the anode particle-electrolyte interfaces. the cathode particle-electrolyte interfaces. The capacitance current is roughly one order of magnitude large than the reaction current at this case (f = 32 Hz and Xc = 0.92). On the anode side, the reaction current is approximately 60 times larger than the capacitance current. See Fig. 6.5(c) and (d). The anode reaction current is about one order of magnitude 146 larger than the cathode reaction current because i0 of graphite at Xa = 0.04 is much larger than i0 of NMC at Xc = 0.92. Most of the total cathode current is matched by the anode reaction current. Thus, the anode capacitance current is much smaller, which compensates the small difference between the total cathode current and anode reaction current. Note that there are noises on the anode capacitance current in Fig. 6.5(c) due to the adjustment processes of ϕa |ac . In some cases, those noises can be large and will be removed before the curve fitting procedure. 147 parameter description NMC graphite A surface area [cm2 ] 1.210×10−4 3.086×10−4 V volume [cm3 ] 5.425×10−8 7.584×10−8 ρLi Li site density [mol/cm3 ] 0.0501 0.312 Csp specific capacitance [F/cm2 ] 1.40×10−4 1.40×10−4 Table 6.1 Electrode property parameters used in simulation 6.3 Result and Discussion Three full-cell simulations are performed to examine the EIS behavior at different conditions. In the first study, all material properties are kept the same as in Chapters 4 and 5. The main focus is to figure out how the state of charge affects the resulting EIS curves. We first simulate the EIS process of a close-to-fully discharged cell. Next, the simulation is performed on the same cell but at a close-to-fully charged state. Hereafter, these two cases are referred to as fully discharged and full charged cases for convenience. The two EIS curves significantly differ. However, on those curves, the two semicircles (each for one of the electrodes) highly overlap. Thus, in the last case, the simulation is conducted with with different values of specific capacitances on the two electrodes such that two distinct semicircles can be observed on the simulated EIS curve. 6.3.1 Fully Discharged Configuration The first full-cell simulation is performed to study the EIS behavior of a fully discharged cell. In this case, most of the lithium in the graphite anode is extracted and is inserted into the NMC cathode. Therefore, the lithium fraction is set to be 0.04 in the graphite anode and to be 0.92 in the NMC cathode. Here, we set the upper utilization limit of NMC particles to be XLi = 0.95, because above which NMC becomes an electrically insulating ceramic metal oxide. This study aims to investigate the effect of different i0 values of the anode and cathode on the full-cell EIS behavior. Thus, the same specific capacitance is used on both the anode and cathode sides. As mentioned in Chapter 5, graphite has four stable single 148 phases. At Xa = 0.04, the lithium-graphite system will stay the Stage-1’ phase. Based on the Li-in-graphite OCV curve (Fig. 5.1(a)), the equilibrium potential for the anode is 0.246 V. Similarly, according to the Li-in-NMC OCV curve (Fig. 4.4(d)), the equilibrium potential for the cathode is 3.350 V. Using the method described in the previous section, we calculate the sinusoidal cell voltage loading and the sinusoidal response total current at each loading frequency. Impedance values for 15 frequencies, ranging from 512 to 0.05 Hz, are calculated from the simulation results. For every frequency point, Li concentration and electrostatic potential are solved in the full cell scale. An example of fully discharged cell is shown on Fig. 6.6 at f = 1 Hz, t = 1.5 s. Figure 6.6 Example images of simulation results in the fully discharged cell at f = 1 Hz and t = 1.5 s: (a) Li fraction in the graphite anode particles, (b) Li fraction in the NMC cathode particles, and (c) salt concentration in the electrolyte. 149 Figure 6.7 (a) Nyquist plot of EIS curve for a fully discharged cell. (b) Magnified view of the EIS curve in the high frequency region. Figure 6.7(a) shows the EIS curve of a fully discharged cell on the Nyquist plot. The shape of the EIS curve appears as a large semicircle, though only part of the circle is obtained from the simulations. The Warburg tail has not been reached in this figure because the frequencies sampled are not low enough. Here, the explicit Euler method for time integration is used in updating the Li fractions in NMC (Eq. (6.1)) and graphite (Eq. (6.4)). Because of Li diffusivity in graphite (at XLi = 0.04) is about 3–4 times larger than that in NMC (at XLi = 0.92), the stable time step size determined by Li diffusivity in graphite is too small for updating of Xc in NMC. Simulations for frequencies lower than those shown above will require significantly longer times. Thus, we stop the simulations at 0.05 Hz. While the EIS curve in Fig. 6.7(a) seems to be a single semicircle, it in fact is consisted of two semicircles: one very large and the other very small. Figure 6.7(b) shows the EIS curve in the high frequency region, from which part of the small semicircle can be observed near the origin. Excluding the mass transport impedance, the analytical formula of impedance for an RC-RC circuit is Ra Rc Z = Za + Zc = + , (6.10) 1 + jf Ca Ra 1 + jf Cc Rc where j is the imaginary unit, Ri is the resistance, Ci is the capacitance, and the subscripts 150 Figure 6.8 Illustration of the RC-RC circuit. a and c indicate anode and cathode, respectively. An illustration of RC-RC circuit is shown on Fig. 6.9. Assuming Ca < Cc , the center of anode semicircle is Ra /2 and the center of cathode semicircle is Ra + Rc /2. Note that since our electrochemical simulations did not reach the Warburg tail yet, we do not include diffusion-impedance components in the circuit model above. The cathode semicircle will be on the right to the anode one because the low frequency region is dominated by the high capacitance component. Ca and Cc determine the impedance point on the EIS curve. If Ca and Cc have the same value, the EIS curve forms a semicircle with a radius of nearly (Ra + Rc ). Here, since the specific capacitance used for the two electrodes is the same, the total capacitance of the two electrodes are proportional to the respective particle surface areas. The surface area of the graphite anode is 3.086×104 µm2 , and the surface area of the NMC cathode is 1.210×104 µm2 . The values are provided in Table 6.1. Thus, the total capacitance of the NMC cathode (4.3204 × 10−2 µF) is roughly only 2.5 times that of the anode total capacitance (1.694 × 10−2 µF). The exchange current density i0 for graphite at XLi = 0.04 is 1.050×10−3 A/cm2 , and that for NMC at XLi = 0.92 is 2.474×10−6 A/cm2 . Using Eq. (5.6), the total charge transfer resistance of the graphite anode is 7.9788 × 104 Ω, and that for the NMC cathode is 8.6312 × 107 Ω. The calculated 151 charge-transfer resistances are provided in Table 6.2. For these values, the semicircles of the anode, the cathode, and the full cell predicted using the circuit model are plotted in Fig. 6.9. The radius of the cathode semicircle is roughly three orders of magnitude larger Figure 6.9 The Nyquist plot of EIS curves calculated from RC-RC circuit model, Eq. (6.10), using the resistance and capacitance data of (a) fully discharged cell and (b) magnified view of (a) in the high frequency range. (c) Fully charge cell and (d) fully charged cell with modified material parameters. than that of the anode semicircle. In this case, because the total capacitances of the two electrodes are similar, the small anode semicircle is mostly merged into the large cathode semicircle. Only a small fraction of the anode semicircle can still be traced on the EIS curve. The EIS curve extracted from microstructure simulations is very similar to that obtained from the RC-RC circuit model. The total resistance of the full cell from the electrochemical simulations is 6.282×107 Ω according to the radius of the fitted circle, which is smaller than 152 the value predicted from the circuit model. The discrepancy may be due to numerical error R of calculating total active surface areas. The surface area calculated using |∇ψ|dΩ in the electrochemical simulation can be different from the value calculated by summing all triangular patches of the isosurface. This error can be reduced by using finer mesh. This set of simulations demonstrate that, in a fully discharged cell, the resistance on the EIS curve is dominated by the cathode property because the i0 of the NMC cathode is very small (which leads to a large total charge-transfer resistance). Moreover, an interesting phenomenon is observed in the simulations. As demonstrated in Section 6.2.3, even though the same specific capacitance is used on both the anode and cathode, the resulting capacitance currents on the two electrodes are significantly different (here, Ica ≪ Icc ); for example, as shown in Figs. 6.5(a) and (c). Because graphite i0 at XLi = 0.04 is much larger than NMC i0 at XLi = 0.92, most of the total cathode current is matched by the anode reaction current. This simulation illustrates an interesting distribution between reaction and capacitance currents in an electrode to satisfy the conservation of total current. The distribution is determined by the electrode’s i0 value if the capacitances are similar in the anode and cathode. 6.3.2 Fully Charged Configuration In the second study, the cell is charged to a fully charged state. In this case, Li is extracted from the NMC cathode and is inserted into the graphite anode. The average Li fraction in the NMC cathode decreases from XLi = 0.92 to 0.25. Here we set the lower utilization limit of NMC to be XLi = 0.25, because below which exfoliation of the NMC layered oxide occurs and the cathode decomposes. For NMC at XLi = 0.25, the equilibrium voltage is 4.169 V (see Fig. 4.4(d)). The average Li fraction in the graphite anode increases from XLi = 0.04 to 0.81 in the fully charged cell according to their Li site densities. As described in Chapter 5, a core-shell phase morphology exhibits in the graphite particles at this Li fraction. The surfaces are in Stage-1 phase with Xa = 0.97 and the particle cores are in Stage-2 phase with 153 Xa = 0.56. As concluded in Chapter 5, the EIS measurement will probe only the particle surface properties. Thus, treating all graphite particles to be uniformly in Stage-1 phase will offer the same simulation result as that obtained from a core-shell phase morphology. For graphite at XLi = 0.97, the equilibrium voltage is 0.0854 V (see Fig. 5.1(a)). In this fully charged case, the graphite i0 is 7.2 × 10−4 A/cm2 and the NMC i0 is 1.448 × 10−4 A/cm2 , which lead to the theoretical total charge-transfer resistances for the graphite anode and for the NMC cathode to be 1.1636 × 105 Ω and 1.475 × 106 Ω, respectively. Note that the cathode resistance is only roughly one order of magnitude larger than that of the anode. This ratio is much smaller compared with the fully discharged case (roughly three orders of magnitude). Because the same Csp is used in the simulations, the theoretical total capacitance is still 1.694 × 10−2 µF and 4.3204 × 10−2 µF for the graphite anode and the NMC cathode, respectively, as in the previous section. All other simulation conditions are the same as in the previous section. The simulations are performed for 15 loading frequencies ranging from 512 to 0.05 Hz. Figure 6.10 (a) Nyquist plot of simulated EIS curve for a fully charged cell. (b) EIS curve from simulated data together with two fitted semi-circles using different segments of the curve. Figure 6.10(a) shows the EIS curve extracted from the simulations of a fully charged cell. The impedance curve seems to exhibit a single semicircle with a Warburg tail. Here, because 154 the Li diffusivity in graphite (at XLi = 0.97) is similar to that in NMC (at XLi = 0.25), the stable time step size determined by the Li diffusivity in graphite is also suitable for updating Li fraction in NMC. Thus, the result can extend to the Warburg region in a reasonable simulation time. It can be noted that the shape of the semicircle is distorted: the region corresponding to the intermediate-to-high frequencies is a little compressed. The impedance data points in the semicircle region cannot be perfectly fitted with one single circle. Figure 6.10(b) shows two circles fitted from the impedance data. The red circle is fitted using the points at frequencies from 4 to 0.5 Hz. Its diameter is 1.276 × 106 Ω. Since the cathode has a larger capacitance, the red circle obtained from the low frequency points is more affected by the cathode resistance. On the other hand, the green circle, fitted using the points at higher frequencies from 512 to 16 Hz, is more affected by the anode resistance. The diameter of the green circle is 1.187 × 106 Ω, which is slightly smaller than that of the red circle. Similar to the fully discharged case in the previous section, because the two electrodes have similar capacitances, the semicircles of the anode and the cathode merge into one EIS arch. Therefore, resistance estimated from the circle fitting may not be quantitatively accurate although the values are reasonably close to that obtained from the RC-RC circuit model (Ra + Rc = 1.5914 × 106 Ω). Again, the deviation can be originating from the numerical error of calculating the surface area in the simulations. Moreover, the EIS curve extracted from the physics-based electrochemical simulations is similar to that obtained from the RC-RC circuit model, see the EIS curve in Fig. 6.9(b). The presented physics-based microstructure electrochemical simulations bridge across the intrinsic material properties, electrode microstructures, and the macroscopic EIS measured quantities. 6.3.3 Fully Charged Cell with Modified Material Parameters In the two previous cases, the semicircles of the anode and the cathode deeply merged. In many experimental measurements, two separate semicircles can be observed. As discussed earlier, the capacitances of the two electrodes need to be sufficiently different such that the 155 two semicircle can be distinguished. In Ref. [66], the value of specific capacitance of double layer on the NMC electrode is roughly a hundred times that on the graphite electrode. As discussed in Chapters 1 and 3, the charge separation in the electrolyte is determined by the dielectric constant of the liquid. Since the two electrodes are immersed in the same electrolyte, the charge separation region should be the same. However, the thicknesses (λs in Eq. (3.4)) of the Helmholtz layer can be significantly different, depending on the electrode surface structures. Thus, specific capacitance can be different on the anode particle surfaces and the cathode particle surfaces. Here, we increase Csp on NMC surfaces by a factor of five to 7 × 10−4 F/cm2 and decrease that on graphite surfaces by ten folds 1.4 × 10−5 F/cm2 . Note that because there is a large range of Csp values reported in the literature data, these values are selected simply to highlight the effect of significantly different capacitances on the two electrodes. Furthermore, in the previous case, the charge-transfer resistance on NMC surface is ten times that on the graphite surface. To result in two semicircles with comparable radii, the exchange current density of NMC is magnified by a factor of three. The relevant quantities are also provided in Table 6.2. Simulations are performed for 16 frequencies ranging from 1024 to 0.1 Hz. Figure 6.11 (a) Nyquist plot of simulated EIS curve for a fully charged cell with modified material parameters. (b) EIS curve from simulated data together with two fitted semi-circles using different segments of the curve. 156 Figure 6.11(a) shows the EIS curve extracted from the simulations. Two distinct semi- circles are observed on the curve: the small one in the high frequency region and the large one in the intermediate frequency region. The Warburg tail appears in the low frequency region. These two arches are fitted with two circles as shown in Fig. 6.11(b). The circles well overlap the data points. The red circle is fitted using points corresponding to frequen- cies from 1024 to 128 Hz. In this regime, the capacitance current stemming from double layer formation is comparable to the reaction current on the anode particle surfaces. The diameter of the red circle is 1.297 × 105 Ω. The green circle is fitted using the data points corresponding to frequencies from 16 to 0.5 Hz, which is the regime where the capacitance and reaction currents on the NMC cathode are comparable. The diamter of the green circle is 4.678 × 105 Ω. These charge-transfer resistances are reasonably close to those calculated using Eq. (5.6): 1.1636 × 105 Ω for the anode and 4.9181 × 105 Ω for the cathode. Since the specific capacitances used in this set of simulations has a 50 folds difference between the cathode and the anode, the semicircles are separate enough from each other. Furthermore, the EIS curve extracted from the simulations is similar to curve predicted using RC-RC circuit model (shown in Fig. 6.9(c)). This demonstrates that the presented microstructure level simulation properly reflects the input material parameters. As ready mentioned, this physics-based electrochemical simulation tool can connect material properties, microstruc- tures, and measured macroscale quantities. As material properties usually contain significant uncertainties, this tool can be utilized to calibrate material parameters. 157 Case Fully Discharged Fully Charged I Fully Charged II Li fraction of anode [%] 0.04 0.97 0.97 Li fraction of cathode [%] 0.92 0.25 0.25 i0 of anode [A/cm2 ] 1.050×10−3 7.200×10−4 7.200×10−4 i0 of cathode [A/cm2 ] 2.474×10−6 1.448×10−4 4.344×10−4 ae Calculated Rct [Ω] 7.979×104 1.164×105 1.164×105 ce Calculated Rct [Ω] 8.631×107 1.475×106 4.918×105 Ca [µF] 1.694 × 10−2 1.694 × 10−2 1.694 × 10−3 Cc [µF] 4.3204 × 10−2 4.3204 × 10−2 2.1602 × 10−1 Table 6.2 Parameters for the different full cell cases. 6.4 Conclusions In this chapter, a NMC cathode microstructure and a graphite anode microstructure are combined to form a full-cell configuration for simulating the EIS process. The EIS curves of the full cell at different states of charge are studied. At the fully discharged state, because the charge transfer resistance of NMC is much higher than that of graphite, the NMC semicircle dominates the entire EIS curve of the full cell. At the fully charged state, the Rct of two electrodes are comparable. However, due to the fact that the Csp used for the two electrodes is the same in the simulations, the semicircles for the anode and cathode deeply merge. With modified values of Rct and Csp , the EIS curve resembles those experimentally observed, show- ing two distinct semicircles and a Warburg tail. The simulations also reveal an interesting voltage distribution between the cathode and anode, as well as an distribution of reaction current and capacitance current in each of the electrodes. These microstructure phenomena is difficult to be directly observed in EIS experiments. This work demonstrates that the physics-based microstructure electrochemical simulations can provide a proper connection between intrinsic material properties, microstructures, and macroscale measured quantities. As material properties usually contain significant uncertainties, this presented simulation framework will be powerful for calibrating intrinsic material properties from measured EIS curves, with experimentally reconstructed electrode microstructures. The simulations on different cathode microstructures illustrate that, for the same material, the charge-transfer 158 resistance on the EIS measurement is also inversely proportional the active surface areas. 159 CHAPTER 7 SUMMARY AND OUTLOOK 160 7.1 Summary Battery operations involve coupled multiphysics electrochemical processes, and these pro- cesses occur in highly complex electrode microstructures. Studying the detailed electrochem- ical processes in very challenging because of the complexities of both mutliphysics phenomena and microstructures. Electrochemical impedance spectroscopy (EIS) technique is widely em- ployed to measure battery electrode properties. However, the obtained quantities, such as capacitance, resistance, and diffusional impedance, are for the whole cell. The underlying connections between the electrode materials’ intrinsic properties, electrode microstructures, and the measured macroscale values, remain poorly understood. This work aims to de- velop a simulation framework that allows us to simulate the detailed process taking place in electrode microstructures during a EIS measurement process. Such that the relationship between material properties, microstructures, and measured EIS values can be elucidated in the simulations. As double layer capacitance plays an important role in the EIS processes, the starting point is to develop a method to simulate double layer formation in complex geometries. The smoothed boundary method (SBM) is employed to reformulate the Nernst-Planck-Poisson (NPP) equations, so that the new governing equations can be solved on mesh non-conformal to the complex geometries. This method greatly accelerates the simulation implementa- tion. Adaptive mesh refinement (AMR) is used to reduce the modeling error stemming from the diffuse interface in the SBM. It also significantly decreases the computational bur- dens. One-dimensional simulations, verified against analytical solutions, demonstrate the accuracy of the presented SBM-AMR method for simulating double layer formation. Multi- dimensional simulations reveal a two-step double layer formation process. First, ions are rapidly adsorbed on (or repelled from) the particle surfaces that have short distances to another particle surfaces. Then, adsorbed ions diffuse around the particle until the double layer uniformly surrounds the whole particle. In the electrolyte examined in the simulations, charge separation (double layer formation) reaches equilibrium in the millisecond scale. 161 To simulate the EIS processes in battery electrodes, we develop a multi-physics electro- chemical model, which solves the coupled governing equations for Li diffusion and current continuity in electrode particles, ionic diffusion and current continuity in the electrolyte, and chemical reaction on the particle surface. The developed SBM-AMR approach is applied to simulate these electrochemical processes in complex electrode microstructures. The double layer capacitance calculated by solving the NPP equations is incorporated into the elec- trochemical modeling. Electrochemical simulations are performed under sinusoidal voltage loading, and impedance values are calculated from the voltage loading and current response. We start with NMC cathode, which is a Li solid solution. Thus, Fick’s diffusion equation is used to model Li transport in the electrode particles. The developed simulation tool is utilized to investigate the impact of exchange current density at different state of charge on the EIS curves. The results demonstrate that the charge-transfer resistance is inversely proportional to the exchange current density. We also study how salt concentration in the electrolyte affects the EIS curves. While the salt concentration affects exchange current den- sity, double-layer capacitance, and ionic diffusivities simultaneously, the simulation results show that the changes of exchange current density has the largest effect on the EIS behavior. In contrast to the Li-solid-solution NMC cathode particles, graphite particles exhibit second-order phase transition upon lithiation/delithiation. Thus, the Cahn-Hilliard phase- field equation is employed to model the phase transition processes in graphite particles. The simulations are performed on experimentally reconstructed graphite electrode microstruc- tures. The results show that the charge-transfer resistance on the EIS curve is inversely proportional to the product of exchange current density and active surface area when the graphite particles are in their four stable single phases. While an electrode with a low tortu- osity of the electrolyte channels can have a better high-rate performance, it likely has a larger resistance on the EIS curve, compared to those having a worse high-rate performance. The results illustrate that this counter-intuitive behavior is due to the fact that a low-tortuosity electrode usually has a smaller active surface area. In graphite particles with a core-shell 162 phase morphology, the EIS measurement only probe the particle surface properties. Our simulations reveal an interesting phenomenon that if phase boundaries intersect particle sur- faces, a low-frequency inductive loop appears on the EIS curve. While Fick’s diffusion model is commonly mistakenly used to model the Li transport during the phase transition processes in graphite, our simulations show that the EIS curves obtained from the phase-field model and from Fick’s diffusion model are indistinguishable. Therefore, EIS measurements cannot detect whether a material can undergo phase transition if this material is in equilibrium. A material undergoes phase transition can only be discerned by whether the EIS measured quantities exhibit a step-wise variation. Lastly, the simulations are performed on a full cell configuration consisting of a NMC cathode microstructure and a graphite anode microstructure. EIS curves of the full cell are extracted from those simulations. While typical experiments only access the cell voltage between the two electrodes, the presented electrochemical simulations show the individual electrostatic potential drop across the cathode particle-electrolyte interfaces and that across the anode particle-electrolyte interfaces. These two quantities impose the conservation of currents in the electrodes. The simulations also reveal a distribution between reaction cur- rent and capacitance current on each of the two electrodes, which is determined by the double-layer capacitance and exchange current density according to the state of charge of that electrode. Ideally, the anode and the cathode each has a semicircle on the EIS curve. The diameter of that semicircle is the charge-transfer resistance of that electrode. If the capacitances of the two electrodes are similar, the two semicircles deeply merge and form one large semicircle. Our simulation results demonstrate this behavior on the resulting EIS curves. If the capacitance of the two electrodes are sufficiently different, the two semicircle separate apart on the EIS curve, which is also illustrated in the simulations. In conclusion, the presented simulation framework incorporates the double-layer capacitance into electro- chemical simulations on electrode microstructures. It successfully provides the bridge to connect the intrinsic material properties, electrode microstructure, and the measured EIS 163 quantities. This open a new door using detailed multiphysics simulations in battery research. 7.2 Outlook Overall, we have successfully developed a physics-based model, which connects electrode microstructures and material properties to their electrochemical performance and resulting EIS curves. However, many material parameters for battery electrode materials, such as exchange current density, double layer capacitance, Li diffusivity, etc., are scarce in the liter- ature data. Even if the data are available, the values contain a high degree of uncertainties due to the difficulties in conducting the relevant measurements. Thus, a potential use of this simulation framework is to calibrate input material parameters based on measured EIS data. Furthermore, additional mechanisms can be incorporated on to the presented electro- chemical model. For instance, the resistance due to the formation of solid electrolyte inter- face (SEI) on the electrode particle surfaces can be included into the model. The resistance increase due to lattice structure change in the electrode particles could also be included in the model. Li loss due to side reaction between plated Li metal and the organic electrolyte can also be considered. In this case, the total active Li in the system could be estimated from the sizes of the two semicircles of the two electrodes. This framework uses the Jacobi relaxation method in solving all the static equations. There are other numerical schemes that can be employed to accelerate the simulations; for instance, using over-relaxation or conjugate gradient methods. The SBM-reformulated governing equations can be solved using finite element methods on non-body-conforming mesh systems. This may be more stable than using finite difference method as in this work. Moreover, there are five coupled equations that are solved using iterative schemes in the full cell simulation code. The iteration usually converges slowly. Particularly, the adjustment procedure in controlling the current conservation frequently leads to noises on the response current. Often, it can lead to divergence of the simulations. A better solver scheme could 164 greatly improve the code performance. Adaptive time steps can also increase the speed of simulations. This framework can be extended to study other different electrochemical systems. For instance, the simulations can be used to simulate the electrochemical processes and EIS behavior of solid-state batteries. This framework can also be applied to simulate the EIS behavior of solid oxide fuel cells. Redox flow cell could be a potential application if computa- tional fluid dynamics is incorporated to the model. Further, the charge separation model for the double layer formation can be combined with computational fluid dynamics to simulate the processes of capacitive desalination for water purification. 165 BIBLIOGRAPHY [1] Sara Taslimi Taleghani, Bernard Marcos, Karim Zaghib, and Gaétan Lantagne. A study on the effect of porosity and particles size distribution on Li-ion battery per- formance. Journal of The Electrochemical Society, 164(11):E3179–E3189, 2017. doi: 10.1149/2.0211711jes. [2] Indrajeet V. Thorat, David E. Stephenson, Nathan A. Zacharias, Karim Zaghib, John N. Harb, and Dean R. Wheeler. Quantifying tortuosity in porous Li- ion battery materials. Journal of Power Sources, 188(2):592–600, 2009. doi: 10.1016/j.jpowsour.2008.12.032. [3] Huanqiao Song, Jiangang Li, Mingsheng Luo, Qiannan Zhao, and Feng Liu. Ultra- thin mesoporous LiV3 O8 nanosheet with exceptionally large specific area for fast and reversible Li storage in lithium-ion battery cathode. Journal of The Electrochemical Society, 168(5):050515, 2021. doi: 10.1149/1945-7111/abfca1. [4] Taira Aida, Takahiro Toma, and Satoshi Kanada. A comparative study of particle size and hollowness of LiNi1/3 Co1/3 Mn1/3 O2 cathode materials for high-power Li-ion bat- teries: effects on electrochemical performance. Journal of Solid State Electrochemistry, 24(6):1415–1425, 2020. doi: 10.1007/s10008-020-04640-z. [5] Yahong Xu, Enyuan Hu, Kai Zhang, Xuelong Wang, Valery Borzenets, Zhihong Sun, Piero Pianetta, Xiqian Yu, Yijin Liu, Xiao Qing Yang, and Hong Li. In situ visu- alization of state-of-charge heterogeneity within a LiCoO2 particle that evolves upon cycling at different rates. ACS Energy Letters, 2(5):1240–1245, 2017. doi: 10.1021/ac- senergylett.7b00263. [6] Mark E. Orazem and Bernard Tribollet. Electrochemical Impedance Spectroscopy. John Wiley & Sons, Inc, Hoboken, New Jersey, 2008. [7] Yuejiu Zheng, Zhihe Shi, Dongxu Guo, Haifeng Dai, and Xuebing Han. A simplification of the time-domain equivalent circuit model for lithium-ion batteries based on low- frequency electrochemical impedance spectra. Journal of Power Sources, 489:229505, 2021. doi: 10.1016/j.jpowsour.2021.229505. [8] Andreas H. Wiedemann, Graham M. Goldin, Scott A. Barnett, Huayang Zhu, and Robert J. Kee. Effects of three-dimensional cathode microstructure on the perfor- mance of lithium-ion battery cathodes. Electrochimica Acta, 88:580–588, 2013. doi: 10.1016/j.electacta.2012.10.104. [9] Bo Yan, Cheolwoong Lim, Leilei Yin, and Likun Zhu. Three dimensional simulation of galvanostatic discharge of LiCoO2 cathode based on X-ray nano-CT images. Journal of The Electrochemical Society, 159(10):A1604–A1614, 2012. doi: 10.1149/2.024210jes. 166 [10] Benedikt Prifling, Daniel Westhoff, Denny Schmidt, Henning Markötter, Ingo Manke, Volker Knoblauch, and Volker Schmidt. Parametric microstructure modeling of com- pressed cathode materials for Li-ion batteries. Computational Materials Science, 169, 2019. doi: 10.1016/j.commatsci.2019.109083. [11] B. E. Conway. Transition from “supercapacitor” to “battery” behavior in electrochem- ical energy storage. Journal of The Electrochemical Society, 138(1539), 1991. [12] John R. Miller. Engineering electrochemical capacitor applications. Journal of Power Sources, 326(15):726–735, 2016. doi: https://doi.org/10.1016/j.jpowsour.2016.04.020. [13] U.S. Department of Energy. Basic research needs for electrical energy storage., 2007. URL https://www.osti.gov/servlets/purl/935429. [14] B. E. Conway. Electrochemical Supercapacitors: scientific fundamentals and technolog- ical applications. Springer US, 1999. [15] P. Simon and Y. Gogotsi. Materials for electrochemical capacitor. Nature Materials, 7(11):845–854, 2008. [16] A. M. Johnson and John Newman. Desalting by means of porous carbon electrodes. Journal of the Electrochemical Society, 118:510, 1971. [17] Lei Zhang, Xiaosong Hu, Zhenpo Wang, Fengchun Sun, and David G. Dorrell. A re- view of supercapacitor modeling, estimation, and applications: A control/management perspective. Renewable and Sustainable Energy Reviews, 81:1868–1878, 2018. [18] A. M. Johnson and John Newman. Desalting by means of porous carbon electrodes. Journal of the Electrochemical Society, 118:510, 1971. [19] Banxian Ruan, Qi You, Jiaqi Zhu, Leiming Wu, Jun Guo, Xiaoyu Dai, and Yuanjiang Xiang. Terahertz biochemical sensor based on strong coupling between waveguide mode and surface plasmons of double layer graphene. IEEE Sensors Journal, 18(18): 7436–7441, 2018. [20] R. E. White G. Sikha and B. N. Popov. A mathematical model for a lithium-ion bat- tery/electrochemical capacitor hybird system. Journal of the Electrochemical Society, 152(8):1682–1693, 2005. [21] Alexander Urban, Dong-Hwa Seo, and Gerbrand Ceder. Computational understanding of li-ion batteries. NPJ Computational Materials, 2:1–13, 2016. [22] Arnulf Latz and Jochen Zausch. Multiscale modeling of lithium ion batteries: thermal aspects. Beilstein Journal of nanotechnology, 6:987–1007, 2015. 167 [23] Gaetan Patry, Alex Romagny, and Sebastien Martinet. Cost modeling of lithium-ion battery cells for automotive applications. Energy Science and Engineering, 3:71–82, 2015. [24] Kevin Leung, Yue Qi, Kevin R. Zavadil, Yoon Seok Jung, Anne C. Dillon, Andrew S. Cavanagh, Se-Hee Lee, and Steven M. George. Using atomic layer deposition to hinder solvent decomposition in lithium ion batteries: First-principles modeling and experi- mental studies. Journal of the American Chemical Society, 133:14741–14754, 2011. [25] H. L. F. von Helmholtz. Studien uer electrische grenzschichten. Annalen der Physik, 243:337–382, 1879. [26] Hainan Wang. Modeling and simulations of electrical energy storage in electrochemical capacitors. PhD thesis, University of California Los Angeles, 2013. [27] J. J. Bikerman. Structure and capacity of electrical double layer. Philosophical Mag- azine, 33(220):384–397, 1942. [28] Martin Z. Bazant, Mustafa Sabri Kilic, Brian D. Storey, and Armand Ajdari. To- wards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions. Advances in Colloid and Interface Science, 152:48–88, 2009. [29] Jacob H. Masliyah and Subir Bhattacharjee. Electrokinetic and colloid transport phe- nomena. John Wiley Sons, Inc., 2006. [30] W. G. Pell and B. E. Conway. Voltammetry at a de levie brush electrode as a model for electrochemical supercapacitor behaviour. Journal of Electroanalytical Chemistry, 500(1-2):121–133, 2001. [31] W. G. Pell and B. E. Conway. Analysis of power limitations at porous supercapacitor electrodes under cyclic voltammetry modulation and dc charge. Journal of Power Sources, 96(1):57–67, 2001. [32] S. Yoon, J. H. Jang, B. H. Ka, and S. M. Oh. Complex capacitance analysis on rate capability of electric-double layer capacitor (edlc) electrodes of different thckness. Electrochemical Acta, 50(11):2255–2262, 2005. [33] Douglas Henderson, Stanislaw Lamperski, Zhehui Jin, and JIanzhong Wu. Density functional study of the electric double layer formed by a high density electrolyte. Journal of Physical Chemistry, 115:12911–12914, 2011. [34] G. M. Torrie, J. P. Valleau, and G. N. Patey. Electrical double layers Monte Carlo and HNC studies of image effects. Journal of Chemical Physics, 76:4615, 1982. [35] Timea Nagy, Douglas Henderson, and Dezso Boda. Simulation of an electrical dou- 168 ble layer model with a low dielectric layer between the electrode and the electrolyte. Journal of Physical Chemistry, 115:11409–11419, 2011. [36] Lutful Bari Bhuiyan, Stanislaw Lamperski, Jianzhong Wu, and Douglas Henderson. Monte Carlo simulation for the double layer structure of an ionic liquid using a dimer model: a comparison with the density function theory. Journal of Chemical Physics, 116:10364–10370, 2012. [37] William Tiedemann and John Newman. Double-layer capacity determination of porous electrodes. Journal of the Electrochemical Society, 122:70–74, 1975. [38] Todd R. Ferguson and Martin Z. Bazant. Nonequilibrium thermodynamics of porous electrodes. Journal of the Electrochemical Society, 159(12):A1967–A1985, 2012. [39] P. M. Biesheuvel, Y. Fu, and M. Z. Bazant. Electrochemistry and capacitive charging of porous electrodes in asymmetric multicomponent electrolytes. Russian Journal of Electrochemistry, 48(6):580–592, 2012. [40] Yuki Kitazumi, Osamu Shirai, Masahiro Yamamoto, and Kenji Kano. Numerical sim- ulation of diffuse double layer around microporous electrodes based on the poisson- boltzmann equation. Electrochimica Acta, 112:171–175, 2013. [41] Hainan Wang and Laurent Pilon. Accurate simulations of electric double layer ca- pacitance of ultramicroelectrodes. Journal of Physical Chemistry, 115:16711–16719, 2011. [42] Martin Z. Bazant, Katsuyo Thornton, and Armand Ajdari. Diffuse-charge dynamics in electrochemical systems. Physical Review E, 70:021506, 2004. [43] Juchuan Li, Fuqian Yang, Xingcheng Xiao, Mark W Verbrugge, and Yang-Tse Cheng. Potentiostatic intermittent titration technique ( pitt ) for spherical particles with finite interfacial kinetics. Electrochimica Acta, 75:56–61, 2012. [44] Minyu Jia, Wenheng Zhang, Xiaoping Cai, Xinju Zhan, Linrui Hou, Changzhou Yuan, and Zaiping Guo. Re-understanding the galvanostatic intermittent titration technique : Pitfalls in evaluation of diffusion coefficients and rational suggestions. Journal of Power Sources, 543:231843, 2022. ISSN 0378-7753. doi: 10.1016/j.jpowsour.2022.231843. [45] Cameron Day, Katie Greig, Alexander Massey, Jennifer Peake, David Crossley, and Robert A W Dryfe. Utilizing cyclic voltammetry to understand the energy storage mechanisms for copper oxide and its graphene oxide hybrids as lithium-ion battery anodes. ChemSusChem, 13:1504–1516, 2020. doi: 10.1002/cssc.201902784. [46] Stefan Skoog and Sandeep David. Parameterization of linear equivalent circuit models over wide temperature and soc spans for automotive lithium-ion cells using electro- 169 chemical impedance spectroscopy. Journal of Energy Storage, 14:39–48, 2017. doi: 10.1016/j.est.2017.08.004. [47] Mohammed Ahmed Zabara, Can Berk Uzundal, and Burak Ulgut. Linear and nonlinear electrochemical impedance spectroscopy studies of Li/SOCl2 batteries. Journal of The Electrochemical Society, 166(6):A811–A820, 2019. doi: 10.1149/2.1231904jes. [48] PalmSens BV. Bode and nyquist plot., 2022. URL https://www.palmsens.com/knowledgebase-article/bode-and-nyquist-plot/. [49] Pietro Iurilli, Claudio Brivio, and Vanessa Wood. On the use of electrochemical impedance spectroscopy to characterize and model the aging phenomena of lithium- ion batteries: a critical review. Journal of Power Sources, 505:229860, 2021. doi: 10.1016/j.jpowsour.2021.229860. [50] Tapesh Joshi, KwangSup Eom, Gleb Yushin, and Thomas F. Fuller. Effects of dissolved transition metals on the electrochemical performance and sei growth in lithium-ion batteries. Journal of The Electrochemical Society, 161(12):A1915–A1921, 2014. doi: 10.1149/2.0861412jes. [51] Xi Zhang, Yizhao Gao, Bangjun Guo, Chong Zhu, Xuan Zhou, Lin Wang, and Jianhua Cao. A novel quantitative electrochemical aging model considering side reactions for lithium-ion batteries. Electrochimica Acta, 343:136070, 2020. doi: 10.1016/j.electacta.2020.136070. [52] Tianhan Gao, Andrew Kim, and Wei Lu. Modeling electrode-level crack and quantify- ing its effect on battery performance and impedance. Electrochimica Acta, 363:137197, 2020. doi: 10.1016/j.electacta.2020.137197. [53] Martin Bettge, Yan Li, Kevin Gallagher, Ye Zhu, Qingliu Wu, Wenquan Lu, Ira Bloom, and Daniel P. Abraham. Voltage fade of layered oxides: Its measurement and impact on energy density. Journal of The Electrochemical Society, 160(11):A2046–A2055, 2013. doi: 10.1149/2.034311jes. [54] Debasish Mohanty, Athena S. Sefat, Jianlin Li, Roberta A. Meisner, Adam J. Rondi- none, E. Andrew Payzant, Daniel P. Abraham, David L. Wood, and Claus Daniel. Correlating cation ordering and voltage fade in a lithium-manganese-rich lithium-ion battery cathode oxide: A joint magnetic susceptibility and tem study. Physical Chem- istry Chemical Physics, 15(44):19496–19509, 2013. doi: 10.1039/C3CP53658K. [55] John Newman and William Tiedemann. Porous-electrode theory with battery appli- cations. AIChE Journal, 21:25–41, 1975. doi: 10.1002/aic.690210103. [56] John S. Newman. Electrochemical Systems. Prentice Hall, Englewood Cliffs, New Jersey 07632, 2nd edition, 1991. 170 [57] David E. Stephenson, Erik M. Hartman, John N. Harb, and Dean R. Wheeler. Model- ing of particle-particle interactions in porous cathodes for lithium-ion batteries. Journal of The Electrochemical Society, 154:A1146, 2007. doi: 10.1149/1.2783772. [58] Venkat Srinivasan and John Newman. Discharge model for the lithium iron- phosphate electrode. Journal of The Electrochemical Society, 151:A1517, 2004. doi: 10.1149/1.1785012. [59] Thanh Son Dao, Chandrika P. Vyasarayani, and John McPhee. Simplification and order reduction of lithium-ion battery model based on porous-electrode theory. Journal of Power Sources, 198:329–337, 2012. doi: 10.1016/j.jpowsour.2011.09.034. [60] Gen Inoue, Hiroki Mashioka, Naoki Kimura, and Yoshifumi Tsuge. Identifying pa- rameters from discharging and relaxation curves of lithium-ion batteries using porous electrode theory. Journal of Chemical Engineering of Japan, 54:207–212, 2021. doi: 10.1252/jcej.20we180. [61] Scott A. Roberts, Victor E. Brunini, Kevin N. Long, and Anne M. Grillet. A frame- work for three-dimensional mesoscale modeling of anisotropic swelling and mechanical deformation in lithium-ion electrodes. Journal of The Electrochemical Society, 161(11): F3052–F3059, 2014. doi: 10.1149/2.0081411jes. [62] Christian Wieser, Torben Prill, and Katja Schladitz. Multiscale simulation process and application to additives in porous composite battery electrodes. Journal of Power Sources, 277:64–75, 2015. doi: 10.1016/j.jpowsour.2014.11.090. [63] S J Cooper, D S Eastwood, J Gelb, G Damblanc, D J L Brett, R S Bradley, P J Withers, P D Lee, A J Marquis, N P Brandon, and P R Shearing. Image based modelling of microstructural heterogeneity in LiFePO 4 electrodes for Li-ion batteries. Journal of Power Sources, 247:1033– 1039, 2014. ISSN 0378-7753. doi: 10.1016/j.jpowsour.2013.04.156. URL http://dx.doi.org/10.1016/j.jpowsour.2013.04.156. [64] Xuekun Lu, Antonio Bertei, Donal P. Finegan, Chun Tan, Sohrab R. Daemi, Julia S. Weaving, Kieran B. O’Regan, Thomas M.M. Heenan, Gareth Hinds, Emma Kendrick, Dan J.L. Brett, and Paul R. Shearing. 3D microstructure design of lithium-ion bat- tery electrodes assisted by X-ray nano-computed tomography and modelling. Nature Communications, 11(1):1–13, 2020. doi: 10.1038/s41467-020-15811-x. [65] Chia-Wei Wang and Ann Marie Sastry. Mesoscale modeling of a Li-ion polymer cell. Journal of The Electrochemical Society, 154:A1035, 2007. doi: 10.1149/1.2778285. [66] Vittorio De Lauri, Lukas Krumbein, Simon Hein, Benedikt Prifling, Volker Schmidt, Timo Danner, and Arnulf Latz. Beneficial effects of three-dimensional structured elec- trodes for the fast charging of lithium-ion batteries. ACS Applied Energy Materials, 4 171 (12), 2021. ISSN 2574-0962. doi: 10.1021/acsaem.1c02621. [67] J. Madison, J. Spowart, D. Rowenhorst, L. K. Aagesen, K. Thornton, and T. M. Pollock. Modeling fluid flow in three-dimensional single crystal dendritic structures. Acta Materialia, 58:2864–2875, 2010. doi: 10.1016/j.actamat.2010.01.014. [68] Abbos Shodiev, Emiliano N. Primo, Mehdi Chouchane, Teo Lombardo, Alain C. Ngandjong, Alexis Rucci, and Alejandro A. Franco. 4D-resolved physical model for electrochemical impedance spectroscopy of Li(Ni1−x−y Mnx Coy )O2 -based cathodes in symmetric cells: Consequences in tortuosity calculations. Journal of Power Sources, 454:227871, 2020. doi: 10.1016/j.jpowsour.2020.227871. [69] Mehdi Chouchane, Emiliano N. Primo, and Alejandro A. Franco. Mesoscale effects in the extraction of the solid-state lithium diffusion coefficient values of battery active materials: Physical insights from 3D modeling. Journal of Physical Chemistry Letters, 11(7):2775–2780, 2020. doi: 10.1021/acs.jpclett.0c00517. [70] Hui Chia Yu, Hsun Yi Chen, and K. Thornton. Extended smoothed boundary method for solving partial differential equations with general boundary conditions on complex boundaries. Modelling and Simulation in Materials Science and Engineering, 20(7): 1–41, 2012. ISSN 09650393. doi: 10.1088/0965-0393/20/7/075008. [71] Christoph Paul Schmidt, Stephan Sinzig, Volker Gravemeier, and Wolfgang A. Wall. A three-dimensional finite element formulation coupling electrochemistry and solid mechanics on resolved microstructures of all-solid-state lithium-ion batteries. 2022. URL https://ssrn.com/abstract=4189627. [72] Robert Termuhlen, Kieran Fitzmaurice, and Hui-Chia Yu. Smoothed boundary method for simulating incompressible flow in complex geometries. Computer Methods in Ap- plied Mechanics and Engineering, 399:115312, 2022. doi: 10.1016/j.cma.2022.115312. [73] Danqi Qu, Robert Termuhlen, and Hui-Chia Yu. Simulation of electrochemical double layer formation with complex geometries. Journal of The Electrochemical Society, 167: 140515, 2020. doi: 10.1149/1945-7111/abc0ab. [74] Chohong Min, Frédéric Gibou, and Hector D. Ceniceros. A supra-convergent fi- nite difference scheme for the variable coefficient Poisson equation on non-graded grids. Journal of Computational Physics, 218(1):123–140, 2006. ISSN 10902716. doi: 10.1016/j.jcp.2006.01.046. [75] Han Chen, Chohong Min, and Frédéric Gibou. A supra-convergent finite difference scheme for the poisson and heat equations on irregular domains and non-graded adap- tive cartesian grids. Journal of Scientific Computing, 31(1-2):19–60, 2007. ISSN 08857474. doi: 10.1007/s10915-006-9122-8. 172 [76] Jessica Luck and Arnulf Latz. The electrochemical double layer and its impedance behavior in lithium-ion batteries. Physical Chemistry Chemical Physics, 21:14753– 14765, 2019. [77] H. Jia, Y. Fu, Y. Zhang, and W. He. Design of hybrid energy storage control system for wind farms based on flow battery and electric double-layer capacitor. 2010 Asia-Pacific Power and Energy Engineering Conference, pages 1–6, 2010. [78] Ning Liu, Ru Chen, and Qing Wan. Recent advances in electric-double-layer transistors for bio-chemical sensing applications. Sensors, 19(15), 2019. ISSN 1424-8220. doi: 10.3390/s19153425. URL https://www.mdpi.com/1424-8220/19/15/3425. [79] John S. Newman. Electrochemical systems. Prentice Hall, Englewood Cliffs, New Jersey 07632, 2nd edition, 1991. [80] Elcap. Simplified illustration of the potential development in the area and in the further course of a helmholtz double layer., 2012. URL https://en.wikipedia.org/wiki/Double layer (surface science). [81] Hainan Wang, Julian Varghese, and Laurent Pilon. Simulation of electric double layer capacitors with mesoporous electrodes: Effects of morphology and electrolyte permit- tivity. Electrochimica Acta, 56:6189–6197, 2011. [82] Jan W. Slotboom. Computer-aided two-dimensional analysis of bipolar transistors. IEEE Transactions on Electron Devices, 20(8):669–679, 1973. [83] Richard S. Varga. Matrix Iterative Analysis. Englewood Cliffs, 1962. [84] Jun Huang. Diffusion impedance of electroactive materials, electrolytic solutions and porous electrodes: Warburg impedance and beyond. Electrochimica Acta, 281:170–188, 2018. [85] C. Capiglia, Y. Saito, H. Kageyama, P. Mustarelli, T. Iwamoto, and T. Tabuchi. Li and F diffusion coefficients and thermal properties of non-aqueous eletrolyte solutions for rechargeable lithium batteries. Journal of Power Sources, 81-82:859–862, 1999. [86] Lars Ole Valoen and Jan N. Reimers. Transport properties of LiPF6 -based li-ion battery electrolytes. Journal of The Electrochemical Society, 152(5):A882–A891, 2005. [87] Bernardo Orvananos Murguia. Modeling and simulation of nanoparticulate lithium ion phosphate battery electrodes. PhD thesis, University of Michigan, 2014. [88] David S. Hall, Julian Self, and J. R. Dahn. Dielectric constants for quantum chemistry and li-ion batteries: solvent blends of ethylene carbonate and ethyl methyl carbonate. The Journal of Physical Chemistry, 119:22322–22330, 2015. 173 [89] H. Zhu, A. Ghoufi, and A. Szymczyk. Anomalous dielectric behavior of nanoconfined electrolytic solutions. Physical review letters, 109:107801 – 5, 2012. [90] Martin Z. Bazant, Katsuyo Thornton, and Armand Ajdari. Diffuse-charge dynamics in electrochemical systems. Physical Review E, 70:021506, 2004. [91] G. Y. Chen, H. K. Lin, and C. W. Lan. Phase-field modeling of twin-related faceted dendrite growth of silicon. Acta Materialia, 115:324–332, 2016. [92] H. K. Lin, C. C. Chen, and C. W. Lan. A simple anisotropic surface free energy function for three-dimensional phase field modeling of multi-crystalline crystal growth. Journal of Crystal Growth, 362:62–65, 2013. [93] Affan Malik and Hui-Chia Yu. Complex electrode microstructure simulations using a smoothed boundary method with adaptive mesh refinement. Journal of the Electro- chemical Society, 169:070527, 2022. doi: 10.1149/1945-7111/ ac7e79. [94] Ming Tang, James F Belak, and Milo R Dorr. Anisotropic phase boundary morphology in nanoscale olivine electrode particles. The Journal of Physical Chemistry C, 115:4922 – 4926, 2011. doi: 10.1021/jp109628m. [95] S Dargaville and T W Farrell. A comparison of mathematical models for phase- change in high-rate LiFePO4 cathodes. Electrochimica Acta, 111:474 – 490, 2013. doi: 10.1016/j.electacta.2013.08.014. [96] Daniel A Cogswell and Martin Z Bazant. Theory of coherent nucleation in phase- separating nanoparticles. Nano Letters, 13:3036 – 3041, 2013. doi: 10.1021/nl400497t. [97] Liang Hong, Linsen Li, Yuchen-Karen Chen-Wiegart, Jiajun Wang, Kai Xiang, Liyang Gan, Wenjie Li, Fei Meng, Fan Wang, Jun Wang, Yet-Ming Chiang, Song Jin, and Ming Tang. Two-dimensional lithium diffusion behavior and probable hybrid phase transformation kinetics in olivine lithium iron phosphate. Nature Communications, 8: 114, 2017. doi: 10.1038/s41467-017-01315-8. [98] Anis Allagui, Ahmed S. Elwakil, and Hichem Eleuch. Highlighting a common confusion in the computation of capacitance of electrochemical energy storage devices. The Jour- nal of Physical Chemistry C, 125(18):9591–9592, 2021. doi: 10.1021/acs.jpcc.1c01288. [99] Bernardo Orvananos, Todd R. Ferguson, Hui-Chia Yu, Martin Z. Bazant, and Katsuyo Thornton. Particle-level modeling of the charge-discharge behavior of nanoparticulate phase-separating Li-ion battery electrodes. Journal of The Electrochemical Society, 161:A535–A546, 2014. doi: 10.1149/2.024404jes. [100] D P Abraham, S Kawauchi, and D W Dees. Modeling the impedance versus voltage characteristics of LiNi0.8 Co0.15 Al0.05 O2 . Electrochimica Acta, 53:2121 – 2129, 2008. doi: 174 10.1016/j.electacta.2007.09.018. [101] Stefan Luding. Introduction to discrete element methods. European Jour- nal of Environmental and Civil Engineering, 12(7-8):785–826, 2008. doi: 10.1080/19648189.2008.9693050. [102] Yrj´’o Jun Huang, Ole J{o}rgen Nydal, Chenhui Ge, and Baodian Yao. An intro- duction to discrete element method: A meso-scale mechanism analysis of granular flow. Journal of Dispersion Science and Technology, 36:1370 – 1377, 2015. doi: 10.1080/01932691.2014.984304. [103] C L Park, P W Voorhees, and K Thornton. Application of the level-set method to the analysis of an evolving microstructure. Computational Materials Science, 85(C), 2014. doi: 10.1016/j.commatsci.2013.12.022. [104] Ruhul Amin and Yet-Ming Chiang. Characterization of electronic and ionic transport in Li1−x Ni0.33 Mn0.33 Co0.33 O2 (NMC333) and Li1−x Ni0.50 Mn0.20 Co0.30 O2 (NMC523) as a function of Li content. Journal of The Electrochemical Society, 163:A1512–A1517, 2016. doi: 10.1149/2.0131608jes. [105] Naoki Nitta, Feixiang Wu, Jung Tae Lee, and Gleb Yushin. Li-ion bat- tery materials: present and future. Materials Today, 18:252–264, 2015. doi: 10.1016/j.mattod.2014.10.040. [106] Ines Baccouche, Sabeur Jemmali, Bilal Manai, Noshin Omar, and Najoua Essoukri Ben Amara. Improved OCV model of a Li-ion NMC battery for online SOC estimation using the extended Kalman filter. Energies, 10:764, 2017. doi: 10.3390/en10060764. [107] Ping-Chun Tsai, Bohua Wen, Mark Wolfman, Min-Ju Choe, Menghsuan Sam Pan, Liang Su, Katsuyo Thornton, Jordi Cabana, and Yet-Ming Chiang. Single-particle measurements of electrochemical kinetics in NMC and NCA cathodes for Li-ion bat- teries. Energy & Environmental Science, 11:860–871, 2018. doi: 10.1039/c8ee00001h. [108] Lars Ole Valøena and Jan N. Reimers. Transport properties of LiPF6 -based Li-ion battery electrolytes. Journal of The Electrochemical Society, 152:A882–A891, 2005. doi: 10.1149/1.1872737. [109] G Horvai. Relationship between charge transfer resistance and exchange current density of ion transfer at the interface of two immiscible electrolyte solutions. Electroanalysis, 3:673–675, 1991. doi: 10.1002/elan.1140030713. [110] H. Zhu, A. Ghoufi, A. Szymczyk, B. Balannec, and D. Morineau. Anomalous dielectric behavior of nanoconfined electrolytic solutions. Physical Review Letters, 109:107801, Sep 2012. doi: 10.1103/PhysRevLett.109.107801. 175 [111] Hui-Chia Yu, Stuart B Adler, Scott A Barnett, and K Thornton. Simulation of the diffusional impedance and application to the characterization of electrodes with complex microstructures. Electrochimica Acta, 354:136534, 09 2020. doi: 10.1016/j.electacta.2020.136534. [112] S. R. Sivakkumar, J. Y. Nerkar, and A. G. Pandolfo. Rate capability of graphite materials as negative electrodes in lithium-ion capacitors. Electrochimica Acta, 55(9): 3330–3335, 2010. ISSN 00134686. doi: 10.1016/j.electacta.2010.01.059. [113] Yulong Fu, Yuqing Jin, Jing Ma, Junhao Liu, Zhi Wang, Bin Wang, and Xuzhong Gong. Lithium-ion transfer strengthened by graphite tailings and coking coal for high- rate performance anode. Chemical Engineering Journal, 442(P1):136184, 2022. ISSN 1385-8947. doi: 10.1016/j.cej.2022.136184. [114] Killian R Tallman, Bingjie Zhang, Lei Wang, Shan Yan, Katherine Thompson, Xiao Tong, Juergen Thieme, Andrew Kiss, Amy C Marschilok, Kenneth J Takeuchi, David C Bock, and Esther S Takeuchi. Anode overpotential control via interfacial modification: Inhibition of lithium plating on graphite anodes. Applied Materials and Interfaces, 2019. doi: 10.1021/acsami.9b16794. [115] Yinsheng Guo, Raymond B. Smith, Zhonghua Yu, Dmitri K. Efetov, Junpu Wang, Philip Kim, Martin Z. Bazant, and Louis E. Brus. Li intercalation into graphite: Direct optical imaging and cahn-hilliard reaction dynamics. Journal of Physical Chemistry Letters, 7(11):2151–2156, 2016. ISSN 19487185. doi: 10.1021/acs.jpclett.6b00625. [116] Nicholas Loeffler, Dominic Bresser, Stefano Passerini, and Mark Copley. Secondary lithium-ion battery anodes: From first commercial batteries to recent research activ- ities. Johnson Matthey Technology Review, 59(1):34–44, 2015. ISSN 20565135. doi: 10.1595/205651314X685824. [117] Wilhelm Pfleging. A review of laser electrode processing for development and manu- facturing of lithium-ion batteries. Nanophotonics, 7(3), 2018. doi: 10.1515/nanoph- 2017-0044. [118] Kuan-Hung Chen, Min Ji Namkoong, Vishwas Goel, Chenglin Yang, Saeed Kazemiab- navi, S.M. Mortuza, Eric Kazyak, Jyoti Mazumder, Katsuyo Thornton, Jeff Sakamoto, and Neil P. Dasgupta. Efficient fast-charging of lithium-ion batteries enabled by laser- patterned three-dimensional graphite anode architectures. Journal of Power Sources, 471:228475, 2020. doi: 10.1016/j.jpowsour.2020.228475. [119] Tao Gao, Yu Han, Dimitrios Fraggedakis, Supratim Das, Tingtao Zhou, Che Ning Yeh, Shengming Xu, William C. Chueh, Ju Li, and Martin Z. Bazant. Interplay of lithium intercalation and plating on a single graphite particle. Joule, 5(2):393–414, 2021. ISSN 25424351. doi: 10.1016/j.joule.2020.12.020. 176 [120] Moses Ender, Jochen Joos, André Weber, and Ellen Ivers-Tiffée. Anode microstruc- tures from high-energy and high-power lithium-ion cylindrical cells obtained by X- ray nano-tomography. Journal of Power Sources, 269(c):912 – 919, 12 2014. doi: 10.1016/j.jpowsour.2014.07.070. [121] Tsutomu Ohzuku, Yasunobu Iwakoshi, and Keijiro Sawai. Formation of lithium- graphite intercalation compounds in nonaqueous electrolytes and their applica- tion as a negative electrode for a lithium ion (shuttlecock) cell. Journal of The Electrochemical Society, 140(9):2490, sep 1993. doi: 10.1149/1.2220849. URL https://dx.doi.org/10.1149/1.2220849. [122] Rachid Yazami, Audrey Martinent, and Yvan Reynier. Some thermodynamics and kinetics aspects of the graphite-lithium negative electrode for lithium-ion batteries. pages 245–258, 2006. doi: 10.1007/1-4020-4812-2 18. [123] Michael Peter Mercer, Chao Peng, Cindy Soares, Harry Ernst Hoster, and Denis Kramer. Voltage hysteresis during lithiation/delithiation of graphite associated with meta-stable carbon stackings. Journal of Materials Chemistry A, 9(1):492–504, 2021. ISSN 20507496. doi: 10.1039/d0ta10403e. [124] Affan Mailk, Kent Snyder, Minghong Liu, and Hui-Chia Yu. Electrochemical dy- namics in hybrid graphite–carbon electrodes. MRS Communications, 2022. doi: 10.1557/s43579-022-00214-4. [125] Affan Mailk, Kent Snyder, Minghong Liu, and Hui-Chia Yu. Microstructure-level phase-field electrochemical simulations of graphite complex electrodes. in preparation, 2022. [126] Simon Müller. X-ray tomography data of four commercial lithium ion bat- tery graphite electrodes: Research data supporting “quantifying inhomogeneity of lithium ion battery electrodes and its influence on electrochemical performance”. https://www.research-collection.ethz.ch/handle/20.500.11850/224851, 2018. [127] Simon Müller, Jens Eller, Martin Ebner, Chris Burns, Jeff Dahn, and Vanessa Wood. Quantifying inhomogeneity of lithium ion battery electrodes and its influence on elec- trochemical performance. Journal of The Electrochemical Society, 165:A339, 2018. [128] Danqi Qu, Affan Mailk, and Hui-Chia Yu. Physics-based simulation of electrochemical impedance spectroscopy of complex electrode microstructures using smoothed bound- ary method. Electrochimica Acta, 2022. doi: 10.1016/j.electacta.2022.141141. [129] Chang-Jun Bae, Can K Erdonmez, John W Halloran, and Yet-Ming Chiang. Design of battery electrodes with dual-scale porosity to minimize tortuosity and maximize perfor- mance. Advanced Materials, 25(9):1254 – 1258, 12 2012. doi: 10.1002/adma.201204055. 177 [130] Jan B Habedank, Joseph Endres, Patrick Schmitz, Michael F Zaeh, and Heinz P Huber. Femtosecond laser structuring of graphite anodes for improved lithium-ion batteries: Ablation characteristics and process design. Journal of Laser Applications, 30(3):032205 – 7, 2018-08. doi: 10.2351/1.5040611. [131] David Linden and Thomas B Reddy. Handbook of Batteries, chapter 2. McGraw-Hill Professional, 3rd edition, 1995. [132] Wei Lai. Electrochemical modeling of single particle intercalation battery materials with different thermodynamics. Journal of Power Sources, 196(15):6534 – 6553, 2011. doi: 10.1016/j.jpowsour.2011.03.055. [133] A Van der Ven and G Ceder. Lithium diffusion in layered Lix CoO2 . Electrochemical and Solid-State Letters, 3(7):301 – 304, 2000. [134] P J Bouwman, B A Boukamp, H J M Bouwmeester, and P H L Notten. Influence of diffusion plane orientation on electrochemical properties of thin film LiCoO2 electrodes. Journal of The Electrochemical Society, 149(6):A699, 2002. doi: 10.1149/1.1471543. [135] Ruhul Amin and Yet-Ming Chiang. Characterization of electronic and ionic transport in Li1−x Ni0.33 Mn0.33 Co0.33 O2 (NMC333) and Li1−x Ni0.50 Mn 0.20 Co0.30 O2 (NMC523) as a function of Li content. Journal of The Electrochemical Society, 163(8):A1512 – A1517, 2016. doi: 10.1149/2.0131608jes. [136] Asif Islam Khan, Korok Chatterjee, Brian Wang, Steven Drapcho, Long You, Claudy Serrao, Saidur Rahman Bakaul, Ramamoorthy Ramesh, and Sayeef Salahuddin. Neg- ative capacitance in a ferroelectric capacitor. Nature Materials, 14(2):182–186, 2015. ISSN 14764660. doi: 10.1038/nmat4148. [137] Wei Wu, Hailong Yin, Hao Zhang, Jia Kang, Yun Li, and Yong Dan. Electrochemical investigation of corrosion of X80 steel under elastic and plastic tensile stress in CO2 environment. Metals, 8(11), 2018. ISSN 20754701. doi: 10.3390/met8110949. [138] Juan Bisquert and Antonio Guerrero. Chemical inductor. Journal of the Ameri- can Chemical Society, 144(13):5996–6009, 2022. doi: 10.1021/jacs.2c00777. PMID: 35316040. [139] Klemen Zelič, Igor Mele, Arghya Bhowmik, and Tomaž Katrašnik. Phase sep- arating electrode materials – chemical inductors. arXiv/cond-mat.mtrl-sci, page arXiv:2206.08089, 2022. doi: 10.1021/jacs.2c00777. 178