# **Oleg Birkholz**

MODELING TRANSPORT PROPERTIES AND ELECTROCHEMICAL PERFORMANCE OF HIERARCHICALLY STRUCTURED LITHIUM-ION BATTERY CATHODES USING RESISTOR NETWORKS AND MATHEMATICAL HALF-CELL MODELS

SCHRIFTENREIHE DES INSTITUTS FÜR ANGEWANDTE MATERIALIEN

**O. BIRKHOLZ**

Modeling hierarchically structured lithium-ion battery cathodes

BAND 105

Oleg Birkholz

**Modeling transport properties and electrochemical performance of hierarchically structured lithium-ion battery cathodes using resistor networks and mathematical half-cell models**

## **Schriftenreihe des Instituts für Angewandte Materialien** *Band 105*

Karlsruher Institut für Technologie (KIT) Institut für Angewandte Materialien (IAM)

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buches.

**Modeling transport properties and electrochemical performance of hierarchically structured lithium-ion battery cathodes using resistor networks and mathematical half-cell models**

by Oleg Birkholz

Karlsruher Institut für Technologie Institut für Angewandte Materialien

Modeling transport properties and electrochemical performance of hierarchically structured lithium-ion battery cathodes using resistor networks and mathematical half-cell models

Zur Erlangung des akademischen Grades eines Doktors der Ingenieurwissenschaften (Dr.-Ing.) von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Oleg Birkholz, M.Sc.

Tag der mündlichen Prüfung: 2. Juni 2021 Hauptreferent: Prof. Dr. Marc Kamlah Korreferent: Prof. Dr. Michael J. Hoffmann

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2022 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 2192-9963 ISBN 978-3-7315-1172-4 DOI 10.5445/KSP/1000143556

*Für meine Eltern*

# **Abstract**

Hierarchically structured active materials in lithium-ion battery (LIB) electrodes with porous secondary particles are promising candidates to increase gravimetric energy density and rate performance of the cell. However, there are still aspects to this technology which are not fully understood. In order to obtain deeper knowledge, the goal of this work is to develop efficient tools to compute effective transport properties of granular cathode structures as well as porous secondary particles and to evaluate their influence on the electrochemical performance of the whole cell.

On the one hand, the resistor network method (RN)—a tool to efficiently compute the effective transport properties of particulate systems—is extended with regard to the transport through the solid and the pore phase of granular media represented by sphere packings with polydisperse size-distributions. As for the solid phase, transport through the volume, via the surface, or a mix of both are considered. For all cases, appropriate analytically derived formulas from literature describing resistance between two single particles are used or combined accordingly. Finally, these single contact models are embedded into the framework of the RN in order to compute effective transport properties. All the proposed models—single contact as well es effective transport models—are verified using finite element methods.

Regarding the pore phase, a novel method concerning the computation of the effective transport properties is developed. By means of the so-called Laguerre tessellation the pore phase of the system is decomposed into cells with each cell surrounding a particle. Consequently, the cell nodes and edges form the basis of equivalent resistor networks. The nodes are identified as pore centers and the edges are the pore throats. As an extension, this model is developed further to account for more than one conducting species in the pore phase. Both methods are either verified using the finite element method or validated with the help of experiments taken from literature.

It is demonstrated that the efficiency of the RN can be used to generate a large database with varying structure combinations. This way, the foundation is created for deriving prediction formulas for the effective conductivity of porous cathode structures represented by sphere packings with overlapping particles as well as the effective resistance of porous secondary particles.

On the other hand, a mathematical model for half-cells with hierarchically structured cathodes with porous secondary particles is proposed. First, the classical half-cell model originating back to Newman and coworkers is revisited and the basic assumptions for the electrochemically based equations are presented. As a next step, the mathematical framework of the volume averaging method is employed to consistently extend the classical to the hierarchically structured half-cell model. For both models, the full set of boundary conditions for the half-cell setup is presented. Finally, the hierarchically structured half-cell model is qualitatively validated by experiments taken from literature.

The validation of this model allows for large-scale parameter studies by varying electronic conductivity and diffusion coefficient of the active material as well as the morphology of the secondary particles. Exemplary results suggest that, while the rate-limiting factor for the classical cathodes is the diffusion coefficient of the active material, in case of the hierarchically structured cathodes, it is the combination of electronic conductivity and inner morphology of the secondary particles.

# **Kurzfassung**

Hierarchisch strukturierte Aktivmaterialien in Elektroden von Lithium-Ionen-Batterien (LIB) mit porösen Sekundärpartikeln sind vielversprechende Kandidaten zur Erhöhung der gravimetrischen Energiedichte und der Ratenabhängigkeit der Zelle. Es gibt jedoch immer noch Aspekte dieser Technologie, die noch nicht vollständig verstanden sind. Um ein tieferes Verständnis darüber zu erlangen, wie die Kathodenstruktur und -morphologie die Zellperformanz beeinflusst, ist das Ziel dieser Arbeit die Entwicklung effizienter Werkzeuge zur Berechnung effektiver Transporteigenschaften für granulare Kathodenstrukturen, die in Zellmodelle importiert werden können, um die elektrochemische Zellleistung von LIBs zu bewerten.

Auf der einen Seite wird die Widerstandsnetzwerkmethode (RN)—ein Werkzeug zur effizienten Berechnung der effektiven Transporteigenschaften von Partikelsystemen—hinsichtlich des Transports durch die Fest- und die Porenphase von granularen Medien, die durch Kugelpackungen mit polydisperser Größenverteilung dargestellt werden, erweitert. Was die Festphase anbelangt, so wird der Transport durch das Volumen der Partikel, über deren Oberfläche oder ein Mix aus beiden betrachtet. Für alle Fälle werden geeignete analytisch hergeleitete Formeln aus der Literatur verwendet oder entsprechend kombiniert, sodass der Widerstand zwischen zwei Einzelpartikeln beschrieben wird. Schließlich werden diese Einzelkontaktmodelle im Rahmen der RN verwendet, um effektive Transporteigenschaften zu berechnen. Alle vorgeschlagenen Modelle—sowohl die Einzelkontakt- als auch die effektiven Transportmodelle—werden mit Finite-Elemente-Methoden (FEM) verifiziert.

Im Hinblick auf die Porenphase wird eine neuartige Methode zur Berechnung der effektiven Transporteigenschaften entwickelt. Mit Hilfe der so genannten Laguerre-Tessellation wird die Porenphase des Systems in Zellen zerlegt, wobei jedes Partikel in ihnen liegt. Die Zellknoten und -kanten bilden die Grundlage für äquivalente Widerstandsnetzwerke. Die Knoten werden als Porenzentren betrachtet und die Kanten sind die Porenhälse. Als Erweiterung wird dieses Modell dahingehend modifiziert, sodass es möglich ist, mehr als eine leitende Spezies in der Porenphase zu berücksichtigen. Beide Methoden werden entweder mit der FEM verifiziert oder mit Hilfe von Experimenten aus der Literatur validiert.

Es wird gezeigt, dass die Effizienz des RN genutzt werden kann, um eine große Datenbank mit unterschiedlichen Strukturkombinationen zu erzeugen. Auf diese Weise lassen sich erfolgreich Vorhersageformeln für den effektiven Widerstand von porösen Sekundärpartikeln sowie die effektive Leitfähigkeit von Kugelpackungen mit überlappenden Partikeln ableiten. Analog zur bekannten Bruggeman-Beziehung können diese Formeln in Zellmodellen verwendet werden, um den Einfluss der effektiven Transporteigenschaften auf die elektrochemische Performanz zu untersuchen.

Auf der anderen Seite wird ein mathematisches Modell für Halbzellen mit hierarchisch strukturierten Kathoden vorgeschlagen. Zunächst wird das klassische, auf Newman und Mitarbeiter zurückgehende Halbzellenmodell rekapituliert und die Grundannahmen für die elektrochemisch basierten Gleichungen vorgestellt. In einem nächsten Schritt wird das mathematische Gerüst der Volumenmittelungsmethode verwendet, um das klassische auf das hierarchisch strukturierte Halbzellenmodell konsequent zu erweitern. Für beide Modelle wird der vollständige Satz von Randbedingungen für den Halbzellenaufbau vorgestellt. Schliesslich wird das hierarchisch strukturierte Halbzellenmodell durch Experimente aus der Literatur qualitativ validiert.

Die Validierung dieses Modells ermöglicht groß angelegte Parameterstudien durch Variation der elektronischen Leitfähigkeit und des Diffusionskoeffizienten des aktiven Materials sowie der Morphologie der Sekundärpartikel. Vorläufige Ergebnisse deuten darauf hin, dass, während der ratenbegrenzende Faktor bei den klassischen Kathoden der Diffusionskoeffizient des aktiven Materials ist, im Falle der hierarchisch strukturierten Kathoden es die Kombination aus elektronischer Leitfähigkeit und innerer Morphologie der Sekundärpartikeln ist.

# **Danksagung**

Die vorliegende Dissertation entstand während meiner Zeit als Doktorand und wissenschaftlicher Mitarbeiter am Institut für Angewandte Materialien - Werkstoff- und Biomechanik (IAM-WBM) des Karlsruher Instituts für Technologie (KIT). Hier durfte ich im Rahmen des Verbundprojekts Werkstoffentwicklung hierarchisch strukturierter Kompositmaterialien für elektrochemische Energiespeicher (HiKoMat) mitwirken.

Das Ziel des Projekts bestand darin, sowohl Partikel- als auch Elektrodenstrukturen bezüglich der elektrochemischen Performanz in Lithium-Ionen-Batterien zu optimieren. Innerhalb eines interdisziplinären Konsortiums gelang es eine Methodenkette zu implementieren, mit der es möglich war, mikrostrukturelle Zusammenhänge zur Zellperformanz in hierarchisch strukturierten Kathoden zu identifizieren.

Ich freue mich, dass ich hierfür mit den in dieser Dissertation entwickelten numerischen Methoden zum Erkenntnisgewinn beitragen und somit eine Hilfe für den Gesamterfolg des Projekts sein konnte. Im Folgenden ist es mir ein Anliegen, den Menschen zu danken, ohne die dies alles nicht möglich gewesen wäre.

Ein besonderer Dank gilt meinem Doktorvater Herrn Prof. Dr. Marc Kamlah, der mir nicht nur wissenschaftlich, sondern auch menschlich stets ein Vorbild war. Den grandiosen Kolleginnen und Kollegen am IAM-WBM möchte ich für die vielen fachlichen und persönlichen Diskussionen und Gespräche auch abseits der Kaffeemaschine danken. Den Projektpartnerinnen und Projektpartnern im HiKoMat-Projekt und allen voran Herrn Dr. Joachim Binder danke ich für die tolle Zusammenarbeit, die nicht zuletzt zum Gelingen der vorliegenden Dissertation beitrug. Überdies möchte ich Herrn Prof. Dr. Michael J. Hoffmann vom Institut für Angewandte Materialien - Keramische Werkstoffe und Technologien (IAM-KWT) für die Übernahme des Korreferats danken.

Abschließend möchte ich meinen Eltern von Herzen danken. Ihr habt mich in all dieser Zeit—in Höhen und Tiefen—liebevoll, geduldig und unterstützend begleitet. Euch widme ich die vorliegende Arbeit.

Karlsruhe, im Juni 2021 *Oleg Birkholz*

# **Contents**




# **1 Introduction**

## **1.1 Lithium-ion batteries**

The continuous development of numerous mobile devices like smartphones, laptops and tablets have a great impact on our daily life. For instance, almost half of the people in Germany own both a smartphone and a tablet [1]. We communicate with each other, schedule appointments, carry out banking transactions or even track our health. Overall, those little helpers have become an important part of our daily routine. Furthermore, fossil fuel independent and electrically driven automobiles continue to enter the market [2]. It can be expected that in some way or the other our mobility will also change. For all these examples lithium-ion battery (LIB) technology plays an important role. It is worth noting that the development and implementation of this technology has come a long way.

### **Development of lithium-ion batteries**

Just recently, in the year 2019, the nobel prize in chemistry was awarded to John Goodenough, M. Stanley Whittingham and Akira Yoshino [3]. The researchers were honored for their contribution to the development and improvement of the lithium-ion battery. The committee recognized that "[t]heir research not only allowed for the commercial-scale manufacture of lithium-ion batteries, but it also has supercharged research into all sorts of new technology, including wind and solar power" [3]. Indeed, over the last three decades the commercial and academic world experienced a large progress and a considerable amount of publications in that field [4].

Whittingham was the first who developed a battery using a titanium disulfide cathode (LiTiS2), a nonaqueous electrolyte and a lithium metal anode. As a result, the cell potential achieved was 2.5V [5]. This was an improvement compared to the lead-acid batteries with a cell potential of 2.0V. Later, Goodenough developed a cathode with an even higher potential. He used lithium cobalt oxide (LiCoO2) where the cell voltage increased up to 4.0V [6]. However, both battery compositions lacked stability due to dendrite growth of the lithium metal anode, which shortens battery lifetime and poses even security risks [7]. Finally, Yoshino used the battery setup from Goodenough, but switched the anode to a carbonaceous material. After running some successful safety tests, in words of Yoshino, this was the birth of the lithium-ion battery [8]. Indeed, modifications of carbonaceous materials were developed ever since [9] and the vast majority of anode material until today is graphite based [10].

### **Design and working principle**

There are two kinds of batteries. Primary batteries can only be discharged once while the secondary batteries can be recharged. Furthermore, a battery is usually built up of so-called cells, where they can be put in series or parallel connection to achieve a desired electric current or voltage. The design and working principle of a secondary battery cell can be observed in Figure 1.1.

The cell has a positive and a negative porous electrode, where both of which consist of so-called active material. This granular material works as a host structure being able to take up or release lithium. Moreover, typically, electronically conductive carbon-black-binder mixture is added to the positive electrode to enhance electronic conductivity. The electrodes are separated by an electronically blocking porous layer, i.e. the separator. The pores of the electrodes as well as the separator are filled with liquid electrolyte to guarantee transport of the lithium ions through the whole cell. A commonly used salt in the electrolyte is lithium hexafluorophosphate (LiPF6), which dissociates into Li<sup>+</sup> and PF<sup>6</sup> <sup>−</sup> in the solvent [10].

Figure 1.1: Design and working principle of a lithium-ion battery cell.

During the discharge process, lithium is oxidized at the active material surface of the anode. Consequently, electrons and lithium ions exit the anode material and move to the cathode via an external circuit and the internal electrolyte phase, respectively. The emerging external electric current can be used to power devices. Together with the oxidation process, lithium ions are reduced at the surface of the cathode material and the lithium atoms enter the crystal structure of the solid. Chemically, the above described process can be expressed as

$$\mathbf{Li} \xleftarrow[\overbrace{\text{reduction}}^{\text{oxidation}} \text{Li}^+ + \mathbf{e}^-,\tag{1.1}$$

which is, in principle, a reversible process. This means that during charging the above described process reverses its direction. Note that the anode and the cathode is defined as the electrode where the oxidation and the reduction occurs, respectively. Technically, only during the discharging the negative electrode is the anode positive electrode is the cathode. Throughout this work, only discharging processes are considered which is why the anode is identical to the negative electrode and the cathode is the positive electrode.

The active material, i.e. the solid particles inside the positive and negative electrode must be capable of taking up and releasing lithium atoms as well as storing them. Commonly, graphite-based materials are used in case of the negative electrodes [10]. Graphite is composed of multiple layers of hexagonal C<sup>6</sup> structures held together by weak van der Waals forces. Lithium atoms can be stored between those layers where a maximum of one lithium atom can be stored per six carbon atoms. Therefore, the notation is Li1−*x*C6, where 0 ≤ *x* ≤ 1 represents the so-called degree of lithiation.

In case of the positive electrode, the used materials are typically oxides of transition metals. The crystal structure must be one that allows lithium to move freely and be stored or released when needed. Additionally, the compound has to be stable for a large range of lithiation. An overview of the possible material structures can be found in [11]. However, new material compositions are continuously under development. Also, commercially used cathode materials can be found in [12]. For example, a commonly used material is lithium manganese oxide (Li*x*(Ni1/3Mn1/3Co1/3)O<sup>2</sup> or NMC), where the capacity is high [13], the rate capability is good [14, 15] and it can operate at high voltages.

A brief look into the future of batteries shows that there is further development in many directions. For instance, in order to tackle stability problems with the liquid electrolyte [16], so-called solid state electrolyte systems are a strong field of research [17]. Moreover, in view of the critical discussions about the resources used in a battery, may it be due to the abundance, geographical distribution or the environmental impact of lithium and other materials, the natural approach is to use sodium-ion batteries instead [18]. However, since those technologies are not market-ready yet, the chosen approach is to tweak the lithium-ion battery technology as far as possible to increase efficiency and thus limit resource consumption.

### **Requirements for lithium-ion batteries**

During the past years, alternative energy sources have been pursued to reduce fossil fuel consumption. Therefore, electrochemical energy conversion and storage is under continuous development [19]. Assuming a renewable energy based electricity production and storage, the promising technologies are batteries, fuel cells and electrochemical capacitors. In this work, the focus is on lithium-ion batteries.

Depending on the field of application the requirements for lithium-ion batteries vary. Generally, the parameters of interest are the specific energy or energy density in the units of W h kg−<sup>1</sup> or W hL−<sup>1</sup> , respectively. The difference between the two is that they are expressed either in terms of mass or in terms of volume. In other words, the best case scenario is having both quantities at high values representing high energy output while keeping weight and space requirements low. Moreover, these quantities can be used to compare different technologies against each other. For instance, the specific energy of gasoline reaches values of 1700W h kg−<sup>1</sup> [20] whereas current lithium-ion batteries yield around 250W h kg−<sup>1</sup> and 440W h kg−<sup>1</sup> against a porous or lithium metal anode, respectively [20].

Apparently, for the lithium-ion battery technology to become a real alternative to, say, fossil fuels, specific energy must approximately double [21]. Since the specific energy is mainly governed by the positive electrode [22], finding better suited positive electrode materials is one solution. Those materials must have greater redox potentials and larger electric capacity [22]. Also, composition and synthesis of the material plays an important role [23, 24]. Additionally, since the materials remain to be poor electronic and ionic conductors, improving the electrode structure can lead to a better exploitation of the applied material. As an example, decreasing the particle size and thus increasing the surface area for the electrochemical reactions leads to higher rate capabilities, which means a more stable discharge capacity for higher currents [25]. As a drawback, however, it was found that cycle stability, i.e. retained capacity over multiple charging and discharging cycles, diminishes. Promising candidates to overcome this problem are the so-called hierarchically structured cathodes where the initial material is processed in such a way that a distinct inner porosity of the active material is achieved. Primary particles, which typically built-up secondary particles of LIB cathodes, now form agglomerates and granular structures inside the secondary particles. Ultimately, in this way created cathodes experience high rate capability and cycle stability at the same time, see [26–30].

## **1.2 State of the art in battery modeling**

There are different methods to model batteries. On the one hand, there are empirical models like equivalent-circuit models. Here, simple systems of cellscale ordinary-difference equations (ODEs) are used, which can be easily and cost-efficiently implemented in battery management systems [31]. On the other hand, there are physics-based models. As the name suggests, they are derived based on laws of physics and they allow predictions whereas the empirical models can only be used within experimentally marked bounds.

Physics-based models can act on different scales [10]. The smallest of which is the molecular scale where the movement of atoms can be tracked and simulated. One scale above is the microscale or the particle-scale where the solid active material particles are spatially resolved, see for example [32–34]. Finally, the macroscale represents the cell level, where the distribution of quantities is not resolved spatially. Instead, the aim is to treat the cell as a continuum while having all inhomogeneities smeared out using appropriate volume averaging methods. This way, the cell can be treated as a homogeneous medium where the microscopic features are represented by effective transport and structure parameters. In this work, the macroscale modeling approach is chosen.

### **Cell modeling**

A mathematical model was developed based on the works by Newman and coworkers [35–37], which is commonly referred to as the Newman cell model. Here, models for electrochemical systems were combined with porous electrodes using the porous electrode theory [38]. In the framework of the porous electrode theory, volume-averaged quantities were used where all the geometrical details are intrinsically accounted for [35]. Ultimately, the electrode is considered to be the superposition of two continua, i.e. the electrolyte and the solid phase. Later, the model was extended to model lithium-ion cells by [36, 39–46]. The key idea of the cell model is that the total electric current inside the cell is the sum of the ionic current, i.e. the current carried by the electrolyte, and the electronic current, i.e. the current carried by the solid phase. Both the currents are linked via electrochemical reactions at the surfaces of active material particles. The reaction is typically described by a charge transfer current or flux using the Butler-Volmer equation [47].

Due to the flexible nature of the cell model, it was extended in many directions [37]. For instance, a model to account for two different particle sizes within one electrode was implemented in [48] and an extension to account for thermal problems was presented in [49]. However, only little effort was taken to model hierarchically structured electrodes [50]. In [51], an impedance model for an agglomerate secondary particle was proposed and in [52], an electrochemicalmechanical model was developed for porous secondary particles. In [53], an agglomerate Newman type cell model was proposed. Here, an additional pore space was introduced on the active material particle level where the transport of species and the electrochemical reactions take place on the internal surfaces. It was possible to simulate better rate capabilities as compared to non-porous secondary particles.

### **Cathode structure modeling**

As mentioned above, apart from material choice, the microstructure of electrodes impacts cell performance [54–60]. Obviously, the particle size distribution influences the microstructure and thus the performance as well [61]. Additionally, the mechanical densification processes, such as calendering of the electrode, changes the microstructure [62–64].

Usually, by densifying the microstructure, the conductivity of the solid phase and therefore the electronic conductivity can be enhanced. However, this also reduces the pore space and thus the ionic conductivity, see [65]. In other words, the effective transport properties, i.e. effective conductivity and effective diffusivity, change with altering microstructure. Ultimately, the composition of the microstructure and the densification process should lead to an optimum balance between the effective transport properties of the solid and the pore phase.

During the past years there have been a couple of theoretically [66–68] and empirically driven prediction formulas [69] to estimate the effective transport properties in porous media. Also, theoretical upper and lower bounds for effective transport properties were derived, as in [70]. Usually, the formulas take the form of functions of the volume fractions of the contributing phases. A wellknown volume fraction based relation to calculate effective transport properties of porous media is the Bruggeman relation [71].

Nowadays, proper tools like the FIB/SEM tomography allow the electrode microstructure to be spatially resolved and virtually reconstructed. Those voxel<sup>1</sup> based data, the effective transport properties can be calculated using finite element methods (FEM) or extended finite element methods (XFEM) [72–74]. This way, the effective transport properties calculated are much more accurate than the above mentioned prediction formulas. However, the tomography, reconstruction and calculation process is expensive in terms of time and resources.

<sup>1</sup> Voxels are 3D extensions of 2D pixels.

When considering the cathode structure as a granular system, a well-established way to model such structures is the Discrete Element Method (DEM) [75]. As compared to spatially resolved FEM or XFEM, in the framework of the DEM, each particle is represented a single element possessing a overall properties only, such as velocity, temperature, potential, and so on. In the most simple case, the elements are spheres, which are geometrically characterized by their coordinates and radii. Modeling granular systems as particulate assemblies is the basis of the Resistor Network Method (RN) [76–80]. For instance, in [76], the effective thermal conductivity is described by a network of contacting spheres where each contact is weighted by a thermal resistance. By extension, this method can be used for non-spherical assemblies, e.g. ellipsoids [81, 82]. As an applied example, in [83], resistor networks were used for the calculation of effective transport properties in solid oxide fuel cells (SOFC). Furthermore, in [84, 85] the RN was employed for the calculation of the solid phase conductivity of both SOFCs and LIBs.

Resistor networks can also be used to calculate the effective transport properties of the pore phase [86]. Among the first was [87–89], where the permeability was calculated using Darcy's law. Pore networks were created based on assemblies of spheres using a so-called Delaunay tessellation. The pore connecting throats were each weighted by a hydraulic conductivity. The hydraulic conductivity was determined based on the geometry based effective throat radius. A similar discretization approach was used in [90]. The so-called Laguerre or generalized Voronoi tessellation was used in [91, 92]. By means of this spatial decomposition technique the pore phase was divided into cells. The cell vertices and edges represent the pore centers and pore throats, respectively, and formed the basis of the resistor network. Pore throat resistances were attributed to the edges based on geometric bottlenecks due to the surrounding particles.

A statistics approach was chosen in [93, 94], where effective transport properties were calculated using a prediction formula based on volume fraction, tortuosity and constrictivity. The latter of which is a quantity resembling the influence of bottlenecks inside the system. To this end, numerous microstructures were generated virtually and the effective transport properties were calculated by simulations. The generation of virtual but realistic microstructure was also done in [95, 96]. Also, recent developments use machine learning techniques on network-based modeling of the effective transport properties of particulate media [97, 98].

## **1.3 Objectives of this work**

The goal of this work is to combine effective transport property delivering tools with well-established cell models. Therefore, the resistor network method is extended in such a way that it becomes capable of providing all necessary quantities needed for cell modeling. Moreover, the Newman cell model is revisited and further developed in order to model hierarchically structured lithium-ion battery electrodes.

The structure of this work is as follows. In Chapter 2, the fundamental theories and tools which are crucial for the whole subsequent work are presented. First, the mathematical equivalence of different transport phenomena is discussed. Second, the volume average method is presented, which is key to the development of the cell models, may it be classical or hierarchically structured cell models. Finally, the theoretical background of the resistor network method is presented and an algorithmic solving scheme is proposed.

Chapter 3 introduces the resistor network method for both the solid and the pore phase. As for the solid phase, the method is extended in particular such that either the transport can be through the volume or via the surface of the particles. Additionally, a combination of both transport mechanisms is presented. Concerning the pore phase, the RN is applied to either the volume of the pore phase or to a mixture of several different species inside the pore phase. Finally, all the extensions are either verified using finite element methods where possible or validated using experimental evidence from literature.

Chapter 4 is dedicated to cell modeling. The derivation of the hierarchically structured cell model necessitates a deeper knowledge of the classical cell model. Therefore, the volume average approach is used to comprehend the standard cell model according to Newman. Later, the same approach is used to consistently derive the hierarchically structured cell model. Additionally, the boundary conditions for both models are mathematically derived and physically motivated.

In Chapter 5, both the resistor network method and the cell models are applied to specific problems. First, structural influences of secondary particles - which are composed of smaller primary particles - on the effective transport is studied by means of the RN. Second, the RN is further used to investigate the structural influence of densely packed sphere assemblies on effective transport via the particle volumes and the surfaces. Finally, the classical and hierarchically structured cell model is applied to real-world cathode structures. Chapter 6 summarizes and concludes this work.

# **2 Fundamentals**

# **2.1 Mathematical equivalence of transport problems**

Under steady-state conditions, the local conservation law can generally be expressed by the continuity equation

$$
\nabla \cdot \vec{F} = 0 \quad , \tag{2.1}
$$

where ~*F* is the flux vector and ∇·(...) is the divergence operator. In case of linear transport, it becomes obvious that different transport phenomena can be treated using the identical mathematical framework, see also [68, 99]. For instance, energy, species and charge conservation is described by

$$
\nabla \cdot \vec{\mathbf{q}} = 0, \quad \nabla \cdot \vec{j} = 0 \quad \text{and} \quad \nabla \cdot \vec{\mathbf{i}} = 0,\tag{2.2}
$$

where ~q, ~*j* and~i is the heat flux density, mass flux density and electric current density vector, respectively.

The flux vectors in Equation (2.2) can be expressed via the constitutive laws dependant on the transport problem under consideration. The thermal, species and charge transport is typically defined by Fourier's, Fick's and Ohm's law

$$\vec{\mathbf{q}} = -\vec{\lambda}\nabla T, \quad \vec{\tilde{j}} = -D\nabla c \quad \text{and} \quad \vec{\tilde{i}} = -\mathbf{x}\nabla\varphi \quad . \tag{2.3}$$

It can be seen here that, in general, the flux is linear dependant on some driving force [100]. Here, λ is the thermal conductivity coefficient, *D* is the diffusion, or diffusivity, coefficient and κ is the electric conductivity coefficient. Furthermore, in the respective constitutive law, the driving force is either the negative gradient of temperature −∇*T* or the negative gradient of concentration −∇*c* or the negative gradient of electric potential −∇ϕ.

## **2.2 Volume averaging method**

In this work, the focus is on porous and heterogeneous materials. In particular, the heterogeneities of the material are considered to be small and homogeneously distributed - in a statistical sense - with respect to the overall dimensions of the system. It is common to refer to the smaller scale as the microscale and the larger scale as the macroscale. Obviously, the structure and the transport processes of the microscale influence the macroscopical transport properties. However, it can be expected that spatially resolved simulations of transport processes on the small heterogeneities level are computationally expensive, if possible at all. Instead, the idea is to subsume information of the microscale and solve the transport equations on the macroscale directly. Typically, a representative volume element (RVE) is defined where its dimensions are sufficiently large to encompass all microscopic phenomena and sufficiently small compared to the system geometry. Inside the RVEs, microscopic equations are transformed to macroscopic ones by employing volume averaging methods [101–105]. In the following, equations defined on the microscale are referred to as microscale equations whereas equations defined on the macroscale are referred to as macroscale equations. The structure on the left-hand side of Figure 2.1 represents a porous and statistically homogeneous material. On the right-hand side of Figure 2.1, the magnified region stands for a representative volume element. The macro- and the microscale is described by the~*x*- and the ~ξ -coordinate vector, respectively, where the representative volume element is defined in terms of ~ξ . Inside the representative volume element, the total volume *V* is occupied by the α- and β-phase such that *V* = *V* <sup>α</sup> +*V* β .

Figure 2.1: Volume average sketch.

Consider a general transport problem in the α-phase described by the continuity equation as

$$\frac{\partial p\_a}{\partial t} + \nabla\_{\vec{\xi}} \cdot \vec{F}\_a = b\_a \,, \tag{2.4}$$

where *<sup>p</sup>*<sup>α</sup> <sup>≡</sup> *<sup>p</sup>*α(~*<sup>x</sup>* <sup>+</sup>~<sup>ξ</sup> ,*t*) is the potential, <sup>~</sup>*F*<sup>α</sup> <sup>≡</sup> <sup>~</sup>*F*α(~*<sup>x</sup>* <sup>+</sup>~<sup>ξ</sup> ,*t*) is the flux vector, and *<sup>b</sup>*<sup>α</sup> <sup>≡</sup> *<sup>b</sup>*α(~*x*+~<sup>ξ</sup> ,*t*) is the source term of this particular phase. The continuity equation must hold for every spatial point ~ξ and time *t* of the RVE which, in turn, is located at a spatial point~*x*.

In order to upscale the problem onto the macroscale, depicted by the *x*coordinate, Equation (2.4) is averaged over the volume of the RVE. Practically, it is integrated over the volume *V* and divided by itself leading to

$$\frac{1}{V} \int\_{V} \left[ \frac{\partial p\_{\alpha}}{\partial t} + \nabla \cdot \vec{F}\_{\alpha} \right] \mathrm{d}V^{\xi} = \frac{1}{V} \int\_{V} b\_{\alpha} V^{\xi},\tag{2.5}$$

or, alternatively,

$$\frac{1}{V} \int\_{V} \frac{\partial p\_{\alpha}}{\partial t} \, \mathrm{d}V^{\sharp} + \frac{1}{V} \int\_{V} \nabla \cdot \vec{F}\_{\alpha} \, \mathrm{d}V^{\sharp} = \frac{1}{V} \int\_{V} b\_{\alpha} V^{\sharp} \, . \tag{2.6}$$

Considering that the integral can be separated as R *V* (...)d*V* <sup>ξ</sup> = R *<sup>V</sup>*<sup>α</sup> (...)d*V* <sup>ξ</sup> + R *<sup>V</sup>*<sup>β</sup> (...)d*V* ξ , where it is defined that α-phase properties are zero in every other phase, Equation (2.6) becomes

$$\frac{1}{V} \int\_{V\_{a}} \frac{\partial p\_{a}}{\partial t} \, \mathrm{d}V^{\xi} + \frac{1}{V} \int\_{V\_{a}} \nabla \cdot \vec{F}\_{a} \, \mathrm{d}V^{\xi} = \frac{1}{V} \int\_{V\_{a}} b\_{a} V^{\xi} \,, \tag{2.7}$$

where *V*<sup>α</sup> ≡ *V*α(~*x*,*t*) is the volume of the α-phase inside the RVE volume.

Take note that volume averages of a property ψ over phase α are defined as

$$
\langle \Psi \mu\_{\alpha}(\vec{x}, t) \rangle = \frac{1}{V} \int\_{V\_{\alpha}(\vec{x}, t)} \Psi\_{\alpha}(\vec{x} + \vec{\xi}, t) \, \mathrm{d}V^{\xi} \tag{2.8}
$$

and

$$
\overline{\Psi}\_{\alpha}(\vec{\mathbf{x}},t) = \frac{1}{V\_{\alpha}(\vec{\mathbf{x}},t)} \int\_{V\_{\alpha}(\vec{\mathbf{x}},t)} \Psi\_{\alpha}(\vec{\mathbf{x}} + \vec{\tilde{\xi}},t) \, \mathrm{d}V^{\tilde{\xi}},\tag{2.9}
$$

where the former of which is called phase average and the latter of which is called intrinsic phase average [106]. Comparing Equations (2.8) and (2.9) it can be concluded that

$$
\langle \Psi\_{\alpha}(\vec{\mathbf{x}},t) \rangle = \underbrace{\frac{V\_{\alpha}(\vec{\mathbf{x}},t)}{V}}\_{\phi\_{\alpha}(\vec{\mathbf{x}},t)} \overline{\Psi}\_{\alpha}(\vec{\mathbf{x}},t) = \phi\_{\alpha}(\vec{\mathbf{x}},t) \overline{\Psi}\_{\alpha}(\vec{\mathbf{x}},t) \,, \tag{2.10}
$$

where φα(~*x*,*t*) is the volume fraction of the α-phase. Consequently, the phase averaged continuity equation in Equation (2.7) can be rewritten as

$$
\langle \frac{\partial p\_a}{\partial t} \rangle + \langle \nabla \cdot \vec{F}\_a \rangle = \langle b\_a \rangle \,. \tag{2.11}
$$

As a next step, each of the terms in Equation (2.11) is transformed using volume averaging theorems [106, 107]. First, the first term on the left-hand side of Equation (2.11) is transformed into

$$
\langle \frac{\partial p\_{a}}{\partial t} \rangle = \frac{\partial \langle p\_{a} \rangle}{\partial t} - \frac{1}{V} \int\_{A\_{a\beta}} p\_{a} \vec{v}\_{a\beta} \cdot \vec{n}\_{a} \, \mathrm{d}A^{\xi} \,, \tag{2.12}
$$

where *A*αβ ≡ *A*αβ (~*x*,*t*) is the interfacial area between the α- and β-phase, <sup>~</sup>*v*αβ <sup>≡</sup> (~*<sup>x</sup>* <sup>+</sup>~<sup>ξ</sup> ,*t*) and <sup>~</sup>*n*<sup>α</sup> <sup>≡</sup>~*n*α(~*<sup>x</sup>* <sup>+</sup>~<sup>ξ</sup> ,*t*) is the velocity and the normal vector of the interfacial area, respectively. Note that the normal vector is pointing from the β- into the α-phase. The integral term represents the change of phase inside the RVE [106]. Assuming non-moving interfacial area and time-independent volume fraction of the α-phase φα, the velocity vector is equal to zero, such that the integral term of Equation (2.12) vanishes, yielding

$$
\langle \frac{\partial p\_a}{\partial t} \rangle = \frac{\partial \langle p\_a \rangle}{\partial t} \tag{2.13}
$$

and the intrinsic volume average expression of Equation (2.13) is

$$\frac{\partial \langle p\_a \rangle}{\partial t} = \frac{\partial \left( \phi\_a \overline{p}\_a \right)}{\partial t} = \phi\_a \frac{\partial \overline{p}\_a}{\partial t},\tag{2.14}$$

respectively. The intrinsic volume average *p*<sup>α</sup> is also called the macroscopic potential.

Second, the second term of the left-hand side of Equation (2.11) is converted using another volume average theorem, leading to

$$
\langle \nabla \cdot \vec{F}\_{a} \rangle = \nabla \cdot \langle \vec{F}\_{a} \rangle - \frac{1}{V} \int\_{A\_{a\vec{\beta}}} \vec{F}\_{a} \cdot \vec{n}\_{a} \, \mathrm{d}A^{\xi} \,. \tag{2.15}
$$

The integral term represents the volume-averaged flux at the interfacial area between the α- and β-phase. Using the average of all surface fluxes *f* αβ , the integral term reduces to

$$\frac{1}{V} \int\_{A\_{a\beta}} \vec{F}\_{\alpha} \cdot \vec{n}\_{\alpha} \, \mathrm{d}A^{\xi} = \underbrace{\frac{A\_{a\beta}}{V} \, \overline{f}\_{a\beta}}\_{a\_{a\beta}} = a\_{a\beta} \overline{f}\_{a\beta} \,, \tag{2.16}$$

where *a*αβ is the specific surface area between the α- and β-phase. In accordance to Section 2.1, the α-phase flux vector

$$
\vec{F}\_a = -k\_a \nabla p\_a \tag{2.17}
$$

is related to the gradient of the potential *p*<sup>α</sup> and conductivity coefficient of *k*<sup>α</sup> of phase α. Before using another volume average theorem on the second term of the right-hand side of Equation (2.15), the constitutive law from Equation (2.17) is inserted first. Extracting the constant conductivity coefficient yields

$$
\langle -k\_a \nabla p\_a \rangle = -k\_a \langle \nabla p\_a \rangle \,. \tag{2.18}
$$

The volume-averaged gradient term in the above equation is

$$
\langle \nabla p\_{\alpha} \rangle = \nabla \langle p\_{\alpha} \rangle - \frac{1}{V} \int\_{A\_{\alpha\beta}} p\_{\alpha} \vec{n}\_{\alpha} \, \mathbf{d} A^{\xi} \,. \tag{2.19}
$$

Assuming that all particles inside the RVE have smooth surfaces and are symmetric with respect to their centers meaning that every normal vector on the surface has a counter part on the other side of the particle which points in the exact other direction. Also *p*<sup>α</sup> is assumed to be uniformly distributed across the surfaces. Eventually, the integral over the surface, i.e. the integral on the right-hand side of Equation (2.19), disappears. This leads to

$$
\left<\nabla p\_{\alpha}\right> = \nabla \left\,,\tag{2.20}
$$

which, expressing in terms of intrinsic volume average, yields

$$
\nabla \langle p\_a \rangle = \nabla \left( \phi\_a \overline{p}\_a \right) = \phi\_a \nabla \overline{p}\_a \,. \tag{2.21}
$$

In spite of the assumptions made before, the volume fraction φ<sup>α</sup> ≡ φα(~*x*) is still dependent on the macroscopic location ~*x*. However, in Equation (2.21), the volume fraction is taken as constant and, thus, extracted from the gradient operator.

Finally, the right-hand side of Equation (2.11) simply is

$$
\langle b\_{\alpha} \rangle = \phi\_{\alpha} \overline{b}\_{\alpha} \,. \tag{2.22}
$$

Summarizing, the macroscopic continuity equation from Equation (2.11) is expressed as

$$
\phi\_a \frac{\partial \overline{p}\_a}{\partial t} + \nabla \cdot \left( -k\_a \phi\_a \nabla \overline{p}\_a \right) - a\_{a\beta} \overline{f}\_{a\beta} = \phi\_a \overline{b}\_a \,. \tag{2.23}
$$

The term *k*αφ<sup>α</sup> is defined using the effective conductivity *k*α,eff, such that Equation (2.23) becomes

$$
\phi\_a \frac{\partial \overline{p}\_a}{\partial t} + \nabla \cdot \left( -k\_{a, \text{eff}} \nabla \overline{p}\_a \right) - a\_{a\beta} \overline{f}\_{a\beta} = \phi\_a \overline{b}\_a \,. \tag{2.24}
$$

The above procedure for the α-phase is performed analogously with respect to the β-phase, which yields

$$
\phi\_{\beta} \frac{\partial \overline{p}\_{\beta}}{\partial t} + \nabla \cdot \left( -k\_{\beta, \text{eff}} \nabla \overline{p}\_{\beta} \right) + a\_{\alpha \beta} \overline{f}\_{\alpha \beta} = \phi\_{\beta} \overline{b}\_{\beta} \,. \tag{2.25}
$$

In accordance to above, φ<sup>β</sup> is the volume fraction, *p*<sup>β</sup> is the macroscopic potential, *k* <sup>β</sup>,eff is the effective conductivity and *b*<sup>β</sup> is the average source term of the β-phase.

## **2.3 Resistor network method**

In this section, the general idea behind the resistor network method is explained and is partially taken from [91]. The method, as used here, follows [108], where the so-called node potential method is introduced. It is being used to calculate electric circuits.

From Section 2.1, irrespective of the underlying physical mechanisms, it becomes clear that the mathematical problems of thermal, diffusive and charge transport has the exact same structure. In the following, the mathematical formulation of the resistor network method is presented using the example of charge transport. When focusing for the calculation of effective transport properties in porous materials on charge transport, may it be electronic or ionic, it is understood in view of the above discussion that the method presented here will also apply to thermal and diffusion processes which plays an important role in the context of lithium-ion batteries, as well.

### **Mathematical formulation of a resistor network**

Consider the example network of nodes and resistors sketched in Figure 2.2.

Figure 2.2: Exemplary electrical circuit to describe the resistor network method. a) Nodes *N I* , *N J* and the currents *I <sup>I</sup>*,*<sup>J</sup>* on the edges between those nodes. b) Potentials ϕ *I* , ϕ *J* at the nodes and resistances *R I*,*J* <sup>κ</sup> at the edges. c) Using potential drop *U* <sup>0</sup>,<sup>1</sup> = ϕ <sup>0</sup> −ϕ 1 in order to calculate the effective transport properties.

The nodes *N I* , *N J* and the currents *I <sup>I</sup>*,*<sup>J</sup>* between those nodes, respectively, are indicated in Figure 2.2a. A yet unknown effective current *I*eff between the nodes *N* 0 and *N* 1 results from a potential drop between those nodes. Note that *N* 0 and *N* <sup>1</sup> have been chosen such that all the other nodes lie between them. In other words, *N* 0 and *N* 1 represent the boundary nodes, i.e. the current collector nodes. At each node *N I* , Kirchhoff's current law,

$$I^I = \sum\_{J}^{n\_{\text{neigh},I}} I^{I,J} = 0 \quad , \tag{2.26}$$

accounting for the conservation of charge, is combined with Ohm's law,

$$I^{I,J} = \frac{U^{I,J}}{R\_{\kappa}^{I,J}} = \frac{\mathfrak{\mathfrak{g}}^I - \mathfrak{\mathfrak{g}}^J}{R\_{\kappa}^{I,J}} \quad , \tag{2.27}$$

representing the constitutive property of the resistors. Here, *n*neigh,I is the number of neighbors of node *N I* and *U I*,*J* is the voltage between the nodes, represented by a potential drop ϕ *<sup>I</sup>* −ϕ *J* . As a result, the equation

$$I^I = \sum\_{J}^{n\_{\text{neigh}, \text{I}}} I^{I,J} = \sum\_{J}^{n\_{\text{neigh}, \text{I}}} \frac{\mathfrak{g}^I - \mathfrak{g}^I}{R\_{\text{K}}^{I,J}} = \sum\_{J}^{n\_{\text{neigh}, \text{I}}} (\mathfrak{g}^I - \mathfrak{g}^J) G^{I,J} = 0 \tag{2.28}$$

can be formulated for each node. Note that the directions of the currents in Figure 2.2a can be chosen arbitrarily as long as the sign of each current in Equation (2.28) is treated consistently. One choice would be that current *I <sup>I</sup>*,*<sup>J</sup>* pointing away from a node *N <sup>I</sup>* has a negative sign.

According to Figure 2.2b, the unknown potentials ϕ *I* , ϕ *J* correspond to the respective nodes *N I* , *N J* and *R I*,*J* <sup>κ</sup> is the resistance between the nodes. Also, the conductance *G <sup>I</sup>*,*<sup>J</sup>* = 1/*R I*,*J* <sup>κ</sup> is used as the reciprocal of the respective resistance. Now, it is possible to assemble a linear system of equations for the given node resistor network.

In order to solve the system of linear equations, the unknown current *I*eff is replaced by an arbitrarily chosen potential drop *U* <sup>0</sup>,<sup>1</sup> = ϕ <sup>0</sup> − ϕ <sup>1</sup> between the boundaries, see Figure 2.2c. After solving for the unknown potentials ϕ *I* , the unknown effective current *I*eff and therefore the effective resistance *R*κ,eff = *U* <sup>0</sup>,1/*I*eff of the system can be calculated.

### **Solution scheme**

In the following, a general scheme for solving resistor networks according to above is presented. Due to its generality, the method described below is well suited to be used in a computer program, see [108].

1) Create the conductivity matrix G as

$$G\_{lj} = \begin{cases} \sum\_{J}^{n\_{\text{neigh},l}} G^{l,J} & \text{if } \quad i = j\\ -G^{l,J} & \text{otherwise} \end{cases} \tag{2.29}$$

and the current vector~*I* as

$$I\_j = \begin{cases} I\_{\rm eff} & \text{if} \quad j = 0 \\ -I\_{\rm eff} & \text{if} \quad j = 1 \\ 0 & \text{otherwise} \end{cases}, \tag{2.30}$$

where *i*, *j* = 0,1,...,*n*, with *n* being the number of nodes. Furthermore, *n*neigh,I is the number of neighbors of the individual node. As shown in Figure 2.2, *I*eff is the current entering the network in node 0 and leaving the network at node 1. The resulting system of linear equations can be set up as

$$
\underbrace{\begin{bmatrix} G\_{00} & G\_{10} & G\_{20} & \cdots & G\_{n0} \\ G\_{01} & G\_{11} & G\_{21} & \cdots & G\_{n1} \\ G\_{02} & G\_{12} & G\_{22} & \cdots & G\_{n2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \underbrace{G\_{0n} & G\_{1n} & G\_{2n} & \cdots & G\_{nn}}\_{\textbf{G}} \end{bmatrix} \underbrace{\begin{bmatrix} \mathfrak{q}\_{0} \\ \mathfrak{q}\_{1} \\ \mathfrak{q}\_{2} \\ \vdots \\ \mathfrak{q}\_{n} \\ \end{bmatrix}}\_{\vec{\Phi}} = \underbrace{\begin{bmatrix} -I\_{\text{eff}} \\ \\ I\_{\text{eff}} \\ \vdots \\ \\ \vdots \\ 0 \\ \end{bmatrix}}\_{\vec{I}}.\tag{2.31}$$

2) Next, the unknown current *I*eff is eliminated from the right-hand side of Equation (2.31). As a result, a modified system of linear equations is created which is smaller by one degree of freedom. To this end, the equality ϕ <sup>0</sup> = ϕ <sup>1</sup> +*U* 0,1 is used and the first row in Equation (2.31) is rewritten as

$$
\begin{bmatrix} G\_{00} & G\_{10} & G\_{20} & \cdots & G\_{n0} \\ & G\_{01} & G\_{11} & G\_{21} & \cdots & G\_{n1} \\ & G\_{02} & G\_{12} & G\_{22} & \cdots & G\_{n2} \\ & \vdots & \vdots & \vdots & \ddots & \vdots \\ & G\_{0n} & G\_{1n} & G\_{2n} & \cdots & G\_{nn} \end{bmatrix} \begin{bmatrix} \Phi\_{1} + U^{0,1} \\ & \Phi\_{1} \\ & & \Phi\_{2} \\ & & \vdots \\ & & \Phi\_{n} \\ & & & \vdots \\ & & & 0 \end{bmatrix} = \begin{bmatrix} -I\_{\text{eff}} \\ I\_{\text{eff}} \\ 0 \\ \vdots \\ 0 \end{bmatrix} \tag{2.32}
$$

Second, the known voltage *U* 0,1 is transferred to the right-hand side which leads to

$$\begin{bmatrix} G\_{00} & G\_{10} & G\_{20} & \cdots & G\_{n0} \\ & G\_{01} & G\_{11} & G\_{21} & \cdots & G\_{n1} \\ & G\_{02} & G\_{12} & G\_{22} & \cdots & G\_{n2} \\ & & \vdots & \vdots & \ddots & \vdots \\ & & & & & & \\ & G\_{0n} & G\_{1n} & G\_{2n} & \cdots & G\_{nn} \end{bmatrix} \begin{bmatrix} \begin{matrix} \boldsymbol{\Psi}\_{1} \\ \hline \begin{matrix} \boldsymbol{\Psi}\_{1} \\ \hline \begin{matrix} \boldsymbol{\Psi}\_{1} \\ \hline \end{matrix} \\ \begin{matrix} \vdots \\ \hline \end{matrix} \\ \end{matrix} = \begin{bmatrix} -I\_{\text{eff}} - G\_{00}U^{0,1} \\ I\_{\text{eff}} - G\_{01}U^{0,1} \\ \vdots \\ -G\_{00}U^{0,1} \\ \end{bmatrix} . \end{bmatrix} . \tag{2.33}$$

Third, the unknown *I*eff is eliminated from the right-hand side of the second line by adding the first two lines. Further, the redundant first line is deleted such that

$$\underbrace{\begin{bmatrix} G\_{01} + G\_{00} + G\_{11} + G\_{10} & G\_{21} + G\_{20} & \cdots & G\_{n1} + G\_{n0} \\ G\_{02} + G\_{12} & G\_{22} & \cdots & G\_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ G\_{0n} + G\_{1n} & G\_{2n} & \cdots & G\_{nn} \end{bmatrix}}\_{\tilde{\mathbf{G}}} \underbrace{\begin{bmatrix} \Phi\_{1} \\ \Phi\_{2} \\ \vdots \\ \Phi\_{n} \\ \vdots \\ \Phi\_{n} \end{bmatrix}}\_{\tilde{\Phi}}$$

$$= \underbrace{\begin{bmatrix} -(G\_{01} + G\_{00})U^{0,1} \\ -G\_{02}U^{0,1} \\ \vdots \\ -G\_{0n}U^{0,1} \\ \tilde{I} \end{bmatrix}}\_{\tilde{I}} \underbrace{\begin{bmatrix} \vdots \\ \vdots \\ \end{bmatrix}}\_{\left(2.34\right)} \end{pmatrix}\_{\left(2.34\right)} \underbrace{\begin{bmatrix} 2.34 \end{bmatrix}}\_{\left(2.35\right)} \underbrace{\begin{bmatrix} \Phi\_{1} \\ \Phi\_{2} \\ \vdots \\ \Phi\_{n} \end{bmatrix}}\_{\left(2.36\right)} \underbrace{\begin{bmatrix} \Phi\_{1} \\ \Phi\_{2} \\ \vdots \\ \Phi\_{n} \end{bmatrix}}\_{\left(2.37\right)} \underbrace{\begin{bmatrix} 2.36 \end{bmatrix}}\_{\left(2.37\right)} \end{bmatrix} \tag{2.34}$$

becomes the modified system of linear equations, where Gˆ is the modified conductivity matrix, <sup>ˆ</sup>~<sup>ϕ</sup> is the modified potential vector and <sup>ˆ</sup>~*<sup>I</sup>* is the modified current vector.

3) Now,

$$
\hat{\mathbf{G}}\,\hat{\vec{\Phi}} = \hat{\vec{I}}\tag{2.35}
$$

can be solved for the modified unknown potential vector ˆ~ϕ.

4) In the next step, the potential values have to be adjusted to the given boundary conditions, that is ϕ<sup>0</sup> = 0 and ϕ<sup>1</sup> = ϕ<sup>0</sup> −*U* 0,1 , by adding the offset potential ϕoff = ϕ<sup>1</sup> −ϕˆ<sup>1</sup> to all ϕ*j*><sup>0</sup>

$$\mathfrak{\mathfrak{q}}\_{j} = \begin{cases} \mathfrak{q}\_{0} & \text{if} \quad j=0\\ \mathfrak{\mathfrak{q}}\_{j} + \mathfrak{q}\_{\text{off}} & \text{otherwise} \end{cases},\tag{2.36}$$

for *j* = 0,...,*n*nodes and reassign potential values ϕ*<sup>j</sup>* to ϕ *<sup>J</sup>* where the value of the index *j* corresponds to node number *J*.

5) In the last step, the unknown current *I*eff can be calculated as the sum of all *n*neigh,I currents entering current collector node *N* 0

$$I\_{\rm eff} = \sum\_{J}^{n\_{\rm nigh}} \frac{\boldsymbol{\varphi}^{0} - \boldsymbol{\varphi}^{J}}{R\_{\rm K}^{0,J}} \quad . \tag{2.37}$$

When the effective current *I*eff is known, the effective resistance *R*κ,eff and therefore the effective conductance *G*eff can be calculated by

$$R\_{\kappa, \text{eff}} = \frac{U^{0,1}}{I\_{\text{eff}}} \quad \text{or} \quad G\_{\text{eff}} = \frac{1}{R\_{\kappa, \text{eff}}} \quad . \tag{2.38}$$

Finally, as in this work the effective conductivity of a representative volume element is considered, domain dimensions have to be taken into account. Usually, rectangular domains are considered of a cross section *A*domain and a length *L*domain such that the effective conductivity can be calculated as

$$\kappa\_{\rm eff} = G\_{\rm eff} \frac{L\_{\rm domain}}{A\_{\rm domain}} = \frac{I\_{\rm eff}}{U^{0,1}} \frac{L\_{\rm domain}}{A\_{\rm domain}} \quad . \tag{2.39}$$

# **3 Modeling effective transport properties in assemblies of spheres**

A large focus of this work is the computation of transport properties in granular materials where the transport can be carried either by solid particles or pores. With modeling those materials using assemblies of spheres, the spheres represent the solid phase and the void space is considered to be the pore phase. In the following, based on the resistor network method scheme, as presented in Section 2.3, models are developed for the computation of effective transport properties in sphere packings. Since the transport can be through the solid phase or the pore phase, models are proposed to create equivalent resistor networks being applicable for both phases, accordingly. Moreover, since this is key to the resistor network method, appropriate models are presented to compute individual resistances inside the network.

# **3.1 Effective transport properties of the solid phase**

To start with, the computation of effective transport properties of the solid phase is considered. Here, the transport can be through the volume or via the surface of the particles. Moreover, a combination of both is possible. Therefore, as a first step, transport through the volume is presented. As a second step, the transport via the surfaces is modeled. Finally, a model is proposed representing a combination of both transport types. In general, the following methods use the geometrical bottleneck effect of contacting particles in order to calculate effective transport properties.

# **3.1.1 Volume resistance of two overlapping solid spheres**

The first transport type considered is via the volume of two overlapping spheres. Therefore, a theory and a model of transport between two overlapping spheres is presented. Afterwards, this model is checked for validity using the finite element method [91].

### **Theory**

Many different cases of heat transfer problems were derived in [109]. In particular, as sketched in Figure 3.1a, a steady-state heat flux φ~ <sup>q</sup> is assumed through a circular aperture between two semi-infinite media. In large distance from the hole, a temperature *T*<sup>0</sup> on the one side and *T*<sup>0</sup> +∆*T* on the other side is applied. In such cases, the thermal resistance can be calculated as

$$R\_{\lambda} = \frac{{1/\lambda^a + 1/\lambda^b}}{4|r\_c|}. \tag{3.1}$$

The thermal bulk conductivities of the semi-infinite media are λ *a* and λ *b* and *rc* is the radius of the circular aperture where the steady-state heat flux φ~ <sup>q</sup> flows through.

If, as a special case, thermal conductivity is equal to λ *<sup>a</sup>* = λ *<sup>b</sup>* = λ *p* , Equation (3.1) reduces to the solution

$$R\_{\lambda} = \frac{1}{2r\_p \lambda^p},\tag{3.2}$$

which is given in [110].

### **Model**

As a next step, Equation (3.1) is employed to calculate the resistance of two overlapping spheres. Therefore, in accordance with Section 2.1, the solution of the thermal problem is applied to a general transport problem. In Figure 3.1b, the two semi-infinite media are interpreted as two half-spheres with the radii *r I* and *r J* , respectively. The bulk conductivities of the *I*'th and *J*'th particle are *k I* and *k J* . Furthermore, the two half-spheres are geometrically overlapping and thus forming a contact radius *r I*,*J <sup>c</sup>* between them. Similar to above, a potential *p*<sup>0</sup> is imposed on the middle-surface of one of the half-spheres and *p*<sup>0</sup> + ∆*p* on the other. The resulting flux is *F*res.

Figure 3.1: Solid-volume resistance. a) Heat flow through a circular aperture between two semiinfinite materials. b) Flow through two overlapping spheres.

Assuming that the potentials are far away from the contact radius and the volume of the half-spheres is comparably large, the analytical Formula in Equation (3.1) can be applied to two overlapping half-spheres, yielding

$$R\_{\text{solid,vol}}^{I,J} = \frac{\mathbf{1}^{\{\mu\}} + \mathbf{1}^{\{\mu\}}}{4 \, r\_c^{I,J}} \,. \tag{3.3}$$

The subscript "solid,vol" refers to the transport through the whole volume of the contacting solid spheres. The herby described resistance shall be called solidvolume resistance which is based on the sphere model.

### **Verification and discussion**

In order to verify Equation (3.3) and thus the solid-volume resistance, finite element simulations were conducted using the commercially available software Abaqus [111]. The FEM provides the spatially resolved solution of the stationary boundary value problem for two half-spheres in contact as shown in Figure 3.1b.

The half-spheres are overlapping to form the contact radius *r I*,*J <sup>c</sup>* . A potential gradient, i.e. temperature gradient ∆*T* in the FEM models, is imposed on the middle surfaces of the spheres. In order to calculate the resulting resistance *R* FEM in the FEM calculation, the total flux, i.e. φ FEM q , at the middle surface of one of the spheres is used to finally obtain

$$R^{\rm FEM} = \frac{\Delta T}{\Phi\_{\rm q}^{\rm FEM}} \,. \tag{3.4}$$

A series of finite element simulations were carried out, where not only the radius ratio *r I* /*r <sup>J</sup>* and the contact radius relative to the smaller sphere *r I* /*r I*,*J <sup>c</sup>* was varied but also the bulk conductivity quotient *k I* /*k J* . In Figure 3.2, the dimensionless representation of the resistance *R*ˆ = *R r<sup>I</sup>* 1/*k I*+1/*k J* is plotted versus the dimensionless contact radius *r I* /*r I*,*J <sup>c</sup>* for three arbitrarily chosen cases but with relatively extreme ratios.

Figure 3.2: Verification of the solid-volume resistance model in Equation (3.3) with varying radii and conductivities using the finite element method.

It can be seen that the values calculated according to Equation (3.3) have a good performance as they agree very well with the results by the finite element simulations. The mean error ¯*e* remains within 4% for relatively large overlaps, e.g. *r I* /*r I*,*J <sup>c</sup>* = 1, as well as for relatively small overlaps, e.g. *r I* /*r I*,*J <sup>c</sup>* = 20. The latter of which means that the contact radius is merely 5% of the smaller of the particle radius.

Furthermore, it has to be noted that the solid-volume resistance from Equation (3.3) even fits well with the FEM reference results when the radius ratio *r I* /*r <sup>J</sup>* was varied from 1/1 to 1/100 and the conductivity ratio *k I* /*k <sup>J</sup>* was varied from 1/1 to 1/300 and 300/1.

In view of the above findings, it can be concluded that Equation (3.3), and thus the solid-volume resistance, is successfully verified to calculate the resistance of two overlapping spheres. It should be noted that the model makes it possible to not only calculate the resistance of two different-sized spheres, but also the bulk conductivities of the particles don't necessarily have to be equal anymore. Therefore, it is an improvement of the formulas proposed in, for example, [76, 84].

# **3.1.2 Effective transport properties of overlapping solid sphere assemblies**

In the above section, it was shown that the volume resistance of two overlapping solid spheres can be described by Equation (3.3), i.e. the solid-volume resistance. For the application to an assembly of spheres, however, it is necessary that the formula also works well inside a network. Therefore, test cases have been produced and the results by RN, where the solid-volume resistance is of crucial importance, are compared to those by FEM analysis. Since there is no analytical solution of the transport problem of random spherical particle assemblies, the results by FEM are taken as the exact solution and thus as the reference [91].

### **Resistor network method for the solid phase**

The following method considers the solid phase of a granular structure to be an assembly of spheres. In contrast to, say, the finite element method, the properties of an individual sphere are not spatially resolved. Instead, each sphere carries overall properties, such as uniform temperature, electric potential, concentration and alike. Overlapping spheres form transport pathways through the assembly and, therefore, increase effective thermal and electric conductivity, diffusivity, etc. of the system.

In order to calculate effective transport properties, the above mentioned transport - or percolated - pathways have to be converted into a network of nodes and resistors. As a first step, as it is highlighted in Figure 3.3a, those clusters of spheres have to be identified which connect the boundaries of the opposing sides. On the topic of cluster identification, refer, for example, to [112–114]. The greycolored particles do not take part in the transport process and can be neglected.

Figure 3.3: Equivalent resistor network to calculate the effective conductivity of the solid volume phase. a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the centers of the particles and edges, i.e. resistances, at the contacts, i.e. overlaps. c) Connection to the current collectors and boundary nodes.

Second, as shown Figure 3.3b, the percolated clusters are converted into equivalent networks by assigning nodes and potentials ϕ *I* , ϕ *J* to the centers of the particles. Also, resistors *R I*,*J* solid,vol are attributed to the edges between those nodes according to the geometric relation from Equation (3.3). Finally, as shown in Figure 3.3c, additional nodes are added to model the boundary nodes where the boundary conditions are imposed on. The resistor network is thus created which is the basis for the calculation of the effective conductivity according to the scheme presented in Section 2.3.

### **Model description**

As mentioned previously, the transport problems in this work, can be described using equivalent mathematical forms, see Section 2.1. Therefore, the methods used for, say, heat transport problems can be used equivalently for electric transport problems. Accordingly, the FEM solution is acquired by using the heat transport module in Abaqus [111], where the 10-node quadratic tetrahedron element DC3D10 is chosen. Similar to above, the finite element simulation serves as the reference since it is the spatially resolved exact solution of the steady-state boundary value problem. Finally, the effective transport properties, achieved by both the FEM and RN methods, are compared to each other. In case the results by the resistor network method don't deviate from the finite element solution, the former of which can be regarded as verified.

There were two steps involved in generating the particle structures. First, the initial structure was created using the random close packing algorithm (RCP) [115]. Generally, the RCP produces randomly distributed, densely packed and overlap-free assemblies of spheres. Following the approach of [116, 117], the algorithm scope of application was extended to cover any given size-distribution of the particles, i.e. the radii of the spheres. In this case, a normal distribution with a given mean radius *r*mean and standard deviation *r*<sup>σ</sup> .

Secondly, the initial structure was further densified using an algorithm which will be called numerical sintering [84, 85]. While keeping the centers fixed in space the radii of the spheres inside the assembly were successively increased until a ceratin threshold was reached. The threshold can be a targeted packing factor, i.e. volume fraction of the solid phase φsolid, or the mean contact angle θ*c*,mean, see [85]. The contact angle θ *I*,*J <sup>c</sup>* of two overlapping spheres is defined as the bigger of the two angles θ *I c* and θ *J c* enclosing the contact radius *r<sup>c</sup>* , see Figure 3.4.

Figure 3.4: Contact angles of two contacting spheres.

Finally, θ*c*,mean is calculated as the mean of the contact angles of all contact pairs in the system. In the following, the threshold was defined using the mean contact angle. It is important to note that during this densification routine any sorts of mechanics were neglected, i.e. contact forces and alike. The geometrical data of the spheres were then imported into a box-shaped simulation domain and in both, the FEM and RN analysis, transport through two opposite conducting surfaces in z-direction was considered.

Concerning the FEM models, as shown in Figure 3.6b, the spheres were cut off at the conducting surfaces and all nodes on one of these surfaces were set to the same temperature, while the potential of the nodes on the opposite surface were dropped by ∆*T* = 1 with respect to this potential. The resulting heat flux φ FEM q was obtained at one of the surfaces where temperature was applied. By using the domain length *L*domain and cross section area *A*domain, it is possible to calculate the effective thermal conductivity as

$$
\lambda\_{\rm eff}^{\rm FEM} = \frac{\phi\_{\rm q}^{\rm FEM}}{\Delta T} \frac{L\_{\rm domain}}{A\_{\rm domain}} \,. \tag{3.5}
$$

As for the RN model, the exact same assembly of spheres can be seen in Figure 3.6a. In the above described framework of the resistor network method, the sphere assembly has to be converted into an equivalent circuit of nodes and resistors. In order to ensure the same electric potential at the conducting surfaces of the simulation domain, additional boundary nodes were created lying exactly on those surfaces. These boundary nodes were constructed as the contact points of those spheres which were overlapping with the conducting surface. The resistance between a boundary node and a node in the simulation domain is calculated as

$$R\_{\text{solid,vol}}^{I,0} = \frac{1/\kappa^l}{4r\_c^{I,0}}\,,\tag{3.6}$$

where *r I*,0 *<sup>c</sup>* is the contact radius of the sphere *I* and the conducting surface, see Figure 3.5.

Figure 3.5: Construction of boundary nodes and computation of boundary solid-volume resistances.

For the resulting network, an electric potential drop of ∆ϕ = 1 was imposed on the boundary nodes of two opposing surfaces. By considering the domain length and cross section area, the effective electric conductivity can be calculated via

$$\mathbf{k}\_{\rm eff}^{\rm RN} = \frac{I\_{\rm eff}}{\Delta \boldsymbol{\wp}} \frac{L\_{\rm domain}}{A\_{\rm domain}},\tag{3.7}$$

where the effective current *I* eff was calculated using the resistor network method scheme from Section 2.3.

### **Test cases**

For the purpose of verification, up to 15 assemblies were created, following the method described above. The cases were divided into certain types of assemblies, see Table 5.1. From the first to the third assembly type, the mean radius was fixed to an arbitrary length unit of 1 and the mean contact angle was fixed to 15°, respectively. The only difference among them was the standard deviation of particle radii, which ranged from 0 to 0.25, such that with increasing assembly type number the degree polydispersity increased, as well.

Figure 3.6: Verification models of the solid-volume resistor network method for an assembly of spheres. a) Potential distribution through the network of contacting spheres. b) Temperature distribution of the spatially resolved finite element mesh.

In Figure 3.6, an example of the test cases is presented. The shown assembly type 3 is the one with the largest polydispersity.

Table 3.1: Structural parameters of the spherical packings for the RN verification of the transport through the volume of solid spheres.


### **Evaluation and discussion**

In order to compare the effective transport properties, i.e. the thermal and electric conductivity, provided by the FEM and the RN simulations, respectively, a socalled effective transport parameter <sup>ˆ</sup>*k*eff is defined. A dimensionless ratio of the effective transport property divided by its bulk conductivity. In case of the thermal transport problem and the electric transport problem, this yields

$$
\hat{k}\_{\rm eff}^{\rm FEM} = \frac{\lambda\_{\rm eff}^{\rm FEM}}{\lambda\_{\rm bulk}} \quad \text{and} \quad \hat{k}\_{\rm eff}^{\rm RN} = \frac{\kappa\_{\rm eff}^{\rm RN}}{\kappa\_{\rm bulk}}, \tag{3.8}
$$

respectively. λbulk is the thermal and κbulk is the electric bulk conductivity. In Figure 3.7, transport parameters resulting from FEM and RN simulations are plotted against each other. On the horizontal axis the results from the FEM solution are shown and on the vertical axis the corresponding values from the RN analysis are plotted. The black solid line in this figure represents a perfect match of both the results whereas values above or below indicate an over- or underestimation of the effective transport parameter by the RN with respect to FEM.

Figure 3.7: Verification of the solid-volume resistor network method using finite element simulations.

It can be seen that the RN results are overestimating the actual effective transport. However, the average error is around *e* = 5% which is within an acceptable range considering the fact that the very complex structure is merely approximated by nodes and resistors. In other words, the exact solution is achieved with a lot less degrees of freedom and within a much shorter time [91].

From the above observations, a very good agreement can be stated for the resistor network approach with the spatially resolved finite element results. Therefore, the RN is considered to be verified to calculate the effective transport properties via the solid volume of assemblies of overlapping spheres with polydisperse sizedistributions.

# **3.1.3 Surface resistance of two overlapping solid spheres**

It was shown that the volume resistance of two overlapping spheres can be calculated by the geometrical contact radius, on the one hand, and the bulk conductivities of the contacting spheres, on the other hand, see Equation (3.3). For some applications, however, the transport doesn't necessarily have to be via the volume. To increase the electronic conductivity, for example, non-conducting active material particles are coated with a high-conductive material [118]. In those cases it may be justified to assume that the transport is solely over the surface of the particles.

### **Theory**

In [119], the electric current flow over the surface of spherical particles was investigated. The electric surface resistance of two overlapping spheres was found to be

$$R\_{\kappa, \text{surf}} = \frac{\rho\_{\text{surf}}}{2\pi} \ln \left( \frac{\tan(\theta/2)}{\tan(\theta/2)} \right) \,. \tag{3.9}$$

Here, ρsurf is the electric surface resistivity, θ*<sup>t</sup>* and θ*<sup>c</sup>* shall be called transport angle and contact angle, respectively.

In Figure 3.8a, exemplarily, the sketch of two spheres in contact can be seen. It is assumed that the electric current~*I* flows completely via the surface and, in addition, between the contact angle θ*<sup>c</sup>* and the transport angle θ*<sup>t</sup>* . The transport angle θ*<sup>t</sup>* must be in the range of [θ*c*,180°) because otherwise the logarithmic expression becomes either negative or infinite.

### **Model**

In order to model the transport via the surface, a suitable shell model is developed. The surface of a sphere is represented by a thin shell, where the bulk conductivity *k*bulk of the material under consideration is applied to. The surface resistivity from above is related to the bulk conductivity as

$$
\rho\_{\rm surf} = \frac{1}{k\_{\rm bulk} \cdot s},
\tag{3.10}
$$

where *s* is the shell thickness. The idea behind this model is that as long as the shell is thin enough, the surface transport is accounted for.

In Figure 3.8b, the model just described is shown. Here, two half-shells with the radii *r I* and *r J* are geometrically overlapping to form the contact radius *r I*,*J <sup>c</sup>* . Now, due to the fact that the former radii can differ in size, the resulting contact angles θ *I c* and θ *J <sup>c</sup>* may also be different. The transport angles are θ *I*,*J <sup>t</sup>* and θ *J*,*I t* . The respective thickness of the shells is given as *s I* and *s J* . Also, the shell bulk conductivities are *k I* and *k J* .

Figure 3.8: Solid-surface resistance. a) Electric current via the surface of two contacting spheres. b) Flow through two overlapping shells.

Finally, the resulting resistance of two overlapping shells is calculated by a series connection of the surface resistance, as in Equation 3.9, where the shellthicknesses and -resistances are accounted for by using the shell model from Equation 3.10.

$$\begin{split} R\_{\text{solid,surf}}^{I,J} &= R\_{\text{solid,surf}}^{I} + R\_{\text{solid,surf}}^{J} \\ &= \frac{1/k^{I}s^{J}}{2\pi} \ln \left( \frac{\tan(\theta\_{\text{r}}^{I,J}/2)}{\tan(\theta\_{\text{c}}^{I}/2)} \right) + \frac{1/k^{J}s^{J}}{2\pi} \ln \left( \frac{\tan(\theta\_{\text{r}}^{I,I}/2)}{\tan(\theta\_{\text{c}}^{I}/2)} \right) . \end{split} \tag{3.11}$$

The subscript "solid,surf" refers to the transport via the surface of the solid particles. Thus, the resistance described here shall be called solid-surface resistance whereas the underlying model is called the shell model.

### **Verification and discussion**

A series of FEM scenarios were generated to verify Equation (3.11). Following the model in Figure 3.8b, the finite element models were created accordingly. Therefore, geometrically overlapping shells were created. To be precise, similar to Section 3.1.1, half-shells were created. A potential gradient was represented by a temperature drop of ∆*T* = 1 between the nodes at the mid-surfaces of the half-shells. The resulting flux is the total heat flux φ FEM q at one of the top or bottom mid-surface. It is used to calculate the resistance of the FEM models by

$$R^{\rm FEM} = \frac{\Delta T}{\Phi\_{\rm q}^{\rm FEM}}\,. \tag{3.12}$$

In Figure 3.9, the results of the investigation can be seen. Here, the dimensionless representation of the resistance *R*ˆ = *R* 1/*k I s I*+1/*k J s J* is plotted versus the ratio *r I* /*r I*,*J c* , where *r I* is the smaller of the two radii. The ratio ranges from around 1 up to around 10, where the former means that the contact radius is almost the same as the smaller radius of the two and the latter means that it is 10% of the size of the smaller radius.

As an example of normal to extreme cases, the radius ratio *r I* /*r <sup>J</sup>* was varied from equal-sized 1/1 to relatively large ratios of 1/10. Furthermore, the conductivity ratio *k I* /*k <sup>J</sup>* was also varied from 1/1 up to 1/300. Lastly, in Figure 3.9a, b and c, the shell thickness ratio *c* = *s I* /*r <sup>I</sup>*−*s <sup>I</sup>* = *s J* /*r <sup>J</sup>*−*s <sup>J</sup>* was varied from 0.05 via 0.1 to 0.15. This stands for the ratio of the shell thicknesses to the inner radii which, in turn, are the radii of the inner surfaces of the shells. For all cases, θ *I*,*J <sup>t</sup>* and θ *J*,*I <sup>t</sup>* were set to 90°.

Generally, evaluating Figure 3.9, it can be noted that *R I*,*J* solid,surf fits best for relatively low values of *r I* /*r I*,*J <sup>c</sup>* , that is for large overlaps. However, in Figure 3.9a, it can be seen that even for a relatively large *r I* /*r I*,*J <sup>c</sup>* ratio the formula performs best for thin shells with *c* = 0.05, where the mean error is around *e* = 3%. Also, for moderate shell thicknesses of *c* = 0.10, the mean error of *e* = 8% is accaptable.

Figure 3.9: Verification of the solid-surface resistance model in Equation (3.11) with varying radii and conductivities using the finite element method. a) Shell thickness ratio *c* = 0.05. b) Shell thickness ratio *c* = 0.10. c) Shell thickness ratio *c* = 0.15.

From what was found above, the solid-surface resistance and, in turn, the shell model is considered to be successfully verified. Using this model, it is possible to calculate the resistance between two thin overlapping shells, where not only the radii can be different, but also the conductivities of the shells. The only restriction is that the shells have to be thin enough. That is to say, the thickness of the shells should not exceed the radii of the sphere which they are covering by 15%.

# **3.1.4 Effective transport properties of overlapping solid shell assemblies**

It was shown, that the resistance of two single contacting shells can be calculated using the solid-surface resistance. Next, it is investigated if the proposed formula in Equation 3.11 holds for an assembly of shells.

### **Resistor network method for shell assemblies**

Recall that the solid phase of a granular structure was taken as collection of spheres. Since now the focus is on the transport via the surface of the spheres, rather, the structure is considered as an assembly of thin shells. Similar to above, overlapping shells form conducting pathways through the system and contribute to the effective conductivity. In order to calculate the effective conductivity of shell assemblies, the resistor network is used and adapted, accordingly.

In Figure 3.10, the general approach on how to apply the RN scheme on the shell model is sketched. First, conducting, i.e. percolating, pathways through the assemblies have to be found, see the highlighted shells in Figure 3.10a. Next, potentials ϕ *I* and ϕ *J* are assigned to the shell centers, i.e. nodes, and the resistances *R I*,*J* solid,surf are assigned to the contact pairs, i.e. edges, as can be seen in Figure 3.10b. Finally, in Figure 3.10c, the boundary conditions are applied to the boundary nodes and the equivalent node network can be solved using the scheme according to Section 2.3.

Figure 3.10: Equivalent resistor network to calculate the effective conductivity via the solid surface. a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the centers of the shells and edges, i.e. resistances, at the contacts, i.e. overlapping shells. c) Connection to the current collectors and boundary nodes.

While using the solid-surface resistance for the contact pairs in Figure 3.10b, setting the transport angle between two contacting shells in the assembly needs special treatment. That is because the presence of more than one contact partner of a shell is likely to influence the transport path along the shell and decreases the transport angle. In Figure 3.11, exemplarily, four shells are sketched, where the *I*'th shell in the middle is in contact with the surrounding *J*'th, *K*'th and *L*'th shell. Now, in order calculate the resistance *R I*,*J* solid,surf between the *I*'th and the *J*'th shell, the adjacent contact pairs have to be accounted for by adjusting the transport angle θ *I*,*J t* in Equation (3.11).

Consider for the solid-shell resistance *R I*,*J* solid,surf between shell *I* and *J* the contribution of the region around the *J*'th and *L*'th shell. Apparently, the transport path becomes shorter, the closer the *L*'th shell approaches *J*. The length of the transport path is represented by the transport angle θ*<sup>t</sup>* and, in particular, θ *I*,*J*−*L <sup>t</sup>* describes the transport path between the contact pair *I* and *J* and the neighborhood of *L*. Note that the transport angle θ *I*,*J*−*L t* starts from the axis between the centroids of *I*'th and *J*'th shell and ends at the contact radius of the *I*'th and *L*'th shell.

Figure 3.11: Calculation of the transport angle for shells with adjecent contacts.

In 3D, the angle is measured on the plane constructed by the axis between the centroids of *I*'th and *J*'th shell and the axis between the centroids of *I*'th and *L*'th shell. Since there can be more than one contacting neighbor to *I*, the other contact pairs have to be considered as well. Therefore, in a similar fashion, the transport angle θ *I*,*J*−*K t* is calculated.

Finally, the transport angle θ *I*,*J t* , entering Equation (3.11), is taken as the average of all transport angles. Thus,

$$\Theta\_{l}^{I,J} = \frac{1}{n\_{\text{neigh},I}} \sum\_{n=1}^{n\_{\text{neigh}}} \Theta\_{l}^{I,J-l(n)} \quad \text{for} \quad i(n) = J, K, L, \dots \\ i(n\_{\text{neigh},I}) \,, \tag{3.13}$$

where *n*neigh,*<sup>I</sup>* is the number of contact neighbors of shell *I* and *i*(*n*) is the index of the *n*'th shell. Note that it is chosen that θ *I*,*J*−*J <sup>t</sup>* = 90° in Equation (3.13) and the maximim transport angle is set to 90°.

### **Model description**

The generation of the shell assemblies was done in a similar fashion to Section 3.1.2. The RCP is used to generate randomly distributed, densely packed and overlap-free collections of spheres. Additionally, the particle size of the generated structures can be a polydisperse distribution, following a normal distribution. Then, the structures were densified until a certain threshold was reached, in this case it was the mean contact angle.

The densified sphere assemblies were the basis for the FEM models. In order to generate shells, a shell-thickness must be provided. With regard to the shell thickness, for every sphere imported from the generated structures, another sphere around the same center was imported, as well, with the shell thickness substracted from its radius. As a next step, the inner sphere was geometricaly substracted from the outer sphere and hence leaving a shell.

In Figure 3.13b, the meshed shells can be seen. Additionally, with respect to the dimensions of the simulation box, the parts of the shells outside the box in z-direction are cut off. The color gradient in the same figure represents the temperature gradient. For the finite element model, a temperature gradient of ∆*T* = 1 was imposed on the nodes of the opposing surface of the domain in zdirection on the cutting plane.

Finally, the resulting heat flux φ FEM <sup>q</sup> was evaluated at the top or bottom surface in z-direction to achieve the effective thermal conductivity as

$$
\lambda\_{\rm eff}^{\rm FEM} = \frac{\Phi\_{\rm q}^{\rm FEM}}{\Delta T} \frac{L\_{\rm domain}}{A\_{\rm domain}} \,, \tag{3.14}
$$

domain length and domain cross-section area is *L*domain and *A*domain, respectively.

The generated sphere assemblies from above were imported and solved using the RN framework. In accordance to the boundary conditions in the finite element model, an electric potential drop of ∆ϕ = 1 is applied to the boundary nodes in zdirection. The boundary nodes were generated as the contact points of the shells and the boundary plane in z-direction und, thus, lying exactly on that plane. The resistance between a boundary node and a node within the simulation domain was calculated as

$$\mathcal{R}^{I,0}\_{\text{solid,surf}} = \frac{1/k^{I}s^{I}}{2\pi} \ln\left(\frac{\tan(\theta^{I,0}\_{\text{t}}/2)}{\tan(\theta^{I,0}\_{\text{c}}/2)}\right),\tag{3.15}$$

47

where θ *I*,0 *<sup>c</sup>* and θ *I*,0 *t* is the contact and transport angle between boundary plane and shell, respectively, see Figure 3.12. Note that θ *I*,0 *t* is computed as described above, where the axis between the *I*'th shell and the boundary 0 is constructed between the shell center and the boundary node coordinates.

Figure 3.12: Construction of boundary nodes and boundary solid-surface resistances.

In Figure 3.14a, the electric potential distribution of an example assembly can be seen. The simulation box is depicted by the black borders.

The resistor network analysis provided the effective electric conductivity

$$
\kappa\_{\rm eff}^{\rm RN} = \frac{I\_{\rm eff}}{\Delta \varphi} \frac{L\_{\rm domain}}{A\_{\rm domain}} \,, \tag{3.16}
$$

where the effective current *I* eff was calculated using the resistor network solving scheme according to Section 2.3.

#### **Test cases**

In the following, up to 15 test cases were produced and, based on the test cases, the RN results were compared to the ones from FEM simulations. Three types of assemblies were generated. As can be seen in Table 3.3, the assemblies differ in their standard deviation of the radii ranging from 0.0, i.e. mono-sized particles, to 0.20, i.e. polydisperse particle size distribution. In every one of the cases, the shell thickness ratio was set to *c* = *s*/*r*−*s* = 0.1. This corresponds to a moderate shell thickness according to the investigation above, which resembles the average case. In Figure 3.13, the RN and FEM models are compare to each other. Here,

Table 3.3: Structural parameters of the spherical packings for the verification of the RN transport via the surface of solid spheres.


the assembly is of the type 3 whith the largest standard deviation *r*<sup>σ</sup> .

Figure 3.13: Verification models of the solid-surface resistor network method for an assembly of shells. a) Potential distribution through the network of discrete elements. b) Temperature distribution of the spatially resolved finite element mesh.

### **Evaluation and discussion**

For comparison reasons, previously mentioned dimensionless effective transport parameters are used. Therefore, the effective thermal and electric conductivities are normalized as

$$
\hat{k}\_{\rm eff}^{\rm FEM} = \frac{\lambda\_{\rm eff}^{\rm FEM}}{\lambda\_{\rm shell}} \text{and} \quad \hat{k}\_{\rm eff}^{\rm RN} = \frac{\kappa\_{\rm eff}^{\rm RN}}{\kappa\_{\rm shell}}, \tag{3.17}
$$

respectively, where λshell is the thermal and κshell is electric conductivity of the shells.

In Figure 3.14, the results can be observerd. The effective transport parameters are plotted against each other. On the horizontal axis, the results from the FEM solution are potted whereas on the vertical axis, the RN reults are provided. The black solid line indicates a perfect match of both results. Points above or below mean over- or underestimation of the finite element results which are taken as reference results.

Figure 3.14: Verification of the solid-surface resistor network method using finite element simulations.

In general, the results provided by the resistor network method overestimate the reference finite element results. On the other hand, the mean error ¯*e* is below 5%. This is a reasonably low error, given that the solid-surface model is a rather rough approximation of the suface transport over a very complex structure.

# **3.1.5 Mixed resistance of overlapping core-shells**

The shell model, as presented above, can be used when the transport properties of the conductive coating is much larger than the core. Due to the large variety of different coatings for the cathode and anode material [120], for instance, there may be cases where the transport property of the core material can not be neglected [121].

### **Model**

To account for such cases, the idea now is to simply superimpose the solidvolume transport, Figure 3.1, and the solid-surface transport, Figure 3.8. As a result, in Figure 3.15, the volume part is represented by the core and the surface part is represented by the shell.

Both the core and the shell radii are represented as *r I* and *r J* . The cores have a conductivity of *k I* core and *k J* core whereas the shells have a conductivity of *k I* shell and *k J* shell. Furthermore, the two resistances, volume and surface, are sharing the same contact radius *r I*,*J <sup>c</sup>* and contact angles θ *I c* and θ *J c* . The transport angles are θ *I*,*J <sup>t</sup>* and θ *J*,*I t* .

Figure 3.15: Flow through two overlapping core-shells.

As a consequence, the superposition of the two transport phenomena leads to a parallel connection of the volume and the shell resistance model, Equations (3.3) and (3.11), respectively, yielding

$$\begin{aligned} R\_{\text{solid,mix}}^{J,J} &= \left( \frac{1}{R\_{\text{solid,surf}}^{J,J}} + \frac{1}{R\_{\text{solid,vol}}^{J,J}} \right)^{-1} \\ &= \left( \frac{1}{\frac{1/k\_{\text{solid}}^J J}{2\pi} \ln \left( \frac{\tan(\theta\_t^{J,J}/2)}{\tan(\theta\_t^{J,j}/2)} \right) + \frac{1/k\_{\text{solid}}^J J}{2\pi} \ln \left( \frac{\tan(\theta\_t^{J,J}/2)}{\tan(\theta\_t^{J,j}/2)} \right)} + \frac{1}{\frac{1/k\_{\text{gas}}^J + 1/k\_{\text{GC}}^J}{4J\_c^J}} \right)^{-1} . \end{aligned} \tag{3.18}$$

Since the equivalent resistance is a mixed form of solid-volume and solid-surface resistance, it will be called solid-mix resistance and the underlying model is called core-shell model.

### **Verification and discussion**

For verification purposes, and as was done above for the previous cases, a series of geometrically overlapping core-shell FEM scenarios were created. The cores and shells were meshed separately but shared the same nodes at the interface between core and shell. Hence, the temperature in the finite element simulation was continiuous at those interfaces. Therefore, the FEM simulation represents a physically more realistic situation than the superposition of the core-shell model.

Exemplarily, in Figure 3.16, those results are shown where the radius ratios *r I* /*r J* take moderate values, i.e. 1/1 and 1/5, and rather extreme values, i.e. 1/10. Also, the shell conductivities *k I* shell and *k J* shell range both from 2 up to 10 whereas the cores conductivities always remain at 1. Finally, for all combinations, the shell thickness ratio *c* = *s I* /*r <sup>I</sup>*−*s <sup>I</sup>* = *s J* /*r <sup>J</sup>*−*s <sup>J</sup>* was varied from 0.05, via 0.1, to 0.15 to examine the influence of the shell thickness.

The results of the investigation are presented in Figure 3.16. Here, the dimensionless representation of the resistances

$$\hat{R} = R\left(\frac{1}{1/k\_{\rm shell}^{I}s^{I} + 1/k\_{\rm shell}^{J}s^{J}} + \frac{r^{I}}{1/k\_{\rm core}^{I} + 1/k\_{\rm core}^{J}}\right) \tag{3.19}$$

is plotted versus the ratio *r I* /*r I*,*J <sup>c</sup>* . *R* represents either the single points provided by the FEM simulations or the dashed lines calculated by the solid-mix resistance from Equation (3.18). Also, for every combination, the relative mean error ¯*e* is shown.

On the one hand, in Figure 3.16, it can be seen that for all cases the numerical results agree the best with the results by the core-shell model, if the overlap is the largest, i.e. *r I* /*r I*,*J <sup>c</sup>* → 1. On the other hand, even for relatively small overlaps, i.e. *r I* /*r I*,*J <sup>c</sup>* → 10, the best agreement is achieved if the shell thickness takes its lowest ratio 0.05, see Figure 3.16a. For example, in case of the largest radius ratio *r I* /*r <sup>J</sup>* = 1/10 and the largest shell conductivities of 10, the mean relative error is around 1%.

Figure 3.16: Verification of the solid-mix resistance model *R I*,*J* solid,mix as presented in Equation (3.18) with varying radii and conductivities using the finite element method. a) Shell thickness ratio *c* = 0.05. b) Shell thickness ratio *c* = 0.10. c) Shell thickness ratio *c* = 0.15.

54

### **Influence of shell and bulk conductivity ratio**

In the following, the limits of the core-shell model, as it is used in Equation (3.18), are discussed. To this end, the shell to core conductivity ratio *k*core/*k* shell has been varied from <sup>1</sup>/<sup>1</sup> up to rather extreme ratios of <sup>1</sup>/100. While on the other hand, the radius and shell-thickness ratios were kept constant as *r I* /*r <sup>J</sup>* = 1/1 and *c* = *s I* /*r <sup>I</sup>*−*s <sup>I</sup>* = *s I* /*r <sup>J</sup>*−*s <sup>J</sup>* = 0.05.

In Figure 3.17, the results of the above described investigation are shown. On top of that, the theoretical limits of the core-shell model are plotted. These are, on the one hand, the solid-volume resistance *R I*,*J* solid,vol from Equation (3.3), if the volume part becomes dominant and, on the other hand, the solid-surface resistance *R I*,*J* solid,surf from Equation (3.11), if the surface part becomes dominant.

Figure 3.17: Influence of conductivity ratio on the solid-mix resistance model *R I*,*J* solid,mix as presented in Equation (3.18) and comparioson to the solid-volume model *R I*,*J* solid,vol taken from Equation (3.3) as well as the solid-surface model *R I*,*J* solid,surf described by Equation (3.11).

Unsurprisingly, in case of the conductivity ratio being 1, the FEM solution tends to the solid-volume resistance curve. Furthermore, the solid-mix resistance is also approaching the solid-volume curve. However, especially for small overlaps, that is *r I* /*r I*,*J <sup>c</sup>* → 10, the solid-mix resistance deviates from the FEM or the solidvolume resistance by around ¯*e* = 24% with respect to the numerical solution. The error originates from the fact that the core-shell model is a superposition of the sphere and the shell model.

Another interesting observation can be made whith regard to the cases where the shell conductivity is 100 times larger than the core conductivity. It seems that this ratio is sufficient to almost exactly match the solid-surface resistance. Also, the relative mean error of the core-shell model to the FEM model is merely *e*¯ = 4%. Finally, the intermediate states are represented by the conductivity ratio *k*core/*k* shell = <sup>1</sup>/10. The curve lies between the solid-volume resistance and the solidsurface resistance and shows a relative mean error of ¯*e* = 3% with respect to the finite element model.

Thus, it can be concluded that the core-shell model, i.e. the solid-mix resistance in Equation (3.18), is capable of covering the limits of solid-volume and solidsurface resistance.

# **3.1.6 Effective transport properties of overlapping solid core-shell assemblies**

It was shown, that the resistance of two contacting core-shells can be calculated using Equation (3.18). Next, it is investigated if the formula holds for an assembly of core-shell particles. Therefore, a series of assemblies of such particles were generated and the effective transport properties provided by RN are compared to the results by FEM analysis.

### **Resistor network method for core-shell assemblies**

In Figure 3.18, the resistor network approached is sketched for an assembly of core-shell particles. First, the percolating clusters are identified, as highlighted in Figure 3.18a. The percolated clusters are then converted to equivalent node resistor networks. This can be seen in Figure 3.18b, where the particle centers are the nodes and the contact pairs are represented by edges. Potentials ϕ *I* and ϕ *J* are assigned to the nodes and the resistances *R I*,*J* solid,mix to the edges. Finally, the boundary conditions are applied to the boundary nodes and the resulting resistor network can be solved using the scheme according to Section 2.3.

Figure 3.18: Equivalent resistor network to calculate the effective conductivity via the solid volume and surface. a) Percolating, i.e. conducting, pathways. b) Nodes, i.e. potentials, at the centers of the core-shells and edges, i.e. resistances, at the contacts, i.e. overlapping core-shells. c) Connection to the current collectors and boundary nodes.

Note that the transport angles in Equation (3.18) have to be determined in the same way as in the shell model from Section 3.1.4.

### **Model description**

Again, similar to Section 3.1.2, randomly distributed, densely packed and overlap-free assemblies of spheres with a polydisperse particle-size distribution were generated using the RCP. Afterwards, the initial structures were further densified using the numerical sintering algorithm until a certain threshold was reached, i.e. the mean contact angle.

As the FEM model is concerned, the generated sphere structures were first imported into the simulation domain. In order to convert the spheres into coreshell particles, a shell thickness was provided. As a result, every sphere was divided into the shell and the core part. In Figure 3.20b, such a core-shell particle assembly can be observed. The parts of the particles outside of the borders in z-direction of the simulation domain were cut off. For the finite element calculation, a temperature gradient of ∆*T* = 1 was applied to the opposing surfaces of the domain in z-direction on the cutting planes.

The resulting heat flux φ FEM <sup>q</sup> was evaluated at the top or bottom surface in zdirection and the effective thermal conductivity was calculated as

$$
\lambda\_{\text{core}} = \frac{\Phi\_{\text{q}}^{\text{FEM}}}{\Delta T} \frac{L\_{\text{domain}}}{A\_{\text{domain}}},\tag{3.20}
$$

where the domain length and domain cross-section area is *L*domain and *A*domain.

The generated sphere assemblies from above were also imported to the RN framework. For the individual core-shell model resistance, as in Equation (3.18), the shell thickness has to be provided. Similar to the FEM model, the potential drop ∆ϕ = 1 was applied to the boundary nodes of the opposing surfaces in zdirection. The boundary nodes were generated as the contact points of the coreshells particles and the boundary plane in z-direction und, thus, lying exactly on that plane. The resistance between a boundary node and a node within the simulation domain was calculated as a parallel connection, again, according to

$$\mathcal{R}\_{\text{solid,mix}}^{I,0} = \left(\frac{1}{\frac{1/l\_{\text{solid}}^I s^I}{2\pi} \ln\left(\frac{\tan(\theta\_l^{I,0}/2)}{\tan(\theta\_c^{I}/2)}\right) + \frac{1/k\_{\text{stone}}^I}{4r\_c^{I,0}}}\right)^{-1} \,. \tag{3.21}$$

where θ *I*,0 *<sup>c</sup>* and θ *I*,0 *t* is the contact and transport angle between boundary plane and shell, respectively, and *r I*,0 *<sup>c</sup>* is the contact radius, see Figure 3.19.

Figure 3.19: Construction of boundary node and boundary resistance.

Finally, the RN provides the effective electric conductivity using

$$\kappa\_{\rm core} = \frac{I\_{\rm eff}}{\Delta \varphi} \frac{L\_{\rm domain}}{A\_{\rm domain}},\tag{3.22}$$

where *I* eff was calculated using the resistor network solving scheme from Section 2.3.

#### **Test cases**

Up to 15 test cases were created and, based on the test cases, the RN results were compared to the ones from FEM simulations. Three types of assemblies were generated. As can be seen in Table 3.5, the assemblies differ in the standard deviation of the radii ranging from 0.0, i.e. mono-sized particles, to 0.20, i.e. polydisperse particle sizes. In every one of the cases, the shell thickness ratio was taken as *c* = *s*/*r*−*s* = 0.1 representing a moderate shell thickness from the above investigation. In Figure 3.20, the RN model is compared to the same by the FEM model. Here, the assembly is of the type 3 with the largest standard deviation *r*<sup>σ</sup> .


Table 3.5: Structural parameters of the spherical packings for the verification of the RN transport through the solid core-shells.

Figure 3.20: Verification models of the solid-mix resistor network method for an assembly of coreshell particles. a) Potential distribution through the network of an assembly of core-shell particles. b) Temperature distribution of the spatially resolved finite element mesh.

#### **Evaluation and discussion**

Before comparing the effective transport properties, they are normalized by th bulk conductivities first. In this case, the bulk conductivities were taken from the core parts. This yields the dimensionless effective transport parameters as

$$
\hat{k}\_{\rm eff,mix}^{\rm FEM} = \frac{\lambda\_{\rm eff}^{\rm FEM}}{\lambda\_{\rm core}} \text{and} \quad \hat{k}\_{\rm eff,mix}^{\rm RN} = \frac{\kappa\_{\rm eff}^{\rm RN}}{\kappa\_{\rm core}}, \tag{3.23}
$$

respectively, where λcore is the thermal and κcore is electric conductivity of the cores. The results of the above investigation can be seen in Figure 3.21. The effective conductivity is overestimated by the RN. However, the mean error is around ¯*e* = 7% which is acceptable given the fact that the core-shell model is a rough simplification. In this sence, the core-shell model can be regarded as verified to a degree sufficient of this work.

Figure 3.21: Verification of the solid-mix resistor network method using finite element simulations.

# **3.2 Effective transport properties of the pore phase**

In this section, the effective transport properties of the pore phase are investigated. The pore phase is considered to be the space between the solid phase, e.g. active material particles in lithium-ion battery cathodes. Here, this space is filled with liquid electrolyte which is needed to provide sufficient lithium-ions and ionic charge for the electrochemical reactions on the active material particle surfaces. In other words, the electrolyte-filled pore phase has to provide a sufficiently large effective diffusivity and conductivity for a lithiumion battery to operate properly.

In addition to the electrolyte, inside the pore phase an electronically conductive carbon-black-binder phase, i.e. filler phase, is distributed. This is in particular necessary since the active material particles are commonly known for their poor electronic conductivity. This, on the other hand, is needed to provide enough electrons for the electrochemical reactions. The pore phase is thus shared by both the liquid electrolyte and the filler phase.

In the following, a resistor network method is developed for the calculation of the effective conductivity of the pore phase. Additionally, an extension is proposed to account for those cases where the pore phase is shared by multiple conducting phases.

# **3.2.1 Volume resistance of pore throats**

The crucial part of the resistor network approach is the identification and computation of individual resistances. In the case of the pore phase, the individual resistances are the pore throats. Therefore, a model is proposed to estimate the resistance of such pore throats based on the geometrical bottleneck effect [91, 92]. For simplicity, the model is presented for the 2D case. The steps necessary for a 3D implementation are mentioned whenever needed.

### **Model**

To start from, inside the pore phase of a granular assembly, it is assumed that transport is hindered by the bottlenecks created by the solid particles. To model such bottlenecks, i.e. pore throats, it is necessary to identify the pores and pore centers first.

Figure 3.22: Pore-volume resistance. a) Flux through a bottleneck between two particles. b) Resulting flux between nodes *N I* and *N J* inside the spatially discretized bottleneck, i.e. pore throat. c) Discretization of throat into sub-throats. Divide sub-throats into small increments in order to calculate wire resistances.

Figure 3.22a shows an example of such a bottleneck where two particles form a narrow region, i.e. a pore throat, between them. Before and after the pore throat are the pore centers which are connected by the pore throat. Between the pore centers, some flux ~*F* is considered to flow through. The exact position of the pore centers is achieved by a spatial decomposition technique, which will be subject to the next part. Next, in Figure 3.22b, the nodes *N I* and *N J* indicate the start and end point of the flux. Those nodes are taken as the pore centers and the gray-shaded region is the spatial decomposition of the pore throat.

A pore throat resistance *R I*,*J* pore,vol is attributed to the pore throat which is expressed as a parallel connection of resistances representing the sides of the node connecting edge. In the 2D example, see Figure 3.22b, these sides are the sub-regions created by triangles of the edge defining nodes and the respective centers of the surrounding particles. Additionally, the overlapping areas of the particles are subtracted from the triangles. Because the sub-regions describe parts of the pore throats, they are called sub-throats.

A sub-throat is sketched in Figure 3.22c. The resistance of such a sub-throat is calculated as a series connection of wire resistances. To this end, each edge is divided into sufficiently small increments ∆*L n* and the resistance ∆*R n <sup>m</sup>* of each increment is calculated by

$$
\Delta R\_m^n = \rho\_{\text{pore}} \frac{\Delta L^n}{A\_{\text{mean}}^n} \,, \tag{3.24}
$$

where ρpore is the resistivity of the pore material. A mean line length *A n* mean of the *n*'th increment has to be passed, in case of 2D. In case of 3D, its a mean cross section area.

The resulting resistance *R I*,*J <sup>m</sup>* of one sub-throat *m* is then calculated as a series connection of the incremental resistances as

$$R\_m^{I,J} = \sum\_{n=1}^{n\_{\text{incr}}} \Delta R\_m^n = \rho\_{\text{pore}} \sum\_{n=1}^{n\_{\text{incr}}} \frac{\Delta L^n}{A\_{\text{mean}}^n} \,, \tag{3.25}$$

where *n*incr is the number of increments between nodes *N I* and *N J* . Here, ∆*L* is estimated individually for each sub-throat. The size of ∆*L* is successively reduced until the resulting resistance calculated for the next increment doesn't differ from the previous one by a certain small threshold, here within 5%.

Finally, as mentioned above, the resulting resistance

$$R\_{\text{pore,vol}}^{I,J} = \left(\sum\_{m=1}^{m\_{\text{stbr}}} \frac{1}{R\_m^{I,J}}\right)^{-1} \tag{3.26}$$

of the edge, i.e. pore throat, which connects the nodes *N I* , *N J* , i.e. the pore centers, is a parallel connection of the resistances *R I*,*J <sup>m</sup>* of all *m*sthr sub-throats belonging to this edge. Thus, a model is derived for pore throat resistances which can be used along the resistor network method when modeling the pore phase effective transport properties. The index "pore,vol" refers to transport through the volume of the pore phase which is why this resistance is called pore-volume resistance.

# **3.2.2 Effective transport properties of pore throat assemblies**

In the above section, the computation of an individual throat resistance was presented, i.e. the pore-volume resistance from Equation (3.26). Next, this kind of resistance is used in the framework of the resistor network method to calculate effective transport properties of the pore phase.

### **Resistor network method for the pore phase**

While it is quite straight forward to discretize the solid phase of a granular system using spheres, as demonstrated in Section 3.1, it is not as obvious to model the pore phase in a similar way, because of the relatively complex shapes of the pores and throats. In Figure 3.23, exemplarily, an equivalent resistor network for the pore phase is sketched.

First, pores and pore connecting throats have to be identified. To this end, the pore phase is discretized using the so-called Laguerre or generalized Voronoi tessellation method [122, 123], where the generators are given by the midpoints of the particles and their weights are given by the corresponding radii. For detailed information regarding this tessellation technique, refer to [124]. The computation is done by employing the software library Voro++ [125].

Figure 3.23: Equivalent resistor network to calculate the effective conductivity of the pore volume phase. a) Discretization of the pore phase using Laguerre cells. b) Potentials at cell vertices and resistances at the edges. c) Connection to the current collectors and boundary nodes.

In general, for isolated points in a domain, which happen in this case to be the sphere centers, the tessellation method assigns a spatial cell to each such point. Each cell contains all points whose distance to the associated cell center is less than or equal to any other sphere or cell center. In case of packings with polydisperse particle sizes, additionally, the distances are adjusted with respect to the particle sizes [122]. In 3D space, the cells take the form of convex polyhedra. Using this method, the whole domain volume can be fully decomposed into cells. Finally, the pore phase is defined as the volume of the cells where the volume of the particles, i.e. solid phase, is subtracted from.

The result of the tessellation is that all particles are wrapped in cells while their borders, i.e. edges and faces, lie in an optimum distance between the cell centers, i.e. particle centers. An exemplary cell discretization can be seen in Figure 3.23a. The spherical particles are wrapped by cells, edges, vertices and, in the 3D-case, faces. Each vertex is chosen as the center of a pore and the corresponding edges as the pore throats. By extension, the network nodes *N I* and *N J* correspond to the vertices where the potentials ϕ *I* , ϕ *J* are associated with and the edges are related to the resistances *R I*,*J* pore,vol, see Figure 3.23b.

In order to compute the resistances using the approach described above, the regions around the individual edges have to be decomposed into throats and subthroats. In case of 2D, the regions are constructed as the areas which are given by triangles defined by the nodes of the considered edge and the associated sphere centers, where the intersection area with the associated spheres is subtracted from. In the 3D case, the surrounding region is the volume given by tetrahedrons, which, as an extension of the 2D case, are additionally defined by the centers of the faces of the cells which meet at the considered edge. The intersecting volumes of the associated spheres are then subtracted from the tetrahedrons.

The boundary nodes have to be identified where the boundary conditions can be applied to, see Figure 3.23c. Finally, the resulting node resistor network is solved by the scheme described in Section 2.3.

### **Verification**

As a next step, the verification of the pore-volume model in the framework of the resistor network method is checked. To this end, similar to Section 3.1.2, different assemblies of varying degree of particle polydispersity were generated using the RCP and densified using the numerical sintering algorithm. The resulting collections of spheres were then imported to both the RN and the FEM framework. Finally, the resulting effective transport properties were computed and compared to each other.

As far as the finite element model is concerned, first, the generated sphere assemblies were imported into a box-shaped domain. Secondly, the spheres were cut out of the domain such that the complement volume remained. The latter of which was identified as the pore phase and can be seen in Figure 3.24b. The color gradient refers to the temperature distribution emerging from the temperature drop of ∆*T* = 1 on the surface nodes of the simulation domain in the z-direction while the other surfaces in x- and y-direction remain insulated. The resulting heat flux φ FEM <sup>q</sup> was evaluated at the top or bottom surface in z-direction and the effective thermal conductivity was calculated as

$$
\lambda\_{\rm eff}^{\rm FEM} = \frac{\Phi\_{\rm q}^{\rm FEM}}{\Delta T} \frac{L\_{\rm domain}}{A\_{\rm domain}},\tag{3.27}
$$

where the domain length and domain cross section area is *L*domain and *A*domain, respectively.

Regarding the resistor network model, the spheres were imported into the boxshaped simulation domain, as well. As a next step, the Laguerre tessellation was used, as described above, to fully discretize the pore phase and which was the basis for the resistor network method. In Figure 3.24a, the potential distribution of such a resistor network can be observed. The vertices on the boundary planes in z-direction were directly used as boundary nodes and an electric potential drop of ∆ϕ = 1 was applied to these boundary nodes.

The effective electric conductivity was calculated by

$$\kappa\_{\rm eff}^{\rm FEM} = \frac{I\_{\rm eff}}{\Delta \varphi} \frac{L\_{\rm domain}}{A\_{\rm domain}} \,, \tag{3.28}$$

where the effective current *I* eff was calculated using the resistor network solving scheme according to Section 2.3.

#### **Test cases**

Next, 15 test cases were created, where the structural parameters upon which the assemblies were generated are presented in Table 3.7. In Figure 3.24, as an example, one assembly of the type 3 with the largest standard deviation is shown.


Table 3.7: Structural parameters of the spherical packings for the verification of the transport through the volume of the pore phase.

Figure 3.24: Verification models of the pore-volume resistor network method for the pore phase of an assembly of spheres. a) Potential distribution through the network of pores and throats. b) Temperature distribution of the spatially resolved finite element mesh.

### **Evaluation and discussion**

In order to compare the RN and FEM results, the effective transport properties have to be converted to effective transport parameters. Therefore, the effective thermal and electric conductivity were normalized by the bulk thermal λbulk and electric conductivity κbulk, respectively, which yields

$$
\hat{k}\_{\rm eff}^{\rm FEM} = \frac{\lambda\_{\rm eff}^{\rm FEM}}{\lambda\_{\rm bulk}} \quad \text{and} \quad \hat{k}\_{\rm eff}^{\rm RN} = \frac{\kappa\_{\rm eff}^{\rm FEM}}{\kappa\_{\rm bulk}}.\tag{3.29}
$$

In Figure 3.25, the comparison of the investigation can be seen. As usual, the finite element results are on the horizontal whereas the resistor network values are on the vertical axis. The solid black line indicates a perfect match of the two methods. Here, the results obtained by the RN are overestimating the reference FEM results by a little bit as the mean error is around *e* = 1%. This, on the other hand, is a low error considering the approximation of the arguably complex pore phase structure by a combination of series and parallel connections of wire resistances. Moreover, concerning resource and time cost, the RN is superior to the FEM calculation [91].

Figure 3.25: Verification of the pore-volume resistor network method using finite element simulations.

Given the rather low error, even for the pore phase of polydispersely sized sphere packings, the resistor network method to calculate the effective conductivity of the pore phase can be considered as verified.

## **3.2.3 Volume resistance of mixed pore throat**

Different other material phases can distributed inside the pore phase. For example, to enhance the electronic conductivity and to fix the active materials of the cathodes of lithium-ion batteries, solid carbon-black (CB) and polymeric binder is added to the pore phase. Thus, the pore phase is shared by liquid electrolyte, on the one hand, and by a conductive filler phase, i.e. carbon-blackbinder phase, on the other hand. In the following, the model of the pore-volume resistance is extended in such a way that it is able to account for such cases. It is assumed that the different phases in the pore phase don't interact with each other.

### **Model**

In Figure 3.26, the bottlenecks, similar to Section 3.2.1, are sketched. The only difference is that another phase is added to the bottleneck space. The added phase is indicated by the black dots in Figure 3.26a. Next, in Figure 3.26b, the pore throat is identified as the gray-shaded area.

Figure 3.26: Pore-mix resistance. a) Flux through a bottleneck between two particles. Additional phase depicted as black dots. b) Resulting flux between nodes *N I* and *N J* inside the spatially discretized bottleneck, i.e. pore throat. c) Discretization of throat into subthroats. Divide sub-throats into small increments in order to calculate wire resistances. Account for volume fraction of the additional phase.

The resistance of such a pore throat with added phases is *R I*,*J* pore,mix,*c* , where *c* refers to the transport through the particular phase *c* of a mixed pore phase. In order to calculate the resistance, again, each edge is divided into sufficiently small increments ∆*L <sup>n</sup>* where an incremental resistance ∆*R n*,*c <sup>m</sup>* of each segment is calculated by

$$
\Delta R\_m^{n,c} = \rho\_{\text{pore}}^c \phi^c \frac{\Delta L^n}{A\_{\text{mean}}^n} \,. \tag{3.30}
$$

Here, ρ *c* pore is the resistivity of the phase *c*, φ *c <sup>m</sup>* is the volume fraction of the phase inside the sub-throat *m* and *A n* mean is the mean cross section area for the *n*'th increment.

The resulting resistance *R I*,*J*,*c <sup>m</sup>* of the added phase *c* inside the sub-throat *m* is then calculated as a series connection of the incremental resistances as

$$R\_{m}^{I,J,c} = \sum\_{n=1}^{n\_{\rm incr}} \Delta R\_{m}^{n,c} = \rho\_{\rm pore}^{c} \phi^{c} \sum\_{n=1}^{n\_{\rm incr}} \frac{\Delta L^{n}}{A\_{\rm mean}^{n}},\tag{3.31}$$

where the number of increments used between Nodes *N I* and *N J* is *n*incr and ∆*L* is chosen individually for each sub-throat. For this, ∆*L* is successively reduced until the resulting resistance of the sub-throat doesn't differ from the previous one by a certain small threshold. In this work, the threshold is set to 5%.

The resulting resistance *R I*,*J* pore,mix,c of the added phase *c* inside the pore throat is a parallel connection of the resistances *R I*,*J*,*c <sup>m</sup>* of all *m*sthr sub-throats.

$$\mathcal{R}^{I,J}\_{\text{pore,mix,c}} = \left(\sum\_{m=1}^{m\_{\text{str}}} \frac{1}{R\_m^{I,J,c}}\right)^{-1} \,. \tag{3.32}$$

Summarizing, for each phase of the pore phase separately, a pore throat resistance was derived where the transport path is shared by more than one phases. The index "pore,mix" refers to transport through the mixed volume of the pore phase which is why this resistance is called pore-mix resistance.

# **3.2.4 Effective transport properties of mixed pore throat assemblies**

In the following, the pore-mix resistance is used inside the framework of the RN. To this end, a pore-mixed resistor network is presented and validated using experimental results from the literature.

### **Resistor network for the pore mix resistance**

The build up of the resistor network for the mixed pore structure is similar to Section 3.2.2. As can be seen in Figure 3.27a, the Laguerre tessellation is used to discretize the pore phase. Recall that the resulting cells are formed by the vertices, edges and faces. The black dots represent an additional randomly distributed phase inside the pore phase.

Figure 3.27: Equivalent resistor network to calculate the effective conductivity of the pore mix phase. a) Discretization of the pore phase using Laguerre cells. b) Potentials at the Laguerre vertices and representative resistances at the edges. c) Account for volume fraction of added phase in each sub-throats. Connection to the current collectors and boundary nodes.

In Figure 3.27b, it can be seen that the vertices are taken as the nodes where the potentials ϕ *I* , ϕ *J* are associated with. The edges stand for the representative poremix resistance *R I*,*J* pore,mix. Its computation is performed using the methodology described in Section 3.2.3. To this end, the distributed additional phase inside the pores is regarded as smeared across the sub-throats. This is depicted by the differently shaded areas, i.e. sub-throats, in Figure 3.27c. Note that the differently shaded areas imply that, if known, each sub-throat can have individual volume fractions of the considered phases. Finally, the boundary conditions can be applied to the boundary nodes and the effective conductivity can be calculated using the scheme described in Section 2.3.

### **Validation**

In the above sections, the FEM was rigorously used to verify the individual resistor network approach. However, in case of lithium-ion batteries, this is not feasible due to the very complex structure of the mixed phases in the pore phase. Alternatively, the resistor network approach for the mixed pore resistance is validated using experiments from literature. In [126], the effective electric conductivity of carbon-black-binder mixtures with and without the presence of active material particles was measured. Several binder compositions were presented. In this case, the binder consisted of polyethylene oxide (PEO), polyvinylidene difluoride-co-hexafluoroproylene (PVdF-HFP) and ethylene carbonate-propylene carbonate (EC-PC). The carbon-black content of the mixtures was varied from 0 to 18.8%.

It can be seen that for a carbon-black volume fraction of 0%, the effective electric conductivity does not vanish. That is because the ionic conductivity of the binder contributes to the electric conductivity, which is around 1 · 10−<sup>7</sup> Sm−<sup>1</sup> . Since this compositions have no relevance in applications, in this investigation, carbonblack-binder mixtures and electrode compositions without carbon-black content are neglected. Therefore, if the carbon-black volume fraction reaches zero, the electric conductivity is also set as zero. This way, the electronic conductivity is extracted from the experimental data. The effective electric conductivity rises orders of magnitudes for carbon-black volume fractions from 1.5 to 6%. This range was identified as the so-called percolation threshold, which resembles the critical volume fraction of conducting species where the likelihood to form conducting pathways increases drastically. For contents larger than that, the effective electric conductivity increases asymptotically until, at the theoretical limit of carbon-black volume fraction of 100%, the bulk electronic conductivity of carbon-black is reached, which is assumed as 1000Sm−<sup>1</sup> [126].

In order to account for the above described behavior, the effective electronic conductivity of the carbon-black-binder mixture provided by the experimental data was fitted by the s-shaped function

$$\kappa\_{\rm CB-B,eff}^{\rm fit} = a \cdot \left(1 - \exp\left(-\left(\frac{\phi\_{\rm CB}}{b}\right)^c\right)\right),\tag{3.33}$$

where the parameters are *a* = 1000, *b* = 0.167 and *c* = 2.358. Here, *a* was chosen as the bulk electronic conductivity of carbon-black. Thus, for large values of the carbon-black volume fraction φCB related to the carbon-black-binder mixture the bulk electronic conductivity would be reached. The fit-function fulfills the criterion that the effective electric conductivity is 0Sm−<sup>1</sup> if the carbon-black content reaches zero, see Figure 3.28b.

As the next step, active material particles were added to the carbon-black-binder mixtures and the effective electric conductivity was measured. In the following, the electrode compositions are divided into two phases, which are the pore and the solid phase. The former of which comprises of the electrolyte and the carbonblack-binder mixture phase and the latter of which is the active material phase. It was reported that the electrode samples were densified such that porosities of 30% were achieved. Since the pores are filled with electrolyte, the resulting phase is identified as the electrolyte phase. The remaining 70% of volume fraction were shared by the carbon-black-binder mixture and the active material phase. This phases contained a volume fraction of 19 − 22% active material and 81 − 78% of carbon-black-binder mixture. For the following investigation, the active material and carbon-black-binder volume fraction with respect to the whole volume was set as the average values, which were around 15% and 55%, respectively.

Experimentally achieved effective electric conductivity values of the above described electrode structures were compared to numerical results. To this end, comparable virtual compositions were created. Similar to Section 3.1.2, the initial sphere packings were generated using the RCP. The structures were further densified by means of the numerical sintering algorithm. Additionally, the algorithm was slightly modified. In this case, the box-shaped simulation domain was isotropically and successively reduced or expanded such that the spheres inside the domain were moved upon or apart from each other. This was done until the targeted porosity was achieved. This way, the particle size was preserved.

The used active material was Li1,2V3O<sup>8</sup> and the size was set to a particle radius of 1 µm, which was deduced from SEM images in [126]. Following the description of the electrode composition from above, the targeted volume fraction of the active material with respect to the whole volume was set to 15%. The remaining pore phase of 85% was shared by 65% and 35% carbon-black-binder mixture and electrolyte, respectively.

Based on this composition, scenarios of sphere assemblies were created which were randomly distributed but structure-wise identical to the electrode in the experiments. Using the pore-mix resistor network method, the presence of added phases, i.e. electrolyte and carbon-black-mixture phase, was accounted for by appropriately setting the volume fraction *c* in the pore throat resistance computation from Equation (3.26). In order to calculate the (effective) bulk conductivity of the carbon-black-binder phase, the fit-function from Equation (3.33) was used.

Figure 3.28: Validation of the pore-mix resistor network method by using experiments from [126]. a) Effective electric conductivity of a cathode structure with active material and carbonblack-binder phase. b) Effective electric conductivity of the carbon-black-binder phase as a function of the carbon-black content, see Equation (3.33).

In Figure 3.28a, the results of the resistor network calculations, i.e. κ RN eff,CB-B, are compared to the measurements from [126]. It can be seen that the values provided by RN match the experiments very good for carbon-black volume fractions above 4%, which is within the observed percolation threshold [126]. However, the numerical results deviate a little bit for carbon-black contents below 4%. This originates from the fact that the fit-function used for the conductivity of the carbon-black-binder phase deviates for lower carbon-black volume fractions. On the one hand, the fit-function was chosen such that it represents the s-shaped percolation phenomenon, where it is difficult to properly describe the behavior at the vicinity of the percolation threshold. On the other hand, usually, carbon-black volume fraction are above 4%, see for instance [127]. In general, the overall quality of the match between the numerical results and the experiments leads to the conclusion that the resistor network approach described above is able to calculate the effective conductivity of added conductive phases in the pore phase.

# **3.3 Summary**

In Sections 3.1 and 3.2, it was shown how the resistor network method can be used to model effective transport properties of both the solid and pore phase in an assembly of spheres with polydisperse size-distribution.

Concerning the solid phase, equivalent resistor networks were created based on boundary connecting clusters of spheres. The transport through this clusters was distinguished between volume and surface type. In order to be used in the framework of the RN, the crucial part is the computation of resistances of individual contact pairs. To this end, regarding the volume transport, the solidvolume resistance was constructed in Equation (3.3) based on the sphere model. The surface transport was computed using the solid-surface resistance from Equation (3.11) where the underlying model was the shell model. Moreover, assuming the transport is partially through the volume and the surface, a solidmix resistance was developed in Equation (3.18). The corresponding model was the core-shell model.

As for the pore phase, equivalent resistor networks were created based on spatial decomposition techniques. Here, the Laguerre tessellation was used where every particle inside an assembly is placed into convex polyhedral cells. This way, the pore phase was geometrically defined as the complement volume of cells and the contained spheres. Additionally, the cell vertices and edges were the basis upon which resistor networks were created. While the vertices resembled the pore centers, the edges were identified as the pore throats. Based on the pore throat geometry, a resistance was calculated, accordingly. Since the transport was through the volume of the pore, the corresponding resistance in Equation (3.26) was named pore-volume resistance. Additionally, since the pore phase can be shared by multiple conducting phases, another model was proposed. The poremix resistance in Equation (3.32) considers the resistance of additional phases using their volume fractions and resistivity.

# **4 Cell modeling**

The working principle of lithium-ion battery cells was described in Section 1.1. In general, electrochemical reactions take place on the surfaces of the active material of the positive and negative electrode. Those electrodes usually are porous structures and comprise of different components. In order to model battery cells, several physics-based approaches on different length-scales can be applied [10]. In the following, a continuum approach is chosen.

A prominent candidate of this modeling approach is the Newman type cell model [41]. While the development of novel cathode structures proceeds, i.e. hierarchically structured electrodes [26–29], the need for novel cell modeling approaches arises, as well. Therefore, in the first section of the following part, the classical cell model is presented first. Next, In the second section, the cell model is consistently extended to account for hierarchically structured electrodes.

## **4.1 Classical half-cell**

In Figure 4.1a, the sketch of the classical half-cell setup is shown. The anode on the left-hand side is a lithium foil and the positive electrode, i.e. cathode during the discharging process, on the right-hand side is a porous composite structure, comprising the active material and carbon-black-binder mixture, which is the conductive filler material. The porous separator is sandwiched between the lithium foil and the positive electrode. In addition, both the porous cathode and separator are completely soaked with liquid electrolyte. The setup is chosen such that one is able to study solely the influence of the porous cathode on the performance of the cell. In contrast, the full-cell setup, with an additional porous anode, both electrodes would influence the cell performance.

Figure 4.1: Half-cell setup of lithium-ion batteries. a) Sandwich structure of a lithium-ion battery cell. The anode is a solid lithium metal. b) The porous structure of the cathode is homogenized.

A suitable model to incorporate battery charging and discharging processes originated from [36, 40–42]. The model basis is the so-called porous electrode theory [47] and [35]. In this theory, the structural details of the actual geometry is smeared out across the model to achieve effective properties. This way, two scales are introduced, where transport defined on the full three dimensional structural details of the electrode on the smaller scale is called microscopic while the corresponding effective transport on the larger smeared scale is referred to as macroscopic. Therefore, transport equations defined on the microscale are referred to as microscale equations whereas transport equations defined on the macroscale are referred to as macroscale equations. Every point of the model represents a superposition of two phases. While the electrolyte represents one phase by itself, typically, the active material and filler material is subsumed as the other phase, i.e. the solid phase. The structural and material properties of each phase are characterized by the volume fraction, specific surface area and effective transport properties. Additionally, the two phases interact with each other via electrochemical reactions on the active material surfaces.

Thus, the particle and pore geometry of the system is not spatially resolved but, rather, can be viewed as smeared across the model. In Figure 4.2b, the smeared region is indicated by the shaded area where, additionally, at every point in the cathode region the active material particles are modeled as spheres. Therefore, the cell model, as presented here, is divided into two levels, i.e. the cell and the particle level. The former of which represents the macroscopic scale, i.e. transport in the porous electrode, whereas the latter describes the transport processes inside the active material.

## **4.1.1 Macroscale equations**

In the framework of the porous electrode theory, at every spatial point of the macroscale model, the solid and the electrolyte phase are superimposed. By extension, the currents and fluxes in the phases are superimposed, as well. It is assumed that, on cell level, the electronic current is carried by the solid phase, whereas the ionic current is carried by the electrolyte phase. Furthermore, while the transport of cations is supposedly via the electrolyte phase on cell level, the insertion process and lithium transport process inside the solid phase is modeled on particle level.

The above transport phenomena are described by four continuity equations representing the conservation of electronic and ionic charge as well as the conservation of mass in the solid and electrolyte phase. In the following, the micro- and macroscale forms of the continuity equations are recalled as originally presented in [35, 41, 47]. As it was shown in [128, 129], the macroscopic forms used in the cell models can be derived by applying volume averaging methods, which were presented in Section 2.2. The same methods are also chosen in the following.

### **Conservation of mass in the solid phase**

First, the transport in the active material is considered. At any point of the electrode, lithium is intercalated into the active material. The inserted lithium is assumed to obey a Fickian type diffusion transport process which can be written as

$$\frac{\partial c\_s}{\partial t} = \nabla \cdot \left( D\_s \nabla c\_s \right) \,, \tag{4.1}$$

where *c<sup>s</sup>* is the lithium concentration and *D<sup>s</sup>* is the corresponding diffusion coefficient of active material. In addition, the active material secondary particles are assumed to be spheres such that spherical coordinates can be employed. Moreover, the transport is modeled as spherically symmetric, such that, using the transformation operation from Appendix A.1, Equation (4.1) can be rewritten in spherically symmetric form yielding

$$\frac{\partial c\_s}{\partial t} = \frac{1}{\mathbf{y}^2} \frac{\partial}{\partial \mathbf{y}} \left( \mathbf{y}^2 D\_s \frac{\partial c\_s}{\partial \mathbf{y}} \right). \tag{4.2}$$

Note that the dimension *y* is the radial coordinate inside the spherical particles.

#### **Conservation of electronic charge in the solid phase**

In the porous electrode, the electronic charge is carried by the solid phase. The microscale electronic current density is described using Ohm's law as

$$
\vec{\mathbf{i}}\_s = -\mathbf{k}\_s^{\text{eon}} \nabla\_\xi \boldsymbol{\varrho}\_s,\tag{4.3}
$$

where κ eon *s* is the electronic conductivity and ϕ*<sup>s</sup>* is the electronic potential of the solid phase. Note that, analogous to Section 2.2, the coordinate ξ refers to the microscale. Due to conservation of charge, the continuity equation reads

$$\nabla\_{\xi} \cdot \vec{\mathbf{i}}\_{s} = \nabla\_{\xi} \cdot \left( -\mathsf{k}\_{s}^{\mathrm{eon}} \nabla\_{\xi} \varphi\_{s} \right) = \mathbf{0}. \tag{4.4}$$

From Section 2.2, using Equation (2.24), the macroscopic form of Equation (4.4) is

$$\nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{eon}} \nabla \overline{\!\! \! \! \! /} \overline{\!\! \! \! /} \right) - a\_{\text{se}} \, \overline{j}\_{\text{se}} \, \mathcal{F} = 0 \,,\tag{4.5}$$

where, comparing to Equation (2.24), *k*α,eff ≡ κ eon *<sup>s</sup>*,eff is the effective electronic conductivity, *p*<sup>α</sup> ≡ ϕ*<sup>s</sup>* is the macroscopic electronic potential and *a*αβ ≡ *ase* is the specific surface area of the solid phase. Electronic charge is produced due to electrochemical reactions at the interface between active material and electrolyte. This is accounted for by the reaction term *b*αβ ≡ *j se*F which represents the exchange current density induced by lithium flux *j se* at the interface between solid and electrolyte phase, where F is the faraday constant.

Rearranging Equation (4.5) leads to

$$\nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{en}} \nabla \overline{\varphi}\_s \right) = a\_{se} \overline{j}\_{se} \mathcal{F}, \tag{4.6}$$

which is the macroscopic continuity equation representing the electronic charge conservation in the solid phase.

### **Conservation of ionic charge in the electrolyte phase**

The ionic charge is carried by the electrolyte phase. Using the concentrated solution theory on binary electrolytes [47], the ionic current density can be written as

$$\overrightarrow{\mathbf{i}}\_{\varepsilon} = -\kappa\_{\varepsilon}^{\rm ion} \nabla\_{\xi} \boldsymbol{\upvarphi}\_{\varepsilon} - \kappa\_{D}^{\rm ion} \nabla\_{\xi} \ln c\_{\varepsilon},\tag{4.7}$$

where κ ion *e* is the ionic conductivity, ϕ*<sup>e</sup>* is the ionic potential and *c<sup>e</sup>* is the lithium concentration of the electrolyte phase. Moreover, the diffusional conductivity is

$$\kappa\_D^{\rm ion} = \frac{\nu \kappa\_e^{\rm ion} \mathcal{R} T}{\mathcal{F}} \left( \frac{s\_+}{n\nu\_+} + \frac{t\_+^0}{z\_+ \nu\_+} - \frac{s\_0 c\_e}{n c\_0} \right) \left( 1 + \frac{\partial \ln f\_\pm}{\partial \ln c\_e} \right), \tag{4.8}$$

where *s*+, *s*<sup>0</sup> and ν<sup>+</sup> are stoichiometric coefficients of the cations in the solute, *z*<sup>+</sup> is the charge number of the cations, *t* 0 <sup>+</sup> is the transference number of the cations with respect to the solvent, *n* is the number of moles of a species and *c*<sup>0</sup> is the concentration of the solvent. In case of LiPF6, which is the electrolyte salt used in this work, the reaction equation is LiPF<sup>6</sup> −−)−−\* Li<sup>+</sup> +PF<sup>6</sup> <sup>−</sup> [10] and, therefore, the parameters are *s*<sup>+</sup> = −1, *s*<sup>0</sup> = 0, *n* = 1, ν<sup>+</sup> = 1, ν = 2 and *z*<sup>+</sup> = 1. Note that ν = ν<sup>+</sup> +ν−, where ν<sup>−</sup> = 1 is the a stoichiometric coefficient of the anions. Equation (4.8) can therefore be simplified as

$$\kappa\_D^{\rm ion} = \frac{2\kappa\_e^{\rm ion}\mathcal{R}\mathcal{T}}{\mathcal{F}} \left(t\_+^0 - 1\right) \left(1 + \frac{\partial \ln f\_\pm}{\partial \ln c\_e}\right). \tag{4.9}$$

Due to conservation of charge, the continuity equation reads

$$
\nabla\_{\xi} \cdot \vec{\mathbf{i}}\_{\varepsilon} = \nabla\_{\xi} \cdot \left( -\kappa\_{\varepsilon}^{\rm ion} \nabla\_{\xi} \varphi\_{\varepsilon} \right) + \nabla\_{\xi} \cdot \left( -\kappa\_{D}^{\rm ion} \nabla\_{\xi} \ln c\_{\varepsilon} \right) = 0 \,, \tag{4.10}
$$

From Section 2.2, using Equation (2.25) and recognizing that the electrolyte represents the β-phase, the macroscopic form of Equation (4.10) is

$$\nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion}} \nabla \overline{\Phi}\_e \right) + \nabla \cdot \left( -\kappa\_{D, \text{eff}}^{\text{ion}} \nabla \ln \overline{c}\_e \right) + a\_{se} \overline{j}\_{se} \mathcal{F} = 0 \,. \tag{4.11}$$

where, comparing to Equation (2.25), the effective conductivity *k* <sup>β</sup>,eff is identified as either the effective ionic conductivity κ ion *<sup>e</sup>*,eff or the effective diffusional conductivity

$$\mathbf{k}\_{D, \text{eff}}^{\text{ion}} = \frac{2\mathbf{k}\_{e, \text{eff}}^{\text{ion}} \mathcal{R} \mathcal{T}}{\mathcal{F}} \left( 1 + \frac{\partial \ln f\_{\pm}}{\partial \ln \mathbb{Z}\_e} \right) \left( t\_+^0 - 1 \right) \,. \tag{4.12}$$

Moreover, the macroscopic potential *p*<sup>β</sup> is either the macroscopic ionic potential ϕ*e* or the macroscopic concentration of lithium *c<sup>e</sup>* of the electrolyte phase. The specific surface area *a*αβ ≡ *ase* is between the solid and the electrolyte phase, where ionic charge is produced due to electrochemical reactions. This is accounted for by the reaction term *j se*F, as well. Rearranging Equation (4.11) yields the macroscopic continuity equation representing conservation of ionic charge in the electrolyte phase as

$$\nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion}} \nabla \overline{\varphi}\_e - \kappa\_{D, \text{eff}}^{\text{ion}} \nabla \ln \overline{c}\_e \right) = -a\_{\text{se}} \overline{j}\_{\text{se}} \mathcal{F} \,. \tag{4.13}$$

#### **Conservation of mass in the electrolyte phase**

Following the above assumed concentrated solution theory of binary electrolytes [47], the flux density of cations is

$$\vec{J}\_{+} = -\mathsf{v}\_{+} D\_{e} \left( 1 - \frac{\mathrm{d} \ln c\_{0}}{\mathrm{d} \ln c\_{e}} \right) \nabla\_{\xi} c\_{e} + \frac{\vec{\mathsf{i}}\_{e} t\_{+}^{0}}{z\_{+} \mathcal{F}} + c\_{+} \vec{\upsilon}\_{0},\tag{4.14}$$

where *c<sup>e</sup>* is the salt concentration of the electrolyte or, in short, electrolyte concentration, *D<sup>e</sup>* is the diffusion coefficient *c*<sup>0</sup> is the concentration of the solvent and ~*v*<sup>0</sup> is the velocity of the solvent. The corresponding continuity equation reads

$$\begin{split} \frac{\partial c\_{+}}{\partial t} &= -\nabla\_{\xi} \cdot \vec{\mathcal{J}}\_{+} \\ &= \nabla\_{\xi} \cdot \left[ \mathbf{v}\_{+} D\_{\epsilon} \left( 1 - \frac{\mathbf{d} \ln c\_{0}}{\mathbf{d} \ln c\_{\epsilon}} \right) \nabla\_{\xi} c\_{\epsilon} - \frac{\vec{\mathcal{I}}\_{\epsilon} t\_{+}^{0}}{z\_{+} \mathcal{F}} - c\_{+} \vec{\mathcal{v}}\_{0} \right] \\ &= \nabla\_{\xi} \cdot \left[ \mathbf{v}\_{+} D\_{\epsilon} \left( 1 - \frac{\mathbf{d} \ln c\_{0}}{\mathbf{d} \ln c\_{\epsilon}} \right) \nabla\_{\xi} c\_{\epsilon} \right] - \frac{\nabla\_{\xi} \cdot \left( \vec{\mathcal{I}}\_{\epsilon} t\_{+}^{0} \right)}{z\_{+} \mathcal{F}} - \nabla\_{\xi} \cdot c\_{+} \vec{\mathcal{v}}\_{0} . \end{split} \tag{4.15}$$

Using electroneutrality conditions, i.e. *c<sup>e</sup>* = *c*+/ν<sup>+</sup> = *c*−/ν<sup>−</sup> or *c*<sup>+</sup> = *ce*ν+, Equation (4.15) can be rearranged as

$$\mathbf{v}\_{+}\frac{\partial c\_{e}}{\partial t} = \mathbf{v}\_{+}\nabla\_{\xi} \cdot \left[ D\_{e} \left( 1 - \frac{\mathbf{d}\ln c\_{0}}{\mathbf{d}\ln c\_{e}} \right) \nabla\_{\xi} c\_{e} \right] - \frac{\nabla\_{\xi} \cdot \left( \overleftarrow{\mathbf{i}}\_{e} t\_{+}^{0} \right)}{z\_{+} \mathcal{F}} - \mathbf{v}\_{+} \nabla\_{\xi} \cdot \left( c\_{e} \overrightarrow{v}\_{0} \right) . \tag{4.16}$$

Dividing by ν<sup>+</sup> on both sides brings

$$\frac{\partial c\_{\epsilon}}{\partial t} = \nabla\_{\xi} \cdot \left[ D\_{\epsilon} \left( 1 - \frac{\text{d} \ln c\_{0}}{\text{d} \ln c\_{\epsilon}} \right) \nabla\_{\xi} c\_{\epsilon} \right] - \frac{\vec{\text{i}}\_{\epsilon} \nabla\_{\xi} \cdot t\_{+}^{0}}{z\_{+} \mathbf{v}\_{+} \mathcal{F}} - \nabla\_{\xi} \cdot \left( c\_{\epsilon} \vec{\text{v}}\_{0} \right), \tag{4.17}$$

85

where the identity ∇<sup>ξ</sup> · ~i *e t* 0 + <sup>=</sup> ✟✟✟ ✟✟✯<sup>0</sup> ∇ξ · ~i *e t* 0 <sup>+</sup> <sup>+</sup>~<sup>i</sup> *<sup>e</sup>*∇<sup>ξ</sup> · *t* 0 <sup>+</sup> was used. The first part of Equation (4.17) is due to the diffusion process, where *D<sup>e</sup>* is the diffusion coefficient and the driving force is the concentration gradient ∇<sup>ξ</sup> *c<sup>e</sup>* . In the second part, the migration part, is governed by the electric current density~i. Finally, the third part is due to the velocity ~*v*<sup>0</sup> of the solvent and is called convection part. It is usually assumed that dln*c*0/dln*c<sup>e</sup>* ≈ 0 and~*v*<sup>0</sup> = 0 [36, 46]. Additionally, *t* 0 <sup>+</sup> is treated as concentration independent, which was observed by experiments [130] and which makes the gradient term vanish. However, in [131], it was indicated that this may not be correct. Finally, the microscopic conservation of mass of the electrolyte phase can be reduced to

$$\frac{\partial c\_{\epsilon}}{\partial t} = \nabla\_{\xi} \cdot \left( D\_{\epsilon} \nabla\_{\xi} c\_{\epsilon} \right) \,. \tag{4.18}$$

From Section 2.2, using Equation (2.25), the macroscopic form of Equation (4.18) is

$$
\phi\_e \frac{\partial \overline{c}\_e}{\partial t} = \nabla \cdot \left( D\_{e, \text{eff}} \nabla \overline{c}\_e \right) - a\_{se} (1 - t\_+^0) \overline{j}\_{se} \,, \tag{4.19}
$$

where, comparing to Equation (2.25), φ<sup>β</sup> ≡ φ*<sup>e</sup>* is the volume fraction, *p*<sup>β</sup> ≡ *c<sup>e</sup>* is the macroscopic concentration, *k* <sup>β</sup>,eff ≡ *De*,eff is the effective diffusivity and *a*αβ ≡ *ase* is the specific surface area of the electrolyte phase. Recall that the reaction term *j se* represents the lithium flux density at the interface between solid and electrolyte phase. Note that it is expected that the anions are not taking part in the electrochemical reaction which is accounted for by the term (1−*t* 0 <sup>+</sup>). Thus, Equation (4.19) represents the macroscopic continuity equation describing mass conservation of the electrolyte phase.

#### **Reaction kinetics**

The electrochemical reaction at the interface of electrolyte and active material is described by a Butler-Volmer type equation [10, 47, 132] which, in case of a lithium-ion cell, can be written as

$$\begin{split} \overline{j}\_{se} = -j(\overline{\eta}, \overline{c}\_{e}, c\_{s, \text{surf}}) &= -\frac{i\_{0}}{\mathcal{F}} \cdot \left\{ \exp\left(\frac{(1-a)\mathcal{F}}{\mathcal{R}T}\overline{\eta}\right) - \exp\left(-\frac{a\mathcal{F}}{\mathcal{R}T}\overline{\eta}\right) \right\} \\ &= -k\_{0}\overline{c}\_{e}^{1-a} (c\_{s, \text{max}} - c\_{s, \text{surf}})^{1-a} c\_{s, \text{surf}}^{a} \\ &\cdot \left\{ \exp\left(\frac{(1-a)\mathcal{F}}{\mathcal{R}T}\overline{\eta}\right) - \exp\left(-\frac{a\mathcal{F}}{\mathcal{R}T}\overline{\eta}\right) \right\}, \end{split} \tag{4.20}$$

where *i*<sup>0</sup> is the exchange current density, *k*<sup>0</sup> is the effective reaction rate constant [10], α or (1 − α) are symmetry factors representing a favouring of cathodic or anodic reaction, respectively. Typically, α is set to 0.5 [47]. Additionally, *cs*,max is the maximum concentration in the solid phase and *cs*,surf is the concentration at the particle surface. Here,

$$\overline{\eta} = (\overline{\varphi}\_s - \overline{\varphi}\_e) - E\_{\text{eq}}(c\_{s, \text{surf}}) \tag{4.21}$$

is the overpotential indicating if the potential difference between solid and electrolyte is above or below the equilibrium potential *E*eq(*cs*,surf). The latter of which is a function of the concentration at the surface of the particle.

### **4.1.2 Classical half-cell model**

In this Section, the half-cell model is presented as it was implemented in this work using COMSOL Multiphysics [133]. Using the software's equation-based modeling module, partial differential equations (PDEs) can be entered directly. In Figure 4.2, the implemented PDEs are summarized, which are the previously presented macroscopic continuity equations.

The transport is considered through thickness direction of the cell only. Based on this assumption, transport in the lateral directions is neglected. Therefore, only one-dimensional transport is considered. Figure 4.2 shows two model levels. Each of the levels contain one-dimensional domains. These are, on the one hand, the separator and positive electrode domain, i.e. cathode, on the cell level. The separator and cathode, i.e. positive electrode, domain is denoted by the superscript "sep" and "pos", respectively. Note that the conservation of charge in the solid phase is omitted in the separator part, since there is no solid phase present. On the other hand, the one-dimensional particle domain can be found on the particle level.

On the top level, the cell level, the macroscopic continuity equations, i.e. Equations (4.6), (4.13) and (4.19), for the respective cell quantities of interest are defined. The respective cell properties of interest are the electronic potential in the solid phase ϕ*<sup>s</sup>* , the ionic potential in the electrolyte phase ϕ*<sup>e</sup>* and the electrolyte phase concentration *c<sup>e</sup>* . On the particle level, the mass conservation in the solid phase, i.e. Equation (4.2), is present. The cell property of interest is the concentration of lithium in the solid phase *c<sup>s</sup>* .

### **Boundary conditions**

In addition to the partial differential equations, the boundary conditions (BCs) are also shown in Figure 4.2. In the special case presented here, the anode is modeled as a lithium metal foil. This assumption leads to different boundary conditions compared to the full-cell model, see [36, 42]. Furthermore, galvanostatic discharge conditions are considered, thus, a constant electric current is applied. In case of galvanostatic charging processes, the sign of the applied current has to be reversed. As was mentioned before, transport in only one direction, i.e. xdirection in Figure 4.2, is considered such that volume integrals become simpler, see Appendix A.2.

### **Electronic charge in the solid phase**

As a first step, the total electronic current in the solid phase is calculated. Therefore, Equation (4.6) is integrated over the cathode volume yielding

$$\begin{aligned} A^{\text{cell}} \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\rho}\_{s}} \right) \text{d\mathbf{x}} &= A^{\text{cell}} \int\_{L^{\text{sep}}}^{L^{\text{tot}}} a\_{s\epsilon} \overline{\boldsymbol{j}}\_{s\epsilon} \mathcal{F} \mathbf{d} \mathbf{x} \\ -A^{\text{cell}} \kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\rho}\_{s}}(L^{\text{tot}}) - A^{\text{cell}} \left( -\kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\rho}\_{s}}(L^{\text{sep}}) \right) &= \\ A^{\text{cell}} \int\_{L^{\text{ep}}}^{L^{\text{tot}}} a\_{s\epsilon} \overline{\boldsymbol{j}}\_{s\epsilon} \mathcal{F} \mathbf{d} \mathbf{x} \,, \end{aligned} \tag{4.22}$$

where

$$-\kappa\_{s, \text{eff}}^{\text{enon}, \text{pos}} \nabla \overline{\varphi}\_s(L^{\text{sep}}) = 0 \tag{4.23}$$

accounts for the fact that the electronic current density must be zero at the separator-cathode interface. Moreover, the electronic current density at the cathode-current-collector interface is identified as

$$\underbrace{-\kappa\_{s,\text{eff}}^{\text{eon}}\nabla\overline{\Phi}\_{s}(L^{\text{tot}})}\_{\mathbb{T}\_{s}(L^{\text{tot}})} = \int\_{L^{\text{sep}}}^{L^{\text{tot}}} a\_{s\epsilon} \overline{j}\_{s\epsilon} \mathcal{F} \, \text{d}x \,. \tag{4.24}$$

Due to galvanostatic boundary conditions, the electronic current density at the cathode-current-collector interface is set as i *s* (*L* tot) = i app. Note that the total cell length is the sum of the separator and the positive electrode length, i.e. *L* tot = *L* sep +*L* pos .

#### **Concentration in the solid phase**

Next, the volume average change in concentration of lithium in the solid phase is calculated. Therefore, Equation (4.2) is integrated over its spherical domain and averaged over its volume, see Appendix A.2, which leads to

$$\begin{split} \frac{3}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \text{y}^{2} \frac{\partial c\_{s}}{\partial t} \, \text{dy} &= \frac{3}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \frac{\partial}{\partial \mathbf{y}} \left( \text{y}^{2} D\_{s} \frac{\partial c\_{s}}{\partial \mathbf{y}} \right) \text{dy} \\ &= \frac{3}{r\_{\text{sec}}} D\_{s} \frac{\partial c\_{s} (r\_{\text{sec}})}{\partial \mathbf{y}} - \frac{3}{r\_{\text{sec}}^{2}} \text{0}^{2} D\_{\frac{\alpha}{s}} \frac{\partial c\_{\text{sc}} (\mathbf{0})}{\partial \mathbf{y}}, \end{split} \tag{4.25}$$

where *r*sec is the radius of the spherical active material particles. Note that

$$D\_s \frac{\partial c\_s(r\_{\text{sec}})}{\partial \mathbf{y}} = \overline{j}\_{\text{se}} \tag{4.26}$$

is set to account for the electrochemical reaction boundary condition at the surfaces of the particles [40, 42], where *j se* is the Butler-Volmer type flux density from Equation (4.20). In other words, by inserting the boundary condition from Equation (4.26) into Equation (4.25), the volume-averaged change in concentration can be expressed as

$$\frac{1}{r\_{\rm sec}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^{2} \frac{\partial c\_{s}}{\partial t} \, \mathbf{dy} = \overline{j}\_{\rm se} \,. \tag{4.27}$$

#### **Ionic current and potential in the electrolyte phase**

The total ionic current in the electrolyte phase is calculated by integrating and volume-averaging Equation (4.13) over the cell volume as

$$A^{\rm cell} \int\_{0}^{L^{\rm sep}} \nabla \cdot \left( -\kappa\_{e, \rm eff}^{\rm ion, \rm sep} \nabla \overline{\rho}\_{e} - \kappa\_{D, \rm eff}^{\rm ion, \rm sep} \nabla \ln \overline{c}\_{e} \right) \mathrm{d}x + $$

$$A^{\rm cell} \int\_{L^{\rm sep}}^{L^{\rm tot}} \nabla \cdot \left( -\kappa\_{e, \rm eff}^{\rm ion, \rm pos} \nabla \overline{\rho}\_{e} - \kappa\_{D, \rm eff}^{\rm ion, \rm pos} \nabla \ln \overline{c}\_{e} \right) \mathrm{d}x = \mathbf{0} \qquad (4.28)$$
 $1 - A^{\rm cell} \int\_{L^{\rm sep}}^{L^{\rm opt}} a\_{\rm se} \overline{J}\_{\rm se} \mathcal{F} \mathbf{d}x \,.$ 

The left-hand side is

*A* cell <sup>Z</sup> *<sup>L</sup>* sep 0 ∇ · −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* −κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* d*x*+ *A* cell <sup>Z</sup> *<sup>L</sup>* tot *L* sep ∇ · −κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* −κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* d*x* = ✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ *A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* sep)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* sep) − *A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (0)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (0) + ✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✿<sup>0</sup> *A* cell −κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* tot)−κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* tot) − ✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ *A* cell −κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* sep)−κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* sep) = −*A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (0)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (0) , (4.29)

where

$$-\kappa\_{e, \text{eff}}^{\text{ion,pos}} \nabla \overline{\varphi}\_e(L^{\text{tot}}) - \kappa\_{D, \text{eff}}^{\text{ion,pos}} \nabla \ln \overline{c}\_e(L^{\text{tot}}) = 0 \tag{4.30}$$

accounts for the fact that no ionic current is allowed at the cathode-currentcollector interface. Also, continuity of ionic current at the separator-cathode interface is represented by

$$\begin{split} -\mathsf{K}\_{e,\mathrm{eff}}^{\mathrm{ion,sep}} \nabla \overline{\boldsymbol{\varphi}}\_{e}(\boldsymbol{L}^{\mathrm{sep}}) - \mathsf{K}\_{\mathrm{D,eff}}^{\mathrm{ion,sep}} \nabla \ln \overline{c}\_{e}(\boldsymbol{L}^{\mathrm{sep}}) &= \\ -\mathsf{K}\_{e,\mathrm{eff}}^{\mathrm{ion,pos}} \nabla \overline{\boldsymbol{\varphi}}\_{e}(\boldsymbol{L}^{\mathrm{sep}}) - \mathsf{K}\_{\mathrm{D,eff}}^{\mathrm{ion,pos}} \nabla \ln \overline{c}\_{e}(\boldsymbol{L}^{\mathrm{sep}}) \,. \end{split} \tag{4.31}$$

Finally, using the results of the left-hand side simplifications, leads to

$$-\underbrace{(-\mathop{\kappa}\_{e,\text{eff}}^{\text{ion,sep}}\nabla\Phi\_e(0) - \mathop{\kappa}\_{D,\text{eff}}^{\text{ion,sep}}\nabla\ln\mathbb{Z}\_e(0)}\_{\mathbb{T}\_e(0)} = -\int\_{L^{\text{sep}}}^{L^{\text{tot}}} a\_{se}\overline{\mathbb{J}}\_{se}\mathcal{F}d\mathbf{x},\tag{4.32}$$

where the total ionic current density at the anode-separator interface is expressed as

$$\overline{\mathbf{i}}\_{\epsilon}(0) = \int\_{L^{\text{sep}}}^{L^{\text{tot}}} a\_{\epsilon e} \overline{\mathbf{j}}\_{\epsilon e} \mathcal{F} \, \text{d}x \,. \tag{4.33}$$

Due to galvanostatic conditions, and by comparing Equation (4.33) to Equation (4.24), the ionic current density at the anode-separator interface is i *e* (0) = i app.

Concerning the ionic potential, usually, the focus is on the potential difference, thus, ϕ*<sup>e</sup>* is arbitrarily set to zero at the anode side [40, 42], which is expressed via

$$
\overline{\varphi}\_e(0) = 0.\tag{4.34}
$$

### **Concentration in the electrolyte phase**

The total flux in the electrolyte phase is achieved by integrating Equation (4.19) and averaging over the cell volume, which results in

$$\begin{split} \frac{1}{L^{\rm tot}} \left( \int\_{0}^{L^{\rm sep}} \theta\_{e}^{\rm sep} \frac{\partial \overline{c}\_{e}}{\partial t} \, \mathrm{d}x + \int\_{L^{\rm sep}}^{L^{\rm pos}} \theta\_{e}^{\rm pos} \frac{\partial \overline{c}\_{e}}{\partial t} \, \mathrm{d}x \right) \\ = \frac{1}{L^{\rm tot}} \left( \int\_{0}^{L^{\rm sep}} \nabla \cdot \left( D\_{e, \rm eff}^{\rm sep} \nabla \overline{c}\_{e} \right) \, \mathrm{d}x \\ \qquad + \int\_{L^{\rm sep}}^{L^{\rm tot}} \nabla \cdot \left( D\_{e, \rm eff}^{\rm pos} \nabla \overline{c}\_{e} \right) - a\_{\rm se} (1 - t\_{+}^{0}) \overline{J}\_{\rm se} \, \mathrm{d}x \right) \\ = \frac{1}{L^{\rm tot}} \left( \int\_{0}^{L^{\rm sep}} \nabla \cdot \left( D\_{e, \rm eff}^{\rm sep} \nabla \overline{c}\_{e} \right) \, \mathrm{d}x + \int\_{L^{\rm sep}}^{L^{\rm tot}} \nabla \cdot \left( D\_{e, \rm eff}^{\rm pos} \nabla \overline{c}\_{e} \right) \, \mathrm{d}x \right) \\ - \frac{1}{L^{\rm tot}} \int\_{L^{\rm sep}}^{L^{\rm sep}} a\_{\rm se} (1 - t\_{+}^{0}) \overline{J}\_{\rm se} \, \mathrm{d}x . \end{split} \tag{4.35}$$

The divergence terms give

$$\begin{split} &\int\_{0}^{L^{\rm sep}} \nabla \cdot \left( D\_{e,\rm eff}^{\rm sep} \nabla \overline{c}\_{e} \right) \mathrm{d}x + \int\_{L^{\rm sep}}^{L^{\rm iso}} \nabla \cdot \left( D\_{e,\rm eff}^{\rm pos} \nabla \overline{c}\_{e} \right) \mathrm{d}x \\ &= D\_{e,\rm eff}^{\rm sep} \nabla \overline{c}\_{e} \langle \mathsf{L}^{\rm sep} \overline{\rangle} - D\_{e,\rm eff}^{\rm sep} \nabla \overline{c}\_{e} \left( 0 \right) + D\_{e,\rm eff}^{\rm pos} \nabla \overline{c}\_{e} \langle \mathsf{L}^{\rm sep} \rangle \overline{\prescript{\rm pos}{}{}{\left( \mathrm{L}^{\rm sep} \right)}} - D\_{e,\rm eff}^{\rm pos} \nabla \overline{c}\_{e} \langle \mathsf{L}^{\rm sep} \rangle \\ &= -D\_{e,\rm eff}^{\rm sep} \nabla \overline{c}\_{e} \left( 0 \right), \end{split} \tag{4.36}$$

where

$$D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(L^{\text{sep}}) = D\_{e, \text{eff}}^{\text{pos}} \nabla \overline{c}\_e(L^{\text{sep}}) \quad \text{and} \quad D\_{e, \text{eff}}^{\text{pos}} \nabla \overline{c}\_e(L^{\text{tot}}) = 0 \tag{4.37}$$

reflects continuity of mass flux between separator and cathode region and implies that no mass flux is allowed at the cathode-current-collector side, respectively.

Using Equation (4.36) in Equation (4.35) yields

$$\begin{split} \frac{1}{L^{\mathrm{tot}}} \left( \int\_{0}^{L^{\mathrm{sep}}} \phi\_{\epsilon}^{\mathrm{sep}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \mathrm{d}x + \int\_{L^{\mathrm{sep}}}^{L^{\mathrm{tot}}} \phi\_{\epsilon}^{\mathrm{pos}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \mathrm{d}x \right) &= \\ -\frac{1}{L^{\mathrm{tot}}} \left( D\_{\epsilon, \mathrm{eff}}^{\mathrm{sep}} \nabla \overline{c}\_{\epsilon}(0) \right) - \frac{1}{L^{\mathrm{tot}}} \left( \int\_{L^{\mathrm{sep}}}^{L^{\mathrm{tot}}} a\_{\epsilon \epsilon} (1 - t\_{+}^{0}) \overline{j}\_{s \epsilon} \, \mathrm{d}x \right) . \end{split} \tag{4.38}$$

Equation (4.38) shows that the volume-averaged temporal change of the electrolyte concentration is due to diffusion at the anode-separator interface and the electrochemical reactions of cations at the electrolyte-solid interface. As was assumed earlier, the applied electroneutrality condition guarantees the same amount of anions and cations in a region. Globally, this is expressed as a constant average electrolyte concentration, i.e. a vanishing temporal change. Therefore, the left-hand side of Equation (4.38) is set to zero yielding

$$-D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(0) = \int\_{L^{\text{sep}}}^{L^{\text{tot}}} a\_{se}(1 - t\_+^0) \overline{j}\_{se} \, \text{d}x \,. \tag{4.39}$$

In other words, the flux of cations due to electrochemical reactions at the surface area of electrolyte and active material is compensated by migration at the anodeseparator interface [38]. Comparing Equation (4.33) to Equation (4.39), the righthand side of the latter can be rewritten in terms of the applied current density as

$$\int\_{L^{\rm sep}}^{L^{\rm tot}} a\_{\rm se}(1 - t\_+^0) \overline{j}\_{\rm se} \, \mathrm{d}x = \frac{(1 - t\_+^0)}{\mathcal{F}} \underbrace{\int\_{L^{\rm sep}}^{L^{\rm tot}} a\_{\rm se} \overline{j}\_{\rm se} \mathcal{F} \, \mathrm{d}x}\_{\overline{i}\_\epsilon(0) = i\_{\rm app}} = \frac{\mathrm{i}\_{\rm app}}{\mathcal{F}} (1 - t\_+^0), \tag{4.40}$$

where, following the assumptions made in Section 2.2, the specific surface area *ase* and, additionally, *t* 0 <sup>+</sup> are assumed constant and, therefore, can be extracted from the integral. Thus, the boundary condition for the flux in the electrolyte phase at the anode-separator interface can be written as

$$-D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(0) = \frac{\mathbf{i}\_{\text{app}}}{\mathcal{F}} (1 - t\_+^0) \,. \tag{4.41}$$

The boundary condition in Equation (4.41) is commonly used to replace the anode in a half-cell model [40, 42, 132]. However, disassembling the continuity equations for the electrolyte phase, as was done above, shows the implications that come along with it. Such as, the volume average lithium concentration in the electrolyte can be assumed to stay constant over time.

### **Cell quantities**

The half-cell model, as presented above, is used to evaluate cell performance. To this end, typically, cell quantities like state of charge, specific capacity, cell voltage and alike are computed and discussed.

### **Specific capacity**

The total capacity, or ampere-hour capacity, or electric capacity, of the cell can refer to either the negative or positive electrode. It is calculated via the initial and maximum lithium concentration in the solid phase of the respective electrode. In case of the half-cell model used in this work, the total capacity of the positive electrode is given as

$$\begin{split} \mathcal{Q}\_{\text{tot}} &= V^{\text{pos}} (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F} \\ &= A^{\text{cell}} L^{\text{pos}} \phi\_s (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}, \end{split} \tag{4.42}$$

where *V* pos *<sup>s</sup>* and φ*<sup>s</sup>* is the volume fraction of the solid phase of the positive electrode, respectively. *A* cell and *L* pos is the volume, cross section area and length of the positive electrode, respectively. *cs*,init and *cs*,max is the initial and maximum allowable concentration of the active material, respectively. The specific capacity is commonly used in the battery community. For this purpose, the capacity is related to the mass of the active material, i.e. *m*AM, which brings

$$\begin{split} \mathcal{Q}\_{\text{spec}} &= \frac{\mathcal{Q}\_{\text{tot}}}{m\_{\text{AM}}} = \frac{A^{\text{cell}} L^{\text{pos}} \phi\_s (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}}{A^{\text{cell}} L^{\text{pos}} \phi\_s \, \mathcal{Q}\_{\text{AM}}} \\ &= \frac{(c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}}{\mathcal{Q}\_{\text{AM}}} \,, \end{split} \tag{4.43}$$

where %AM is the density of the underlying active material.

### **Cell voltage**

The cell voltage is defined as the voltage drop between positive and negative electrode-current-collector interface. Since the anode is represented as a lithium metal foil, the anode potential is assumed to be zero. Thus, the cell voltage is equal to the cathode-current-collector interface as

$$U\_{\rm cell} = \overline{\varphi}\_s(L^{\rm tot})\,. \tag{4.44}$$

#### **State of charge and depth of discharge**

Typically, the cell voltage is drawn versus the so-called state of charge (SOC) or depth of discharge (DOD) of the cell. SOC is defined as the releasable capacity relative to its maximum allowable, whereas DOD is defined as the released capacity relative to its maximum allowable [134, 135]. In case of a full-cell setup, the state of charge or depth of discharge can be related to the capacities of either the positive or the negative electrode. As for the half-cell model in this work, the state parameters are defined using the positive electrode capacity. Since the electrode capacity is directly related to the amount of lithium stored inside the active material, the DOD can be calculated as

$$\text{DOD}(t) = \frac{\overline{c}\_{s, \text{avg}}(t) - c\_{s, \text{init}}}{c\_{s, \text{max}} - c\_{s, \text{init}}},\tag{4.45}$$

where *cs*,avg(*t*) is the cell's volume-averaged lithium concentration of the solid phase at time *t*. During discharge, as the lithium concentration increases in the positive electrode, the value of Equation (4.45) increases from 0 to 1, representing an initial concentration state and a maximum concentration state where, ideally, all sites of the active material crystal structure are filled with lithium [10]. In other words, the cell's depth of discharge increases while, by the same time, its state of charge decreases from 1 to 0. Therefore, the cell's SOC and DOD are related via

$$\text{SOC} = 1 - \text{DOD} \,. \tag{4.46}$$

### **Calculating the** *nC* **current**

In the battery community, it is very common to use the concept of C rate. The C rate represents an electrical current, where a C rate of 1 corresponds to a current that is needed to charge or discharge a LIB within one hour. Higher C rates stand for higher currents whereas lower C rates mean lower currents. In reality, the value of the C rate is highly dependent on the structure, geometry, electrochemistry, etc., such that, usually, experiments are needed, in order to find the electric current representing a C rate of 1.

In the cell model, however, the 1*C* rate representing electric current can be calculated directly from the given structural parameters. It was shown that the state of charge can be calculated by knowing the state of lithium concentration of the electrode. In case of the presented half-cell model, the volume-averaged temporal change of lithium concentration of the positive electrode is generally calculated as

$$\frac{\partial \overline{c}\_{s, \text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \frac{\mathfrak{J}}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} \frac{\partial c\_{s}}{\partial t} \, \text{dy} \, \text{d}x. \tag{4.47}$$

Recall that, from Equation (4.25), the temporal change of lithium of a single particle is known as

$$\frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 \frac{\partial c\_s}{\partial t} \, \mathbf{dy} = \frac{3}{r\_{\rm sec}} \overline{j}\_{\rm se} \,. \tag{4.48}$$

Volume-averaged integration of Equation (4.48) over the electrode volume leads to

$$\frac{1}{L^{\rm pos}} \int\_{L^{\rm sep}}^{L^{\rm tot}} \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 \frac{\partial c\_s}{\partial t} \mathbf{d} \mathbf{y} \, \mathbf{dx} = \frac{1}{L^{\rm pos}} \int\_{L^{\rm sep}}^{L^{\rm tot}} \frac{3}{r\_{\rm sec}} \overline{j}\_{s\varepsilon} \, \mathbf{dx} . \tag{4.49}$$

Inserting Equation (4.49) into Equation (4.47) brings

$$\frac{\partial \overline{\sigma}\_{s, \text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \frac{\mathfrak{Z}}{r\_{\text{sec}}} \overline{j}\_{s\epsilon} \, \text{d} \mathbf{x} \,. \tag{4.50}$$

Additionally, from Equation (4.24) and the boundary condition iapp = i *s* (*L* tot), the applied current density is represented as

$$\mathbf{i}\_{\rm app} = \int\_{L^{\rm sep}}^{L^{\rm tot}} a\_{\rm se} \overline{\mathbf{j}}\_{\rm se} \mathcal{F} \, \mathrm{d}x \,. \tag{4.51}$$

Under the assumption that the interfacial area *ase* and the active material particle radii *r*sec are constant, they can be extracted form the integrals of Equations (4.50) and (4.51), converting them into

$$\frac{\partial \overline{c}\_{s, \text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \frac{\mathfrak{J}}{r\_{\text{sec}}} \int\_{L^{\text{exp}}}^{L^{\text{tot}}} \overline{j}\_{\text{se}} \, \text{d}x \tag{4.52}$$

and

$$\mathbf{i}\_{\rm app} = a\_{\rm se} \mathcal{F} \int\_{L^{\rm sep}}^{L^{\rm tot}} \mathbf{\overline{j}}\_{\rm se} \, \mathrm{d}x \,. \tag{4.53}$$

Finally, comparing Equations (4.52) and (4.53), brings

$$\frac{\partial \overline{c}\_{s, \text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \frac{\mathfrak{Z}}{r\_{\text{sec}}} \frac{\mathfrak{i}\_{\text{app}}}{a\_{\text{se}} \mathcal{F}} \,. \tag{4.54}$$

From Equation (4.54), it can be seen that the volume-averaged temporal change of lithium concentration of the positive electrode can be computed by the applied current density iapp and the structural properties *L* pos , *r*sec and *ase*. Note that those properties are taken as constants in this work, which also makes Equation (4.54) constant.

In order to calculate *cs*,avg(*t*) at time *t*, Equation (4.54) is multiplied by *t* and the initial concentration *cs*,init is added such that

$$\begin{split} \overline{c}\_{s,\text{avg}}(t) &= \frac{\partial \overline{c}\_{s,\text{avg}}}{\partial t} \cdot t + c\_{s,\text{init}} \\ &= \frac{1}{L^{\text{pos}}} \frac{\mathfrak{Z}}{r\_{\text{sec}}} \frac{\mathrm{i}\_{\text{app}}}{a\_{\text{sec}} \mathcal{F}} \cdot t + c\_{s,\text{init}} \cdot \mathrm{d} \end{split} \tag{4.55}$$

Inserting Equation (4.55) in Equation (4.45), the depth of discharge at time *t* is

$$\text{DOD}(t) = \frac{\frac{1}{D^{\text{src}}} \frac{3}{r\_{\text{src}}} \frac{\text{i}\_{\text{app}}}{a\_{\text{se}} \mathcal{F}} \cdot t}{c\_{s, \text{max}} - c\_{s, \text{init}}}. \tag{4.56}$$

A fully discharged cell is represented by setting DOD(*t*) = 1 within one hour *t* = 3600 s in Equation (4.56). Rearranging yields the 1*C*-current as

$$\mathbf{i}\_{\rm app}^{\rm IC} = (c\_{s,\rm max} - c\_{s,\rm init}) \frac{L^{\rm pos} r\_{\rm sec} a\_{\rm se} \mathcal{F}}{\mathbf{3} \cdot \mathbf{3} \mathbf{6} \mathbf{0} \mathbf{0}}.\tag{4.57}$$

Generally, the *nC*-current is

$$\mathbf{i}\_{\rm app}^{\rm nC} = n(c\_{s,\rm max} - c\_{s,\rm init}) \frac{L^{\rm pos} \, r\_{\rm sec} \, a\_{\rm se} \, \mathcal{F}}{\mathfrak{3} \cdot \mathfrak{3}600} \, . \tag{4.58}$$

Note that the above derivations do only apply for cases where the specific surface area can be calculated using *ase* = 3φ*<sup>s</sup> r*sec , which implies the assumption that the active material particles are of spherical shape and detached from each other [136]. Morphological variations from this formula necessitate special treatment of the above presented derivations, which shall not be considered in the present work.

Figure 4.2: Mathematical model of the half-cell setup.

# **4.2 Hierarchically structured half-cell**

As mentioned in Section 1, the performance of LIBs can be increased by socalled hierarchically structured cathodes where, by postprocessing the initial active material, the porosity of secondary particles is increased. As a result, the rate capability and cycle stability can be improved [26–29]. In order to simulate such structures, in the following, the previously presented half-cell model is extended for hierarchically structured cathodes. Consequently, the model is called hierarchically structured half-cell model [137].

In Figure 4.3a, the sketch of the hierarchically structured half-cell is presented. Similar to the classical half-cell setup from Section 4.1, on the left-hand side, the anode is a lithium metal foil and, on the right-hand side, the positive electrode is a porous composite structure which is composed of the active material and the filler material. The latter of which is a carbon-black-binder-mixture. Also, a separator is placed between anode and positive electrode. In the hierarchically structured electrode case, however, the active material secondary particles which are built up by primary particles—show a distinct porosity. Both the separator and composite electrode structure—including the secondary particle pores—are filled with liquid electrolyte.

During charge or discharge, electrochemical reactions take place only at the interfacial area between primary particles and electrolyte inside the secondary particles. Due to the porous secondary particles, an additional transport region is introduced into the system. Therefore, as an extension of the classical halfcell model [36, 40–42], the hierarchically structured half-cell model is divided into three levels. The cell, the secondary particle and the primary particle level. Following the classical half-cell model, the porous electrode theory [35, 47] does also apply to the hierarchically structured half-cell model.

Figure 4.3: Hierarchically structured half-cell setup of lithium-ion batteries. a) Sandwich structure of a lithium-ion battery cell. The anode is a solid lithium metal and the cathode has porous active materials, i.e. secondary particles. b) The porous structure of both the cathode and the secondary particles are homogenized.

On cell level, transport is carried by the electrolyte and the solid phase, where, similar to the half-cell model, the solid phase accounts for the active material and the filler material. On secondary particle level, all transport is via the electrolyte or the primary particles, where the latter of which is identified as the solid phase on this level. Finally, the intercalation and transport process of lithium inside the solid phase is modeled on the primary particle level.

A model for porous secondary particles was also presented in [53]. In the present work, however, the secondary particle and electrolyte share the same specific surface area for exchange flux. Additionally, to avoid double counting of exchange flux density, the electrochemical reactions are only between primary particles and the surrounding electrolyte within the secondary particle and electrochemical reactions between secondary particle and electrolyte are neglected. Moreover, transport of lithium inside the secondary particles is only via the internal electrolyte phase and is intercalated via electrochemical reactions into the additionally introduced primary particle level.

### **4.2.1 Macroscale equations**

In this section, the mathematical formulation of the hierarchically structured half-cell model is derived. The electronic and ionic charge transport, as well as cationic flux has to be accounted for on both, the cell and the secondary particle level. Therefore, in total, six continuity equations have to be defined and solved. Similar to Section 4.1.1, microscopic transport equations are defined for the respective phases and volume average theorems are used to convert them to macroscopic forms. Additionally, on primary particle level, the lithium transport is modeled by a Fickian type diffusion. A summarized version of the presented PDEs can be found in Figure 4.4 and 4.5. In the following, the superscript "sec" and "prim" refers to the considered level, i.e. secondary particle or primary particle level. If no superscript is given, the properties refer to the cell level.

### **Primary particle level**

On the primary particle level, lithium moves due to diffusion. Spherical symmetry is assumed for the primary particles.

### **Conservation of mass in the solid phase**

The transport process inside the primary particles is taken as one-dimensional Fickian diffusion by

$$\frac{\partial c\_s^{\text{prim}}}{\partial t} = \frac{1}{z^2} \frac{\partial}{\partial z} \left( z^2 D\_s \frac{\partial c\_s^{\text{prim}}}{\partial z} \right) . \tag{4.59}$$

Here, the concentration of lithium is *c* prim *<sup>s</sup>* and the diffusion coefficient in the solid phase is *D<sup>s</sup>* . The dimension *z* is introduced as the radial coordinate inside the primary particle.

### **Secondary particle level**

On the secondary particle level, the derivation of the macroscale forms of the conservation equations is quite similar to the derivations in Section 4.1.1. Additionally, and as can be observed by experiments [29], hierarchically structured secondary particles are of spherical shape. Therefore, spherical symmetry is employed and the derived macroscale equations are converted into spherically symmetric forms according to Appendix A.1.

### **Conservation of charge in the solid phase**

Similar to Equation (4.6), the macroscopic form of the conservation of charge in the solid phase of a secondary particle is

$$\nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{coon}, \text{sec}} \nabla \overline{\mathfrak{P}}\_s^{\text{sec}} \right) = a\_{\text{se}}^{\text{sec}} \overline{J}\_{\text{se}}^{\text{sec}} \mathcal{F} \,, \tag{4.60}$$

where κ eon,sec *<sup>s</sup>*,eff is the effective electronic conductivity of the solid phase, ϕ sec *s* is the macroscopic electronic potential and *a* sec *se* is the specific interfacial area of the primary particle surfaces and the electrolyte inside the secondary particles. Note that the macroscopic surface flux of the secondary particles *j* sec *se* can be calculated using the Butler-Volmer relation according to Equation (4.69).

Converting to the spherical coordinate system, and employing spherical symmetry, yields

$$\frac{1}{\mathcal{Y}^2} \frac{\partial}{\partial \mathbf{y}} \left( -\mathbf{y}^2 \kappa\_{s, \text{eff}}^{\text{con,sec}} \frac{\partial \overline{\mathbf{g}}\_s^{\text{sec}}}{\partial \mathbf{y}} \right) = a\_{\text{se}}^{\text{sec}} \overline{\mathbf{j}}\_{s\epsilon}^{\text{sec}} \mathcal{F}. \tag{4.61}$$

Here, the dimension *y* is introduced as the radial coordinate inside the secondary particle.

### **Conservation of charge in the electrolyte phase**

Analogous to Equation (4.13), the macroscopic form of the conservation of charge in the electrolyte phase of a secondary particle is

$$\nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion,sec}} \nabla \overline{\boldsymbol{\varphi}}\_{e}^{\text{sec}} - \kappa\_{D, \text{eff}}^{\text{ion,sec}} \nabla \ln \overline{c}\_{e}^{\text{sec}} \right) = -a\_{\text{se}}^{\text{sec}} \overline{j}\_{\text{se}}^{\text{sec}} \mathcal{F}, \tag{4.62}$$

where κ ion,sec *<sup>e</sup>*,eff and κ ion,sec *<sup>D</sup>*,eff is the effective ionic and diffusional conductivity of the electrolyte phase, respectively, ϕ sec *e* is the macroscopic ionic potential and *c* sec *e* is the macroscopic electrolyte concentration inside the secondary particles.

The spherically symmetric version of Equation (4.62) is

$$\frac{1}{\mathcal{Y}^2} \frac{\partial}{\partial \mathbf{y}} \left( -\mathbf{y}^2 \kappa\_{e, \text{eff}}^{\text{ion,sec}} \frac{\partial \overline{\mathbf{g}}\_e^{\text{sec}}}{\partial \mathbf{y}} - \mathbf{y}^2 \kappa\_{D, \text{eff}}^{\text{ion,sec}} \frac{\partial \ln \overline{\mathbf{c}}\_e^{\text{sec}}}{\partial \mathbf{y}} \right) = -a\_{\text{se}}^{\text{sec}} \overline{\mathbf{J}}\_{\text{se}}^{\text{sec}} \mathcal{F}. \tag{4.63}$$

#### **Conservation of mass in the electrolyte phase**

Finally, on the secondary particle level, the conservation of mass in the electrolyte is derived as follows. Adapting the macroscopic form of the conservation of mass in the electrolyte phase from Equation (4.19) as

$$\phi\_e^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sec}}{\partial t} = \nabla \cdot \left( D\_{e, \rm eff}^{\rm sec} \nabla \overline{c}\_e^{\rm sec} \right) - a\_{\rm se}^{\rm sec} (1 - t\_+^0) \overline{j}\_{\rm se}^{\rm sec}, \tag{4.64}$$

where *D* sec *<sup>e</sup>*,eff is the effective diffusion coefficient of the electrolyte phase inside the secondary particle.

Converting Equation (4.64) to the spherical coordinate system and employing spherical symmetry yields

$$\phi\_e^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sec}}{\partial t} = \frac{1}{\mathbf{y}^2} \frac{\partial}{\partial \mathbf{y}} \left( \mathbf{y}^2 D\_{e, \rm eff}^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sec}}{\partial \mathbf{y}} \right) - a\_{\rm se}^{\rm sec} (1 - t\_+^0) \overline{f}\_{\rm se}^{\rm sec} \,. \tag{4.65}$$

### **Cell level**

In order to derive macroscopic continuity equations for the cell level of the hierarchically structured half-cell model, volume averaging methods are applied according to Section 2.2.

### **Charge conservation in the solid phase**

The microscopic continuity equation regarding the conservation of electronic charge on the cell level is analogous to Equation (4.4). Using the volume average approach by employing Equation (2.24) brings

$$\nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{enon}} \nabla \overline{\varphi}\_s \right) = \phi\_{\text{sec}} \overline{\mathbf{i}}\_{s, \text{sec}} \,, \tag{4.66}$$

where *k*α,eff ≡ κ eon *<sup>s</sup>*,eff is the effective electronic conductivity, *p*<sup>α</sup> ≡ ϕ*<sup>s</sup>* is the macroscopic electronic potential. In contrast to the common half-cell model, the electrochemical reactions at the interface between the secondary particle surfaces and the surrounding electrolyte, i.e. the surface term in Equation (2.24), is neglected in this work. Rather, the production of electronic charges on the cell level is due to electrochemical reactions from within the secondary particles on the surface area of the primary particles. This is expressed as the volumeaveraged source term being identified as φα*b*<sup>α</sup> ≡ φseci*s*,sec, where φsec is the volume fraction of secondary particles without internal porosity and i*s*,sec is the volume-averaged production term of electronic charges inside the secondary particles.

### **Charge conservation in the electrolyte phase**

The microscopic continuity equation regarding the conservation of ionic charge is identical to Equation (4.10). Using the volume-average approach, by applying Equation (2.25), yields

$$\nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion}} \nabla \overline{\otimes}\_e - \kappa\_{D, \text{eff}}^{\text{ion}} \nabla \ln \overline{c}\_e \right) = \phi\_{\text{sec}} \overline{\mathbf{i}}\_{e, \text{sec}} \,. \tag{4.67}$$

Comparing to Equation (2.25), the effective conductivity, *k* <sup>β</sup>,eff, is either the effective ionic conductivity, κ ion *<sup>e</sup>*,eff, or the effective diffusional conductivity κ ion *D*,eff on the cell level. Moreover, the macroscopic potential, *p*<sup>β</sup> , is either the macroscopic ionic potential, ϕ*<sup>e</sup>* , or, *c<sup>e</sup>* , the macroscopic concentration of lithium of the electrolyte phase. Additionally, the generation of ionic charges is due to electrochemical reactions on the surfaces of the primary particles inside the secondary particles. Therefore, the volumetric production term on the right-hand side of Equation (4.67) is defined by the volume-averaged production term of ionic charges from within the secondary particles i*e*,sec. Note that the reaction term in Equation (2.25) on the surfaces of the secondary particles is neglected.

#### **Mass conservation in the electrolyte phase**

The derivation of the macroscopic continuity equation regarding the mass conservation of the electrolyte phase starts with the microscopic form provided by Equation (4.18). Using volume-averageing, by employing Equation (2.25), results in

$$
\phi\_e \frac{\partial \overline{c}\_e}{\partial t} = \nabla \cdot \left( D\_{e, \text{eff}} \nabla \overline{c}\_e \right) + \phi\_{\text{sec}} \overline{J}\_{e, \text{sec}} \,. \tag{4.68}
$$

where, comparing to Equation (2.25), φ<sup>β</sup> ≡ φ*<sup>e</sup>* is the volume fraction, *p*<sup>β</sup> ≡ *c<sup>e</sup>* is the macroscopic concentration and *k* <sup>β</sup>,eff ≡ *De*,eff is the effective diffusion coefficient on cell level. The electrochemical reaction on the surfaces between the secondary particles and the electrolyte is neglected. Instead, the production of cations *je*,sec is governed by electrochemical reaction from within the secondary particles.

### **Reaction kinetics**

Similar to Section 4.1.1, the electrochemical reaction at the interface between primary particles and electrolyte is described by a Butler-Volmer type equation on secondary particle level. In case of the hierarchically structured half-cell model, it reads

$$\begin{split} \overline{f}\_{sc}^{\text{sec}} &= -j(\overline{\eta}^{\text{sec}}, \overline{c}\_{e}^{\text{sec}}, c\_{s, \text{surf}}^{\text{sec}}) \\ &= -\frac{\overline{\overline{\eta}}\_{0}}{\mathcal{F}} \cdot \left\{ \exp\left(\frac{(1-\alpha)\mathcal{F}}{\mathcal{R}T} \overline{\eta}^{\text{sec}}\right) - \exp\left(-\frac{\alpha \mathcal{F}}{\mathcal{R}T} \overline{\eta}^{\text{sec}}\right) \right\} \\ &= -k\_{0} (\overline{c}\_{e}^{\text{sec}})^{1-\alpha} (c\_{s, \text{max}} - c\_{s, \text{surf}}^{\text{sec}})^{1-\alpha} (c\_{s, \text{surf}}^{\text{sec}})^{\alpha} \\ &\cdot \left\{ \exp\left(\frac{(1-\alpha)\mathcal{F}}{\mathcal{R}T} \overline{\eta}^{\text{sec}}\right) - \exp\left(-\frac{\alpha \mathcal{F}}{\mathcal{R}T} \overline{\eta}^{\text{sec}}\right) \right\}, \end{split} \tag{4.69}$$

where, compared to Equation (4.69), i sec 0 is the exchange current density on secondary particle level and *c* sec *<sup>s</sup>*,surf is the lithium concentration at the primary particle surfaces. Moreover, η sec refers to the overpotential at the surface of a primary particle inside a secondary particle by

$$
\overline{\eta}^{\rm sec} = (\overline{\rho}\_s^{\rm sec} - \overline{\rho}\_e^{\rm sec}) - E\_{\rm eq}(c\_{s,\rm surf}^{\rm sec}) \,, \tag{4.70}
$$

where the equilibrium potential, *E*eq(*c* sec *<sup>s</sup>*,surf), is a function of *c* sec *<sup>s</sup>*,surf.

### **4.2.2 Hierarchically structured half-cell model**

As can be seen in Figures 4.4 and 4.5, the hierarchically structured half-cell is divided into three levels, each of which containing one-dimensional domains. On the top level, the cell level, the separator and the positive electrode domain, i.e. cathode, is modeled. Here, the superscripts "sep" and "pos" indicate the respective region.

### **Boundary conditions**

Similar to the half-cell model from Section 4.1.2, the anode is modeled as a lithium metal foil. Again, galvanostatic charge or discharge is considered, thus the applied current is assumed to be constant. In the following, boundary conditions are derived by integrating over domains in the cartesian and spherical coordinate systems. The mathematical basics of the integration operations performed can be found in Appendix A.2.

### **Electronic charge and potential in the solid phase**

As a first step, concerning the right-hand side of Equation (4.66), the volumeaveraged source term i*s*,sec is evaluated by volume-averaged integration of Equation (4.61), yielding

$$\begin{split} & \frac{3}{r\_{\rm sec}} \int\_{0}^{r\_{\rm sec}} \frac{\partial}{\partial \mathbf{y}} \left( -\mathbf{y}^{2} \, \mathbf{x}\_{s, \rm eff}^{\rm cons,sec} \frac{\partial \overline{\Phi}\_{s}^{\rm sc}}{\partial \mathbf{y}} \right) \mathbf{d} \mathbf{y} = \frac{3}{r\_{\rm sc}^{-3}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^{2} \, \mathbf{a}\_{\rm sc}^{\rm sec} \overline{f}\_{\rm sc}^{\rm sc} \mathcal{F} \mathbf{d} \mathbf{y} \\ & \qquad \qquad \qquad \qquad \qquad \frac{3}{r\_{\rm sec}} r\_{\rm sec} \, ^{2} \mathbf{x}\_{s, \rm eff}^{\rm cons,sec} \frac{\partial \overline{\Phi}\_{s}^{\rm sc} (r\_{\rm sec})}{\partial \mathbf{y}} - \\ & \qquad \qquad \qquad \qquad \left( -\frac{3}{r\_{\rm sec}} 0^{2} \, \mathbf{x}\_{\rm r, \rm eff}^{\rm cons,sec} \frac{\partial \overline{\Phi}\_{s}^{\rm sc} (\mathbf{0})}{\partial \mathbf{y}} \right) = \frac{3}{r\_{\rm sec}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^{2} \, \mathbf{a}\_{\rm sc}^{\rm sec} \overline{f}\_{\rm sc}^{\rm sec} \mathcal{F} \mathbf{d} \mathbf{y}, \end{split} \tag{4.71}$$

where *r*sec is the secondary particles radius. Consequently, by recognizing that the electronic current density at the surface of the secondary particle is

$$\underbrace{-\mathsf{K}\_{s,\text{eff}}^{\text{con,sec}}\frac{\partial\overline{\mathsf{op}}\_{s}^{\text{sec}}(r\_{\text{sec}})}{\partial\mathbf{y}}}\_{\overline{\mathsf{i}}\_{s,\text{surf}}^{\text{sec}}}=\frac{1}{r\_{\text{sec}}^{2}}\int\_{0}^{r\_{\text{sec}}}\mathcal{y}^{2}a\_{\text{se}}^{\text{sec}}\overline{\mathcal{J}}\_{\text{se}}^{\text{sec}}\mathcal{F}\,\mathrm{d}\mathbf{y},\tag{4.72}$$

i*s*,sec can be expressed as

$$\overline{\mathbf{i}}\_{s, \text{sec}} = \frac{\mathfrak{Z}}{r\_{\text{sec}}} \overline{\mathbf{i}}\_{s, \text{surf}}^{\text{sec}} = \frac{\mathfrak{Z}}{r\_{\text{sec}}} \int\_0^{r\_{\text{sec}}} \mathbf{y}^2 a\_{\text{sec}}^{\text{sec}} \overline{\mathbf{j}}\_{\text{se}}^{\text{sec}} \mathcal{F} \mathbf{d} \mathbf{y}, \tag{4.73}$$

109

which, in turn, can be used for the right-hand side of Equation (4.66) yielding

$$\boldsymbol{\phi}\_{\text{sec}}\mathbf{\tilde{i}}\_{s,\text{sec}} = \underbrace{\boldsymbol{\phi}\_{\text{sec}}\,\mathbf{\tilde{i}}\_{s,\text{surf}}}\_{a\_{st}}\mathbf{\tilde{i}}\_{s,\text{surf}}^{\text{sec}} = a\_{se}\mathbf{\tilde{i}}\_{s,\text{surf}}^{\text{sec}},\tag{4.74}$$

where, in case of equal-sized and detached spherical secondary particles, *ase* = 3φsec *r*sec [136] is identified as specific surface area between secondary particle and electrolyte.

As the next step, the electronic current over the cell domain is calculated by volume-averaged integration of Equation (4.66) as

$$\begin{aligned} A^{\text{cell}} \int\_{L^{\text{sp}}}^{L^{\text{tot}}} \nabla \cdot \left( -\kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\varphi}}\_{s} \right) \text{d\mathbf{x}} &= A^{\text{cell}} \int\_{L^{\text{sp}}}^{L^{\text{tot}}} \boldsymbol{\phi}\_{\text{sec}} \overline{\boldsymbol{\mathfrak{i}}}\_{s, \text{sec}} \, \text{d\mathbf{x}} \\ -A^{\text{cell}} \kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\mathfrak{p}}}\_{s} (L^{\text{tot}}) - A^{\text{cell}} \left( -\kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\boldsymbol{\mathfrak{p}}}\_{s} (L^{\text{sep}}) \right) &= \\ A^{\text{cell}} \int\_{L^{\text{sp}}}^{L^{\text{tot}}} \boldsymbol{\phi}\_{\text{sec}} \frac{3}{r\_{\text{sec}}^{-3}} \int\_{0}^{r\_{\text{sec}}} \boldsymbol{\chi}^{2} a\_{s\text{e}}^{\text{sec}} \overline{J}\_{s\text{e}}^{\text{sec}} \mathcal{F} \mathbf{d} \boldsymbol{\chi} \mathbf{x}, \end{aligned} \tag{4.75}$$

where *A* cell is the cross section area of the cell. Note that the result from Equation (4.73) was used. Additionally,

$$-\kappa\_{s, \text{eff}}^{\text{con,pos}} \nabla \overline{\varphi}\_s(L^{\text{sep}}) = 0 \tag{4.76}$$

accounts for the fact that electrons are not allowed to enter the separator domain. Finally, the electronic current density at the cathode-current-collector side is

$$\underbrace{-\kappa\_{s,\text{eff}}^{\text{en},\text{pos}}\nabla\overline{\varphi}\_{s}(L^{\text{tot}})}\_{\widetilde{\mathbb{F}}\_{s}(L^{\text{tot}})} = \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{3}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} y^{2} a\_{\text{se}}^{\text{sec}} \overline{j}\_{\text{se}}^{\text{sec}} \mathcal{F} \, \text{dyd} \,\text{x},\tag{4.77}$$

where i *s* (*L* tot) is set to be the applied electric current density iapp. Continuity of the electronic potential between secondary particle and the cell level is realized by a Dirichlet boundary condition

$$
\overline{\mathfrak{P}}\_s^{\rm sec}(r\_{\rm sec}) = \overline{\mathfrak{P}}\_s. \tag{4.78}
$$

### **Concentration in the solid phase**

Next, the average change in concentration of lithium in the solid phase is considered. Therefore, volume-averaged integration over the spherical primary particle level of Equation (4.59) brings

$$\begin{split} \frac{3}{r\_{\text{prim}}^{\text{-}3}} \int\_{0}^{r\_{\text{prim}}} z^{2} \frac{\partial c\_{s}^{\text{prim}}}{\partial t} \, \mathrm{d}z &= \frac{3}{r\_{\text{prim}}^{\text{-}3}} \int\_{0}^{r\_{\text{prim}}} \frac{\partial}{\partial z} \left( z^{2} D\_{s} \frac{\partial c\_{s}^{\text{prim}}}{\partial z} \right) \mathrm{d}z \\ &= \frac{3}{r\_{\text{prim}}^{\text{-}3}} r\_{\text{prim}}^{\text{-}2} D\_{s} \frac{\partial c\_{s}^{\text{-} \text{prim}}(r\_{\text{prim}})}{\partial z} - \frac{3}{r\_{\text{prim}}^{\text{-}3}} 0 \underline{^{2}} D\_{s} \frac{\partial c\_{s}^{\text{-} \text{prim}}(\theta)}{\partial z}, \end{split} \tag{4.79}$$

which describes the volume-averaged temporal change of lithium concentration inside a primary particle. Note that

$$D\_s \frac{\partial c\_s^{\text{prim}}(r\_{\text{prim}})}{\partial z} = \overline{j}\_{sc}^{\text{sec}} \tag{4.80}$$

accounts for the electrochemical reaction boundary condition at the surfaces of the primary particles. Here, *j* sec *se* is the Butler-Volmer type reaction term from Equation (4.69).

### **Ionic current and potential in the electrolyte phase**

The volume-averaged source term on the right-hand side of Equation (4.67) is tackled next. Therefore, Equation (4.63) is volume-averaged over its spherical domain, which brings

3 *r*sec 3 Z *<sup>r</sup>*sec 0 ∂ ∂ *y* −*y* 2 κ ion,sec *e*,eff ∂ϕ sec *e* ∂ *y* −*y* 2 κ ion,sec *D*,eff ∂ ln*c* sec *e* ∂ *y* d*y* = − 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se j* sec *se* F d*y* 3 *r*sec 3 −*r*sec 2 κ ion,sec *e*,eff ∂ϕ sec *e* (*r*sec) ∂ *y* −*r*sec 2 κ ion,sec *D*,eff ∂ ln*c* sec *e* (*r*sec) ∂ *y* − ✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✿<sup>0</sup> 3 *r*sec 3 −0 2 κ ion,sec *e*,eff ∂ϕ sec *e* (0) ∂ *y* −0 2 κ ion,sec *D*,eff ∂ ln*c* sec *e* (0) ∂ *y* = − 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se j* sec *se* F d*y* . (4.81)

Identifying that

$$\begin{split} \underbrace{-\kappa\_{e,\text{eff}}^{\text{ion,sec}} \frac{\partial \overline{\partial}\_{e}^{\text{sec}} (r\_{\text{sec}})}{\partial \mathbf{y}} - \kappa\_{D,\text{eff}}^{\text{ion,sec}} \frac{\partial \ln \overline{c}\_{e}^{\text{sec}} (r\_{\text{sec}})}{\partial \mathbf{y}}}\_{\Gamma\_{e,\text{surf}}^{\text{sec}}} &= \\ & - \frac{1}{r\_{\text{sec}}^{\text{c}}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} a\_{\text{se}}^{\text{sec}} \overline{J}\_{\text{se}}^{\text{sec}} \mathcal{F} \mathbf{dy} \end{split} \tag{4.82}$$

is the ionic current density at the surface of the secondary particle, the volumeaveraged ionic current density is

$$\tilde{\mathbf{i}}\_{e,\text{sec}} = \frac{3}{r\_{\text{sec}}} \overline{\mathbf{i}}\_{e,\text{surf}}^{\text{sec}} = -\frac{3}{r\_{\text{sec}}} \int\_0^{r\_{\text{sec}}} \mathbf{y}^2 d\_{\text{se}}^{\text{sec}} \overline{\mathbf{j}}\_{\text{se}}^{\text{sec}} \mathcal{F} \mathbf{d} \mathbf{y},\tag{4.83}$$

where i sec *<sup>e</sup>*,surf is the ionic current density at the surfaces of the secondary particles. The right-hand side of Equation (4.67) is rewritten as

$$\boldsymbol{\phi}\_{\rm sec} \mathbf{\tilde{i}}\_{s,\rm sec} = \underbrace{\boldsymbol{\phi}\_{\rm sec} \, \frac{\mathfrak{J}}{r\_{\rm sec}} \, \mathbf{\tilde{i}}\_{e,\rm surf}^{\rm sec}}\_{a\_{\rm sc}} = a\_{\rm sc} \mathbf{\tilde{i}}\_{e,\rm surf}^{\rm sec} \, . \tag{4.84}$$

The total ionic current, i.e. volume-averaged integration of Equation (4.67), is

$$\begin{split} A^{\text{cell}} \int\_{0}^{L^{\text{sep}}} \nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion,sep}} \nabla \overline{\boldsymbol{\varphi}}\_{e} - \kappa\_{D, \text{eff}}^{\text{ion,sep}} \nabla \ln \overline{c}\_{e} \right) d\mathbf{x} + \\ A^{\text{cell}} \int\_{L^{\text{sep}}}^{L^{\text{cut}}} \nabla \cdot \left( -\kappa\_{e, \text{eff}}^{\text{ion,pos}} \nabla \overline{\boldsymbol{\varphi}}\_{e} - \kappa\_{D, \text{eff}}^{\text{ion,pos}} \nabla \ln \overline{c}\_{e} \right) d\mathbf{x} = \\ A^{\text{cell}} \int\_{L^{\text{sep}}}^{L^{\text{cut}}} \phi\_{\text{sec}} \overline{\boldsymbol{i}}\_{e, \text{sec}} \mathbf{dx}. \end{split} \tag{4.85}$$

The left-hand side is

*A* cell <sup>Z</sup> *<sup>L</sup>* sep 0 ∇ · −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* −κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* d*x*+ *A* cell <sup>Z</sup> *<sup>L</sup>* tot *L* sep ∇ · −κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* −κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* d*x* = ✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ *A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* sep)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* sep) − *A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (0)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (0) + −*A* cell ✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✘✿<sup>0</sup> κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* tot)−κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* tot) − ✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭✭ *A* cell −κ ion,pos *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (*L* sep)−κ ion,pos *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (*L* sep) = −*A* cell −κ ion,sep *<sup>e</sup>*,eff ∇ϕ*<sup>e</sup>* (0)−κ ion,sep *<sup>D</sup>*,eff ∇ln*c<sup>e</sup>* (0) , (4.86)

where

$$-\kappa\_{e,\text{eff}}^{\text{ion,pos}}\nabla\overline{\varphi}\_e(L^{\text{tot}}) - \kappa\_{D,\text{eff}}^{\text{ion,pos}}\nabla\ln\overline{c}\_e(L^{\text{tot}}) = 0\tag{4.87}$$

accounts for the fact that no ionic current is allowed to enter or leave the cathodecurrent-collector interface. Also, continuity of ionic current at the separatorcathode interface is represented by

$$\begin{split} -\mathsf{x}\_{e,\mathrm{eff}}^{\mathrm{ion,sep}} \nabla \overline{\mathsf{\varPhi}}\_{e}(L^{\mathrm{sep}}) - \mathsf{x}\_{D,\mathrm{eff}}^{\mathrm{ion,sep}} \nabla \ln \overline{\mathsf{c}}\_{e}(L^{\mathrm{sep}}) &= \\ -\mathsf{x}\_{e,\mathrm{eff}}^{\mathrm{ion,pos}} \nabla \overline{\mathsf{\varPhi}}\_{e}(L^{\mathrm{sep}}) - \mathsf{x}\_{D,\mathrm{eff}}^{\mathrm{ion,pos}} \nabla \ln \overline{\mathsf{c}}\_{e}(L^{\mathrm{sep}}) \,. \end{split} \tag{4.88}$$

Using Equation (4.83), the right-hand side of Equation (4.85) becomes

$$A^{\rm cell} \int\_{L^{\rm sep}}^{L^{\rm tot}} \phi\_{\rm sec} \overline{\mathbf{i}}\_{\rm e,sec} \, \mathrm{d}x = -A^{\rm cell} \int\_{L^{\rm sep}}^{L^{\rm tot}} \phi\_{\rm sec} \frac{\mathfrak{J}}{r\_{\rm sec}^{-3}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^2 d\_{\rm se}^{\rm sec} \overline{J}\_{\rm se}^{\rm sec} \mathcal{F} \mathrm{d}y \, \mathrm{d}x. \tag{4.89}$$

Finally, the total ionic current density at the anode-separator interface is

$$\begin{split} -\underbrace{\left(-\kappa\_{e,\text{eff}}^{\text{ion,sep}}\nabla\overline{\rho}\_{e}(0) - \kappa\_{D,\text{eff}}^{\text{ion,sep}}\nabla\ln\overline{c}\_{e}(0)\right)}\_{\mathbb{T}\_{e}(0)} &= \\ -\int\_{L^{\text{sep}}}^{L^{\text{out}}}\phi\_{\text{sec}}\frac{3}{r\_{\text{sec}}^{-3}}\int\_{0}^{r\_{\text{sec}}}\mathbf{y}^{2}d\_{\text{se}}^{\text{sec}}\overline{f}\_{\text{se}}^{\text{sec}}\mathcal{F}\mathbf{d}\mathbf{y}d\mathbf{x},\end{split} \tag{4.90}$$

where, comparing to Equation (4.77), i *e* (0) is the applied current density iapp. Note, the continuity of the ionic potential between secondary particle and the cell level is realized by a Dirichlet boundary condition

$$
\overline{\mathfrak{P}}\_e^{\mathrm{sec}}(r\_{\mathrm{sec}}) = \overline{\mathfrak{P}}\_e. \tag{4.91}
$$

### **Concentration in the electrolyte phase**

As a final step, the volume-averaged source term on the right-hand side of Equation (4.68) is dealt with by volume-averaged integration of Equation (4.65), which results in

$$\begin{split} \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 \phi\_e^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sc}}{\partial t} \, \mathrm{d}y &= \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \frac{\partial}{\partial \mathbf{y}} \left( \mathbf{y}^2 D\_{e,\rm eff}^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sc}}{\partial \mathbf{y}} \right) \\ &- \mathbf{y}^2 a\_{\rm sc}^{\rm sc} (1 - t\_+^0) \overline{f}\_{\rm se}^{\rm sc} \, \mathrm{d}y \\ &= \frac{3}{r\_{\rm sec}} r\_{\rm sec} \, ^2 D\_{e,\rm eff}^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sc} (r\_{\rm sc})}{\partial \mathbf{y}} \\ &- \frac{3}{r\_{\rm sec}} \overline{\mathbf{J}}^2 D\_{e,\rm eff}^{\rm sec} \frac{\partial \overline{c}\_e^{\rm sc} (\mathbf{0})}{\partial \mathbf{y}} \\ &- \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sc}} \mathbf{y}^2 a\_{\rm sc}^{\rm sc} (1 - t\_+^0) \overline{f}\_{\rm se}^{\rm sc} \, \mathrm{d}y. \end{split} \tag{4.92}$$

The flux density at the secondary particle surface, *j* sec *<sup>e</sup>*,surf, is identified as

$$\begin{split} \underbrace{-D\_{e,\text{eff}}^{\text{sec}} \frac{\partial \overline{c}\_{e}^{\text{sec}}(r\_{\text{sec}})}{\partial \mathbf{y}}}\_{\overline{f\_{e}}, \text{surf}} &= -\frac{1}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} \phi\_{e}^{\text{sec}} \frac{\partial \overline{c}\_{e}^{\text{sec}}}{\partial t} \, \text{dy} \\ & \qquad -\frac{1}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} a\_{s\text{e}}^{\text{sec}} (1 - t\_{+}^{0}) \overline{f\_{s}}^{\text{sec}} \, \text{dy} \end{split} (4.93)$$

thus

$$\begin{split} \overline{J}\_{e,\text{sec}} = \frac{3}{r\_{\text{sec}}} \overline{J}\_{e,\text{surf}}^{\text{sec}} = -\frac{3}{r\_{\text{sec}}} \int\_0^{r\_{\text{sec}}} \chi^2 \Phi\_e^{\text{sec}} \frac{\partial \overline{c}\_e^{\text{sec}}}{\partial t} \,\text{dy} \\ & - \frac{3}{r\_{\text{sec}}} \int\_0^{r\_{\text{sec}}} \chi^2 a\_{\text{se}}^{\text{sec}} (1 - t\_+^0) \overline{j}\_{\text{se}}^{\text{sec}} \,\text{dy}, \end{split} \tag{4.94}$$

where *j* sec *<sup>e</sup>*,surf is the mass flux density at the secondary particle surfaces, which renders the corresponding term in Equation (4.68) as

$$\phi\_{\rm sec} \overline{j}\_{e,\rm sec} = \underbrace{\phi\_{\rm sec} \frac{3}{r\_{\rm sec}}}\_{a\_{\rm sc}} \overline{j}\_{e,\rm surf}^{\rm sec} = a\_{\rm se} \overline{j}\_{e,\rm surf}^{\rm sec} \,. \tag{4.95}$$

Applying volume-averaged integration on Equation (4.68) and using the result from Equation (4.94) leads to

1 *L* tot <sup>Z</sup> *<sup>L</sup>* sep 0 φ sep *e* ∂ *c<sup>e</sup>* ∂*t* d*x*+ Z *<sup>L</sup>* tot *L* sep φ pos *e* ∂ *c<sup>e</sup>* ∂*t* d*x* ! = 1 *L* tot <sup>Z</sup> *<sup>L</sup>* sep 0 ∇ · *D* sep *<sup>e</sup>*,eff∇*c<sup>e</sup>* d*x*+ Z *<sup>L</sup>* tot *L* sep ∇ · *D* pos *<sup>e</sup>*,eff∇*c<sup>e</sup>* +φsec *je*,sec d*x* ! = 1 *L* tot <sup>Z</sup> *<sup>L</sup>* sep 0 ∇ · *D* sep *<sup>e</sup>*,eff∇*c<sup>e</sup>* d*x*+ Z *<sup>L</sup>* tot *L* sep ∇ · *D* pos *<sup>e</sup>*,eff∇*c<sup>e</sup>* d*x* ! + 1 *L* tot <sup>Z</sup> *<sup>L</sup>* tot *L* sep <sup>φ</sup>sec − 3 *r*sec Z *<sup>r</sup>*sec 0 *y* 2 φ sec *e* ∂ *c* sec *e* ∂*t* d*y* − 3 *r*sec Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se* (1−*t* 0 <sup>+</sup>)*j* sec *se* d*y* d*x* ! = 1 *L* tot <sup>Z</sup> *<sup>L</sup>* sep 0 ∇ · *D* sep *<sup>e</sup>*,eff∇*c<sup>e</sup>* d*x*+ Z *<sup>L</sup>* tot *L* sep ∇ · *D* pos *<sup>e</sup>*,eff∇*c<sup>e</sup>* d*x* ! − 1 *L* tot <sup>Z</sup> *<sup>L</sup>* tot *L* sep φsec 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 φ sec *e* ∂ *c* sec *e* ∂*t* d*y*d*x* − 1 *L* tot <sup>Z</sup> *<sup>L</sup>* tot *L* sep φsec 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se* (1−*t* 0 <sup>+</sup>)*j* sec *se* d*y* (4.96)

The divergence terms read

$$\begin{split} & \int\_{0}^{L^{\text{sep}}} \nabla \cdot \left( D\_{\varepsilon,\text{eff}}^{\text{sep}} \nabla \overline{c}\_{\varepsilon} \right) \text{d\mathfrak{x}} + \int\_{L^{\text{sep}}}^{L^{\text{ion}}} \nabla \cdot \left( D\_{\varepsilon,\text{eff}}^{\text{pos}} \nabla \overline{c}\_{\varepsilon} \right) \text{d\mathfrak{x}} \\ &= D\_{\varepsilon,\text{eff}}^{\text{sep}} \nabla \overline{c}\_{\varepsilon} \langle \mathsf{L}^{\text{sep}} \overline{\rangle} - D\_{\varepsilon,\text{eff}}^{\text{sep}} \nabla \overline{c}\_{\varepsilon}(0) + D\_{\varepsilon,\text{eff}}^{\text{pos}} \Sigma \overline{c}\_{\varepsilon} \langle \mathsf{L}^{\text{tor}} \rangle^{\bullet} - D\_{\varepsilon,\text{eff}}^{\text{sep}} \Sigma \overline{c}\_{\varepsilon} \langle \mathsf{L}^{\text{sep}} \rangle \\ &= -D\_{\varepsilon,\text{eff}}^{\text{sep}} \nabla \overline{c}\_{\varepsilon}(0), \end{split} \tag{4.97}$$

where

$$D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(L^{\text{sep}}) = D\_{e, \text{eff}}^{\text{pos}} \nabla \overline{c}\_e(L^{\text{sep}}) \quad \text{and} \quad D\_{e, \text{eff}}^{\text{pos}} \nabla \overline{c}\_e(L^{\text{tot}}) = 0 \tag{4.98}$$

reflects continuity of mass flux between separator and cathode region and implies that no mass flux is allowed at the cathode-current-collector side.

Using Equation (4.97) in Equation (4.96) yields

$$\begin{split} &\frac{1}{L^{\text{tot}}} \left( \int\_{0}^{L^{\text{sep}}} \phi\_{\epsilon}^{\text{sep}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \mathrm{d}x + \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\epsilon}^{\text{pos}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \mathrm{d}x \right) = -\frac{1}{L^{\text{tot}}} \left( D\_{\epsilon, \text{eff}}^{\text{sep}} \nabla \overline{c}\_{\epsilon}(0) \right) \\ & \qquad \quad - \frac{1}{L^{\text{tot}}} \left( \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{\overline{3}}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \chi^{2} \phi\_{\epsilon}^{\text{sec}} \frac{\partial \overline{c}\_{\epsilon}^{\text{sec}}}{\partial t} \, \mathrm{d}y \, \mathrm{d}x \right) \\ & \qquad \quad - \frac{1}{L^{\text{tot}}} \left( \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{\overline{3}}{r\_{\text{sec}}} \int\_{0}^{r\_{\text{sec}}} \chi^{2} a\_{se}^{\text{sec}} (1 - t\_{+}^{0}) \overline{f}\_{se}^{\text{sec}} \, \mathrm{d}y \, \mathrm{d}x \right) . \end{split} \tag{4.99}$$

Finally, the volume-averaged temporal change of concentration in the electrolyte phase is

$$\begin{split} \frac{1}{L^{\text{tot}}} \left( \int\_{0}^{L^{\text{sep}}} \phi\_{\epsilon}^{\text{sep}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \text{d}\mathbf{v} + \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\epsilon}^{\text{pos}} \frac{\partial \overline{c}\_{\epsilon}}{\partial t} \, \text{d}\mathbf{x} \right) \\ &+ \frac{1}{L^{\text{tot}}} \left( \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{3}{r\_{\text{sec}}^{-3}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} \phi\_{\epsilon}^{\text{sec}} \frac{\partial \overline{c}\_{\epsilon}^{\text{sec}}}{\partial t} \, \text{d}\mathbf{y} \, \text{d}\mathbf{x} \right) = \\ &- \frac{1}{L^{\text{tot}}} \left( D\_{\epsilon, \text{eff}}^{\text{sep}} \nabla \overline{c}\_{\epsilon}(0) \right) - \frac{1}{L^{\text{tot}}} \left( \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{3}{r\_{\text{sec}}^{-3}} \int\_{0}^{r\_{\text{sec}}} \mathbf{y}^{2} a\_{s\text{e}}^{\text{sec}} (1 - t\_{+}^{0}) \overline{J}\_{s\epsilon}^{\text{sec}} \, \text{d}\mathbf{y} \, \text{d}\mathbf{x} \right). \end{split} \tag{4.100}$$

Equation (4.100) shows that the volume-averaged temporal change of electrolyte concentration is due to diffusion at the anode-separator interface on cell level and electrochemical reactions of cations at the electrolyte-primary-particle interface on secondary particle level. As was assumed earlier, the applied electroneutrality condition enforces the same amount of anions and cations. On cell and secondary particle level, this is expressed as a constant volume-averaged concentration in the electrolyte, i.e. a vanishing temporal change. Therefore, the left-hand side of Equation (4.100) is set to zero yielding

$$-D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(0) = \int\_{L^{\text{sep}}}^{L^{\text{tot}}} \phi\_{\text{sec}} \frac{3}{r\_{\text{sec}}} \int\_0^{r\_{\text{sec}}} \mathbf{y}^2 a\_{\text{se}}^{\text{sec}} (1 - t\_+^0) \overline{f}\_{\text{se}}^{\text{sec}} \, \text{dyd} \, \text{x} \,. \tag{4.101}$$

By comparing Equations (4.90) and (4.101), the latter of which is rewritten as

Z *<sup>L</sup>* tot *L* sep φsec 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se* (1−*t* 0 <sup>+</sup>)*j* sec *se* d*y*d*x* = (1−*t* 0 +) F <sup>Z</sup> *<sup>L</sup>* tot *L* sep φsec 3 *r*sec 3 Z *<sup>r</sup>*sec 0 *y* 2 *a* sec *se j* sec *se* F d*y*d*x* ! | {z } i *e* (0)=iapp = i app(1−*t* 0 +) F , (4.102)

where *t* 0 <sup>+</sup> is assumed to be constant. Using the boundary condition from above, i.e. setting the ionic current density at the anode-separator interface equal to the applied current density, yields

$$-D\_{e, \text{eff}}^{\text{sep}} \nabla \overline{c}\_e(0) = \frac{\mathbf{i}\_{\text{app}}}{\mathcal{F}} (1 - t\_+^0) \,. \tag{4.103}$$

Interestingly, Equation (4.103) is the same as for the half-cell model [40, 42, 132]. In order to satisfy continuity of concentration, the concentration between secondary particle surface and the cell level is prescribed as Dirichlet boundary condition

$$
\overline{c}\_e^{\text{sec}}(r\_{\text{sec}}) = \overline{c}\_e \,. \tag{4.104}
$$

### **Cell quantities**

Analogous to the half-cell model from Section 4.1.2, the evaluation of the cell performance is done by computing the appropriate cell quantities. Note that the cell voltage and state of charge or depth of discharge is computed similar to the half-cell model from Section 4.1.2 and are therefore omitted here.

### **Specific capacity**

The specific capacity is calculated by first evaluating the total capacity as

$$\begin{split} \mathcal{Q}\_{\text{tot}} &= V\_s^{\text{pos}} (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F} \\ &= A^{\text{cell}} L^{\text{pos}} \phi\_s \phi\_s^{\text{sec}} (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}, \end{split} \tag{4.105}$$

where *V* pos *<sup>s</sup>* is the volume of the solid phase inside the positive electrode, and *cs*,init and *cs*,max is the initial and maximum allowable concentration of the active material, respectively. The former of which, can be calculated using the positive electrode's cross section area *A* cell, length *L* pos, volume fraction of secondary particles φ*<sup>s</sup>* and volume fraction of the secondary particles' solid phase φ sec *s* .

Finally, the specific capacity is

$$\begin{split} \mathcal{Q}\_{\text{spec}} &= \frac{\mathcal{Q}\_{\text{tot}}}{m\_{\text{AM}}} = \frac{A^{\text{cell}} L^{\text{pos}} \phi\_s \phi\_s^{\text{sec}} (c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}}{A^{\text{cell}} L^{\text{pos}} \phi\_s \phi\_s^{\text{sec}} \,\mathcal{Q}\_{\text{AM}}} \\ &= \frac{(c\_{s,\text{max}} - c\_{s,\text{init}}) \mathcal{F}}{\mathcal{Q}\_{\text{AM}}} \,, \end{split} \tag{4.106}$$

where %AM is the density of the active material.

### **Calculating the** *nC* **current**

In case of the presented hierarchically structured half-cell model, the volumeaveraged temporal change of lithium concentration in the positive electrode is calculated as

$$\frac{\partial \overline{c}\_{s, \text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \int\_{L^{\text{seq}}}^{L^{\text{tot}}} \frac{3}{r\_{\text{sec}}^{-3}} \int\_{0}^{r\_{\text{sec}}} y^2 \frac{3}{r\_{\text{prim}}^{-3}} \int\_{0}^{r\_{\text{prim}}} z^2 \frac{\partial c\_s^{\text{prim}}}{\partial t} \, \text{d}z \, \text{dy} \, \text{d}x. \tag{4.107}$$

From Equations (4.79) and (4.80), it is known that

$$\frac{3}{r\_{\rm prim}} \int\_0^{r\_{\rm prim}} z^2 \frac{\partial c\_s^{\rm prim}}{\partial t} \, \mathrm{d}z = \frac{3}{r\_{\rm prim}} \overline{j}\_{se}^{\rm sec} \,. \tag{4.108}$$

Importing Equation (4.108) into Equation (4.107) yields

$$\frac{\partial \overline{c}\_{s,\text{avg}}}{\partial t} = \frac{1}{L^{\text{pos}}} \int\_{L^{\text{exp}}}^{L^{\text{tot}}} \frac{\mathfrak{J}}{r\_{\text{sec}}^{-\mathfrak{J}}} \int\_{0}^{r\_{\text{sec}}} y^2 \frac{\mathfrak{J}}{r\_{\text{prim}}} \overline{j}\_{\text{se}}^{\text{sec}} \, \text{dy} \, \text{dx} \, . \tag{4.109}$$

Moreover, from Equation (4.77), it was found that

$$\int\_{L^{\rm sep}}^{L^{\rm tot}} \Phi\_{\rm sec} \frac{\mathfrak{J}}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 a\_{\rm se}^{\rm sec} \overline{J}\_{\rm se}^{\rm sec} \mathcal{F} \mathbf{d} \mathbf{y} \, \mathbf{dx} = \mathbf{i}\_{\rm app} \,, \tag{4.110}$$

if galvanostatic boundary conditions apply and i *s* (*L* tot) is set to a constant applied electric current density iapp. Under the assumption that φsec, *r* prim and *a* sec *se* are constant properties, Equations (4.109) and (4.110) are rewritten as

$$\underbrace{\frac{1}{L^{\rm pos}} \int\_{L^{\rm sep}}^{L^{\rm tot}} \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 \frac{3}{r\_{\rm prim}} \overline{f}\_{\rm se}^{\rm sec} \, \mathbf{dy} \, \mathbf{dx}} = \frac{1}{L^{\rm pos}} \frac{3}{r\_{\rm prim}} \int\_{L^{\rm sep}}^{L^{\rm post}} \frac{3}{r\_{\rm sec}} \int\_0^{r\_{\rm sec}} \mathbf{y}^2 \overline{f}\_{\rm se}^{\rm sec} \, \mathbf{dy} \, \mathbf{dx} \tag{4.111}$$

and

$$\underbrace{\int\_{L^{\rm{spc}}}^{L^{\rm{tot}}} \Phi\_{\rm{sec}} \frac{3}{r\_{\rm{sec}}^{-3}} \int\_{0}^{r\_{\rm{spc}}} \text{y}^{2} a\_{\rm{se}}^{\rm{sec}} \overline{j}\_{\rm{se}}^{\rm{sec}} \mathcal{F} \text{dy} \,\mathrm{d}x} = \Phi\_{\rm{ssc}} d\_{\rm{se}}^{\rm{sec}} \mathcal{F} \int\_{L^{\rm{spc}}}^{L^{\rm{s\alpha}}} \frac{3}{r\_{\rm{sec}}^{-3}} \int\_{0}^{r\_{\rm{s\alpha}}} \text{y}^{2} \overline{j}\_{\rm{se}}^{\rm{sec}} \text{dy} \,\mathrm{d}x. \tag{4.112}$$

Comparing Equations (4.111) and (4.112) yields

$$\begin{split} \frac{\partial \overline{c}\_{s, \text{avg}}}{\partial t} &= \mathbf{i}\_{\text{app}} \frac{1}{L^{\text{pos}}} \frac{3}{r\_{\text{prim}}} \frac{1}{a\_{se}^{\text{sec}}} \frac{1}{\phi\_{\text{sec}}} \frac{1}{\mathcal{F}} \\ &= \frac{\mathbf{i}\_{\text{app}}}{\mathcal{F}} \frac{3}{r\_{\text{sec}} a\_{se}} \frac{3}{r\_{\text{prim}} a\_{se}^{\text{sec}}}, \end{split} \tag{4.113}$$

where by φsec = *ase r*sec 3 it was assumed that the secondary particles are equal-sized and detached [136]. *cs*,avg(*t*) at time *t* is calculated analogous to Equation (4.55), which is inserted into Equation (4.45) such that the depth of discharge is calculated as *cs*,avg(*t*) at time *t* is calculated by multiplying Equation (4.113) by *t* and adding the initial concentration *cs*,init yielding

$$\begin{split} \overline{c}\_{s,\text{avg}}(t) &= \frac{\overline{\partial}\overline{c}\_{s,\text{avg}}}{\overline{\partial}t} \cdot t + c\_{s,\text{init}} \\ &= \frac{\mathrm{i}\_{\mathrm{app}}}{\mathcal{F}} \frac{\mathfrak{Z}}{r\_{\mathrm{sec}}a\_{\mathrm{se}}} \frac{\mathfrak{Z}}{r\_{\mathrm{prim}}a\_{\mathrm{se}}^{\mathrm{sec}}} \cdot t + c\_{s,\text{init}} \cdot \mathrm{d} \end{split} \tag{4.114}$$

Inserting Equation (4.114) in Equation (4.45), the depth of discharge at a time *t* is

$$\text{DOD}(t) = \frac{\frac{\text{i}\_{\text{app}}}{\mathcal{F}} \frac{3}{r\_{\text{scc}} a\_{\text{sc}}} \frac{3}{r\_{\text{prim}} a\_{\text{sc}}^{\text{sec}}} \cdot t}{c\_{s, \text{max}} - c\_{s, \text{init}}} \,. \tag{4.115}$$

The 1*C*-current can be computed by rearranging Equation (4.115) as well as setting *t* = 3600 s and DOD(*t*) = 1 representing a fully discharged cell within one hour as

$$\mathbf{i}\_{\rm app}^{\rm lC} = (c\_{s,\rm max} - c\_{s,\rm init}) \frac{L^{\rm pos} \, r\_{\rm sec} \, a\_{\rm sc} \, r\_{\rm prim} \, a\_{\rm sc}^{\rm sec} \, \mathcal{F}}{9 \cdot 3600} \,. \tag{4.116}$$

In general, the *nC*-current is

$$\mathbf{i}\_{\rm app}^{\rm nC} = n(c\_{s,\rm max} - c\_{s,\rm init}) \frac{L^{\rm pos} r\_{\rm sec} a\_{\rm se} r\_{\rm prim} a\_{\rm se}^{\rm sec} \mathcal{F}}{9 \cdot 3600}. \tag{4.117}$$

Note that the above derivations do only apply for cases where the specific surface area can be calculated using *ase* = 3φ*<sup>s</sup> r*sec and *a* sec *se* = 3φ sec *s r* prim , which implies the assumption that the active material particles, either secondary or primary particles, are of spherical shape and detached from each other [136]. Morphological variations from these formulas necessitate special treatment of the above presented derivations, which shall not be considered in the present work.

# **4.3 Validation of the hierarchically structured half-cell model**

The following section aims at validating the previously presented hierarchically structured half-cell model [137]. To this end, geometry, structure and material properties from [30] were imported into the model and the simulation results were compared to measurements reported in [30].

### **4.3.1 Geometry, structure and transport properties**

In [30], both classical and hierarchically structured electrodes were prepared and electrochemically characterized. Moreover, the morphology of hierarchically structured electrodes was varied in terms of primary particle sizes, inner porosity and secondary particle size. Additionally, statistical image analysis based on synchrotron tomography was used to investigate the structural properties of the electrodes.

Concerning classical electrodes, dense pristine powder (p-NMC) was used, where LiNi1/3Mn1/3Co1/3O<sup>2</sup> (NMC) was the active material. The hierarchically structured electrodes were prepared using nanostructured powder (n-NMC), which was obtained by grinding, spray drying and calcinating the p-NMC powder. See [30] for a more detailed description of the process. In both cases, slurries were produced by adding conductive filler material to the powders. The conductive filler material comprised of polyvinylidene difluoride (PVDF) binder, carbon black and graphite. Finally, the electrodes were produces by casting the slurries onto an aluminium foil. In Tables 4.1 and 4.2, the geometry and structure properties of the resulting electrodes are summarized. Note that the specific surface area is computed as *a* = 3φ *r* [136], where—depending on model level—the respective value is calculated by

$$a\_{se} = \frac{3\,\phi\_s^{\rm pos}}{r\_{\rm sec}} \quad \text{or} \quad a\_{se}^{\rm sec} = \frac{3\,\phi\_s^{\rm sec}}{r\_{\rm prim}},\tag{4.118}$$

respectively.

Figure 4.4: Mathematical model of the hierarchically structured half-cell setup.

4.3 Validation of the hierarchically structured half-cell model

Figure 4.5: Mathematical model of the hierarchically structured half-cell setup. (continued)


Table 4.1: Geometry properties of the classical and hierarchically structured half-cells taken from [30].

In the following, the classical electrode setup named p-NMC is compared to hierarchically structured electrodes denoted by n-NMC, which are based on the same material. Moreover, the secondary particle's radii are comparable among all electrodes. Additionally, the nanostructured secondary particles were treated with different calcination temperatures of 850°C and 900°C, which led to larger primary particle sizes and lower porosities in case of the higher temperature. In order to distinguish between those two, the n-NMC electrodes are denoted by n-NMC-F850 and n-NMC-F900.

The effective transport properties are calculated using the bulk material property multiplied by an effective transport parameter <sup>ˆ</sup>*k*eff, which incorporates the morphology of the transport paths of the corresponding conducting phase. On the one hand, it can be calculated by employing the resistor network method for both the solid and pore networks on virtual particle assemblies, as presented in Sections 3.1 and 3.1. On the other hand, the statistically derived socalled M-factor can be used equivalently which accounts for volume fraction, tortuosity and constrictivity of the transport phase [92, 93]. One prominent way of calculating effective transport parameters, however, is the Bruggeman correlation [71, 138], which is a function of the volume fraction of the transport phase and the Bruggeman transport coefficient brugg in the form of <sup>ˆ</sup>*k*eff <sup>=</sup> <sup>φ</sup> brugg . As for the cell level in this investigation, this type of relation shall be used to compute the effective transport parameters of the cathode as follows.


Table 4.2: Structure properties of the classical and hierarchically structured half-cells taken from [30].

*a*) Assumed. *b*) Calculated using Equation (4.118).

In [92], the RN was extensively used to generate a large database of sphere assembly scenarios with randomly distributed and slightly overlapping particles following a polydisperse size-distribution. By means of statistical methods, prediction formulas for the transport parameters of the considered packings were developed in terms of the mean contact angle and the standard deviation of the particle radii. Notably, it was found that a Bruggeman type correlation can be employed to compute effective transport parameters of the pore phase. In contrast to the transport through the solid phase, it seemed that neither the considerably large polydispersity nor the degree of the investigated densification had an influence on the transport through the pore phase. There, a Bruggeman coefficient of 1.342 was found to achieve the best approximation [92]. Considering that the underlying degree of polydispersity and densification corresponds to the cathode structures in [30], the transport coefficient is therefore applied in this investigation.

In addition to that, as a simple model assumption, the electronic and ionic transport carried by the filler phase and electrolyte phase, respectively, is modeled by splitting the above-mentioned Bruggeman-based effective transport parameter of the total pore phase according to the volume fractions. This way, both phases can be regarded as smeared homogeneously over the pore phase. Recall from before that the filler phase comprises of binder, carbon black and graphite. In this work, the filler phase is equipped with an effective electronic conductivity resulting from all three components.

In this investigation, using the volume fractions φ pos *f* of the filler and φ pos *<sup>e</sup>* of the electrolyte phase from Table 4.2, the volume fraction of the pore phase is obtained as

$$
\phi\_{\text{pore}}^{\text{pos}} = \phi\_e^{\text{pos}} + \phi\_f^{\text{pos}}, \tag{4.119}
$$

which is used to calculate the effective transport parameter of the pore phase as

$$
\hat{k}\_{\rm eff,pore}^{\rm pos} = \phi\_{\rm pore}^{\rm pos\ 1.342}\,. \tag{4.120}
$$

In the framework of the cell models presented in this work, electronic and ionic transport is via the filler and electrolyte phase, respectively. Using the additive split of the pore phase from Equation (4.119), the effective electronic and ionic transport properties are assigned according to the volumetric share of the filler and electrolyte phase as

$$
\hat{k}\_{\rm eff}^{\rm ion,pos} = \phi\_{\rm pore}^{\rm pos\ 1.342} \cdot \left(\frac{\phi\_f^{\rm pos}}{\phi\_{\rm pore}^{\rm pos}}\right) \quad \text{and} \quad \hat{k}\_{\rm eff}^{\rm ion,pos} = \phi\_{\rm pore}^{\rm pos\ 1.342} \cdot \left(\frac{\phi\_e^{\rm pos}}{\phi\_{\rm pore}^{\rm pos}}\right), \tag{4.121}
$$

respectively. In case of the separator on the cell in this investigation, its porosity was assumed as 0.5 , which is in the range of the typical characteristics of commercial separators, and its Bruggerman transport coefficient was taken as 3.0, which is within the range reported in literature [139].

Finally, effective transport parameters of the primary particle and pore networks on secondary particle level for the hierarchically structured half-cell model were used in form of the M-factors reported in [30], which are based on the real morphology. The effective transport parameters are summarized in Table 4.3.


Table 4.3: Transport properties of classical and hierarchically structured electrodes taken from [30].

*a*) Calculated using Equation (4.121).

## **4.3.2 Electrochemistry and material parameters**

From [30], the equilibrium potential curves, i.e. open circuit voltage curves (OCV), were taken as the ones measured at the lowest C rate of *C*/20, see Figure 4.6.

Figure 4.6: OCV curves taken from [30].

The half-cell models, as presented here, need information on the minimum and maximum allowable lithium concentration inside the active cathode material. Those quantities are computed here as follows. During the discharging process, lithium is intercalated into the active material starting from an initial depth of discharge DODinit. Arriving at the maximum capacity, the final maximum depth of discharge DODmax is reached. The initial and maximum depth of discharge is computed by assuming that—due to the stoichiometry of LiNi1/3Mn1/3Co1/3O2—the theoretical maximum capacity of NMC being around *Q* th.max spec = 278mA h g−<sup>1</sup> , see Appendix B.1. In [30], it was found that the reversibly accessible capacity in case of p-NCM is 158mA h g−<sup>1</sup> and in case of n-NCM-F900 is *Q* tot spec = 161mA h g−<sup>1</sup> . In this investigation, the same capacity of 161mA h g−<sup>1</sup> was used for n-NCM-F850. As a model assumption, here the cell is discharged until the maximum theoretical capacity of NMC is reached, which means DODmax = 1. This is justified by discussing Figure 4.6. It can be seen that at a cut-off voltage of 3V the gradients of the discharge curves tend towards infinity. This implies that the material's maximum capacity is almost reached. Additionally, due to the lithium metal as anode, it can be assumed that the supply of lithium ions is enough for the cathode material to fully lithiate. It should be noted that, during the initial charging and discharging processes, the capacity of the NMC material might not be fully exploited because of the formation of a lithium-rich phase at the surfaces of the active material particles which hinders lithium diffusion [140]. However, there it was also shown that by an extended voltage relaxation step at deep discharge, this phase disappears and the capacity efficiency can be fully recovered. In turn, theoretically, this means that the provided voltage window is enough to fully lithiate the given material, i.e. reaching the maximum theoretical capacity.

Subsequently, the onset of the simulation, i.e. initial depth of discharge DODinit , is calculated relative to DODmax. By using the reversibly accessible capacity from before, the initial depth of discharge is computed as

$$\text{DOD}^{\text{init}} = \frac{\mathcal{Q}\_{\text{spec}}^{\text{th.max}} - \mathcal{Q}\_{\text{spec}}^{\text{tot}}}{\mathcal{Q}\_{\text{spec}}^{\text{th.max}}},\tag{4.122}$$

which is 0.42 and 0.43 in cases of n-NCM-F850/n-NMC-F900 and p-NCM, respectively. Furthermore, the maximum and minimum concentration in the solid phase is computed by

$$\begin{split} c\_{s, \text{max}} &= \text{DOD}^{\text{max}} \frac{\mathcal{Q}\_{\text{spec}}^{\text{th.max}} \, \mathcal{Q}\_{\text{NMC}}}{\mathcal{F}} \text{ and } \\ c\_{s, \text{init}} &= \text{DOD}^{\text{init}} \frac{\mathcal{Q}\_{\text{spec}}^{\text{th.max}} \, \mathcal{Q}\_{\text{NMC}}}{\mathcal{F}} \,, \end{split} \tag{4.123}$$

respectively, where the NMC density of %NMC= 4770 kgm−<sup>3</sup> [141] was used, which can be calculated according to Appendix B.2.

In literature, values of NMC lithium diffusion coefficients can vary from 1 · 10−<sup>15</sup> m<sup>2</sup> s −1 to 1 · 10−<sup>13</sup> m<sup>2</sup> s −1 [142–146] and the electronic conductivity can vary from 2.2 · 10−<sup>4</sup> Sm−<sup>1</sup> to 5.2 · 10−<sup>6</sup> Sm−<sup>1</sup> [145, 147]. Additionally, in [142, 148], electronic conductivity was measured as a function of temperature and lithiation state of NMC. Here, the lithiation state is depicted by γ = *<sup>c</sup>s*/*cs*,max in LiγNi1/3Mn1/3Co1/3O2, where—theoretically—γ = 1 refers to a fully lithiated and γ = 0 to a fully delithiated state.

In Figure 4.7, the electronic conductivity of NMC is plotted over the lithiation state at 30°C, where the experimental values from [142] were approximated by the fit-function

$$\lg(\text{kg}^{\text{en}}\_{\text{NMC}}) = -2.988 \,\text{y}^{10.98} - 2.629 \,\text{y} + 0.624.\tag{4.124}$$

It can be observed that the value of conductivity ranges from ≈ 1 to ≈ 1 · 10−<sup>5</sup> Sm−<sup>1</sup> for lithiation states of γ = 0.25 to γ = 1.00, respectively. This result by [142] was acknowledged by [143]. There, additionally, the diffusion coefficient and reaction rate for NMC was measured. Therefore, in the following, material parameters were taken from [143] and [142] because a complete and, thus, more likely consistent set of diffusion coefficient, electronic conductivity and reaction rate constant for NMC can be expected.

Figure 4.7: Electronic conductivity of NMC over lithiation. Experimental values taken from [142] and fitted by Equation (4.124).

The kinetic reaction rate constant is often treated as a fitting parameter [149]. However, in case of different electrode structures, the fitting process could yield different reaction rate constants for each variation. In contrast to that, the approach in this work is to approximate the effective reaction rate constant *k*<sup>0</sup> purely based on the underlying material. By taking the exchange current density from experiments on NMC in [143], the exact same material based constant is used for every simulation in this work. This way, the comparability between different cathode and cell structures can be improved.

To compute the effective reaction rate constant, first, the Butler-Volmer type reaction equation from Equation (4.20) is employed. Based on this equation, the reaction rate constant can be calculated as

$$\begin{split}-\frac{i\_0}{\mathcal{F}} &= -k\_0 \overline{c}\_\epsilon^{1-\alpha} (c\_{s,\text{max}} - c\_{s,\text{surf}})^{1-\alpha} c\_{s,\text{surf}}^{\alpha} \\ k\_0 &= \frac{i\_0}{\mathcal{F}} \frac{1}{\overline{c}\_\epsilon^{1-\alpha} (c\_{s,\text{max}} - c\_{s,\text{surf}})^{1-\alpha} c\_{s,\text{surf}}^{\alpha}}. \end{split} \tag{4.125}$$

Secondly, the experimental value of 7.2Am−<sup>2</sup> is used for the exchange current density *i*0, which is taken from [143]. From Table 4.4, the mean value of maximum and minimum lithium concentration in the solid phase for p-NMC of 35402molm−<sup>3</sup> and the mean electrolyte concentration of 1000molm−<sup>3</sup> is used for *cs*,surf and *c<sup>e</sup>* , respectively. α is taken as 0.5. As a result, the effective reaction rate constant is recalculated as

$$k\_0 = \frac{7.2}{\mathcal{F}} \frac{1}{1000^{0.5} (49477 - 35402)^{0.5} 49477^{0.5}} \approx 1.0 \cdot 10^{-10} \,\text{s}.\tag{4.126}$$

Therefore, in all simulations in this work, *k*<sup>0</sup> is taken as 1.0 · 10−<sup>10</sup> mol−1m2.5/s.

The LP30 electrolyte (1 M LiPF<sup>6</sup> v/v: 1:1 DMC:EC) was considered, where the concentration dependent material parameters were taken from [131]. The initial concentration is 1000molm−<sup>3</sup> . In accordance with Section 4.1.1, *t* 0 <sup>+</sup> is set to a constant value of 0.23, which is the average value of *t* 0 <sup>+</sup> for the concentration *c<sup>e</sup>* between the assumed extreme values 0 and 3000molm−<sup>3</sup> .


Table 4.4: Material properties.

## **4.3.3 Results and Discussion**

The above described material, geometry and structural parameters were imported into the classical and the hierarchically structured half-cell model from Sections 4.1 and 4.2, respectively. In the following, discharging currents corresponding to different C rates were applied in each model, where the C rates ranged from rather low to rather large values. As in [30], C rates of *C* = 1/20, 1/10, 1/5,1/2,1,2,3,5,7,10 were used, which ultimately aims at revealing rate capability of the underlying cell composition.

The comparison of experiments and simulation results can be observed in Figure 4.8. Recall that the electrode structure with unmodified dense cathode particles is modeled by using the classical half-cell model and is denoted by p-NCM. Moreover, electrodes built-up by nanostructured secondary particles are modeled by means of the hierarchically structured half-cell model and are referred to as n-NCM-F850 and n-NCM-F900. The experimental values in Figure 4.8b were obtained by taking the mean value of the corresponding measurements in [30]. Notably, n-NCM-F850 and n-NCM-F900 show better rate capabilities as compared to p-NCM. Clearly, the same behavior can be observed in Figure 4.8a, where the simulation results also predict better performance for the nanostructured particle electrodes. Especially for rather high C rates of 7*C* and 10*C*, the hierarchically structured half-cell model results are very close to the measured values. While the simulations estimate around 70mA h g−<sup>1</sup> and 50mA h g−<sup>1</sup> retained specific capacity for these C rates, experiments measure a little bit higher values of around 75mA h g−<sup>1</sup> and 55mA h g−<sup>1</sup> . Note that in case of the classical half-cell model, the simulated value is a little bit higher in case of 10*C* than the measurements. The simulation predicts approximately 30mA h g−<sup>1</sup> , while experiments found approximately 25mA h g−<sup>1</sup> . Moreover, in case of both the classical and hierarchically structured half-cell model the retained specific capacity only decreases gradually for C rates between 1/20*C* and 3*C*. Apparently, the simulation curves show a distinct drop between 3*C* and 7*C*, whereas the experimental curves decrease steadily from low C rates to higher C rates.

Figure 4.8: Rate capabilities of p-NCM, n-NCM-F850 and n-NCM-F900. a) Simulation results. b) Experiments from [30].

Given that the results clearly show a qualitatively good agreement with experiments, the proposed hierarchically structured half-cell model can be regarded as validated in a sense that it allows for investigations to predict at least qualitatively the influence of geometry, structure and material on electrochemical performance. Moreover, as is shown in the following, the model offers detailed information on how those parameters affect the cell quantities. For instance, the drop behavior mentioned before can be understood by investigating the local distribution of DOD, denoted here as dod. It is calculated similar to Equation (4.45), however, the volume-averaged lithium concentration of the solid phase at time *t* is not evaluated for the whole cell, rather, the average concentration is computed at every point *x* along the cathode on cell level and at every point *y* along the secondary particle level.

On cell level, the local distribution is computed as

$$\text{dod}(t, \mathbf{x}) = \frac{\overline{c}\_s(t, \mathbf{x}) - c\_{s, \text{init}}}{c\_{s, \text{max}} - c\_{s, \text{init}}}. \tag{4.127}$$

The volume-averaged lithium concentration of the solid phase at time *t* at the position *x* is calculated as

$$\overline{c}\_{s}(t, \mathbf{x}) = \frac{\mathfrak{J}}{r\_{\rm sec} \, ^{\mathfrak{J}}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^{2} c\_{s}(t, \mathbf{x}) \, \mathbf{dy} \tag{4.128}$$

in case of the classical half-cell model, whereas this value is

$$\overline{\sigma}\_{s}(t, \mathbf{x}) = \frac{\mathfrak{Z}}{r\_{\rm sec}^{-3}} \int\_{0}^{r\_{\rm sec}} \mathbf{y}^{2} \frac{\mathfrak{Z}}{r\_{\rm prim}^{-3}} \int\_{0}^{r\_{\rm prim}} z^{2} c\_{s}^{\rm prim}(t, \mathbf{x}) \, \mathbf{dz} \, \mathbf{dy},\tag{4.129}$$

in case of the hierarchically structured half-cell model.

Figure 4.9 shows dod over the normalized cathode position for both the classical and hierarchically structured half-cell model. *t* is set as the end time of the simulation, which is determined by the cut-off voltage 3V. The cathode positions (*x* − *L* sep)/*L* pos = 0 and (*x* − *L* sep)/*L* pos = 1 refer to the separatorcathode interface and the cathode-current-collector interface, respectively, while (*x*−*L* sep)/*L* pos = 0.5 is the cathode middle position. Additionally, in Figure 4.9a - c, the state of charge distribution is shown for different C rates, namely 3*C*, 5*C* and 7*C*.

It can be observed that in case of 3*C*, the distribution is constant over the cathode thickness. This applies to p-NMC as well as n-NMC-F850 and n-NMC-F900. While dod is almost 1 in case of n-NMC-F850/n-NMC-F900, the depth of discharge is lower in case of p-NMC. At 5*C*, all curves show decreasing behavior towards the cathode-current-collector interface. In case of n-NMC-F850 and n-NMC-F900, a small gradient is visible in the first half of the cathode whereas in the second half the gradient is more pronounced. In case of p-NMC, the decreasing behavior can be observed directly at the separator-cathodeinterface. Overall, the dod values of the hierarchically structured half-cells are constantly larger than 0.60 whereas the classical half-cell model ranges around 0.45. Finally, at 7*C*, the depth of discharge distribution of n-NMC-F850 and n-NMC-F900 looks similar to the one of p-NMC at 5*C*. However, for the nanostructured electrodes, a gradient is visible right at the separator-cathode interface and the curve flatten towards the cathode-current collector interface. Moreover, the dod level of n-NMC-F850 and n-NMC-F900 is around the same value as well. The distribution in case of p-NMC is almost constant over the cathode region and is around 0.25.

By comparing the dod distribution of the different cells, an important aspect considering cathodes with nanostructured secondary particles can be extracted. For increasing C rates, the cathode active material is more efficiently utilized in case of hieararchically structured cathodes, which yields to superior rate performance as compared to cathodes with monolithic secondary particles. In addition to the dod distribution along the *x*-direction, the depth of discharge distribution can be investigated along the radius in *y*-direction inside the secondary particle, which is referred to as dodsec(*t*, *y*).

On secondary particle level, the depth of discharge distribution along the radial position is calculated as

$$\mathbf{d}\mathbf{d}^{\rm sec}(t,\mathbf{y}) = \frac{\overline{c}\_{s}^{\rm sec}(t,\mathbf{y}) - c\_{s,\rm init}}{c\_{s,\rm max} - c\_{s,\rm init}}.\tag{4.130}$$

The volume-averaged lithium concentration of the solid phase at time *t* at the position *y* is calculated as

$$
\overline{c}\_s^{\rm sec}(t, \mathbf{y}) = c\_s(t, \mathbf{y})\,,\tag{4.131}
$$

in case of the classical half-cell model, whereas this value is

$$\overline{c}\_s^{\rm sec}(t,\mathbf{y}) = \frac{3}{r\_{\rm prim}^{-3}} \int\_0^{r\_{\rm sec}} z^2 c\_s^{\rm prim}(t,\mathbf{y}) \,\mathrm{d}z. \tag{4.132}$$

in case of the hierarchically structured half-cell model.

Figure 4.9: Local depth of discharge along cell direction at different C rates. a) 3*C*. b) 5*C*. c) 7*C*.

In Figure 4.10, dodsec(*t*, *y*) is plotted versus both the normalized *x* and *y* position inside the secondary particle. The cathode position (*x* − *L* sep)/*L* pos = 0 refers to a secondary particle sitting at the separator-cathode interface, whereas (*x* − *L* sep)/*L* pos = 0.5 and (*x* − *L* sep)/*L* pos = 1 are secondary particles in the middle of the electrode and at the cathode-current-collector interface, respectively. Additionally, *y*/*r*sec = 1 denotes to the secondary particle surface, whereas *y*/*r*sec = 0 is the center. Again, the depth of discharge distribution is shown for the distinct C rates of 3*C*, 5*C* and 7*C*, respectively, and *t* is set as the simulation time at the cut-off voltage of 3V.

In general, the curves corresponding to n-NCM-F850 and n-NCM-F900 show less pronounced gradients as compared to the p-NCM curves. At 3*C*, it can be observed that the surfaces of the secondary particles in both classical and hierarchically structured half-cell model show high states of discharge at around 1.0. However, from the surface to the center of the secondary particles a more drastic decrease can be seen in case of p-NMC as compared to n-NMC-F850 and n-NMC-F900, where the latter two of these remain stable at a high level. At higher C rate of 5*C*, the differences in dodsec distribution in the secondary particles become more visible. While n-NMC-F850 and n-NCM-F900 show a rather constant distribution, a clear gradient can be observed for p-NCM. Additionally, following the observations made for Figure 4.9, the state of discharge distributions at the secondary particle surfaces along the cathode direction decline from the separator-cathode interface to the cathode-currentcollector interface. An interesting observation can be made for the case of 7*C*. Again, the gradients in the dodsec distributions inside the secondary particles are more pronounced for p-NCM than in the cases of n-NMC-F850 and n-NMC-F900. However, the depth of discharge at the secondary particle surface is actually higher at the separator-cathode interfaces in case of p-NMC than in the cases of n-NMC-F850 and n-NCM-F900.

From the above observations of the dodsec distribution, the conclusion from above can be extended. In addition to the enhanced utilization of the cathode's secondary particles due to an internal porosity, the depth of discharge distribution

Figure 4.10: Local depth of discharge along cell and secondary particle direction at different C rates. a) 3*C*. b) 5*C*. c) 7*C*.

and, thus, the concentration level of the active material is homogeneous. Even for larger C rates, the homogeneity remains stable inside the nanostructured secondary particles, which leads to a better rate performance as compared to the non-porous secondary particles.

As was reported in [30], the investigated electrodes n-NNC-F850, n-NMC-F900 and p-NMC show comparable loadings, i.e. active material mass per unit cross section. As was discussed above, the electrodes with the porous secondary particles show better utilization of the active material. This way, the energy density, i.e. energy per unit mass, is improved as well.

From the above observations, the superior performance of hierarchically structured electrodes as compared to classical electrodes can be reproduced by means of the presented models. This means, for instance, the hierarchically structured half-cell model can be used to estimate sensitivity of geometry, structural and material parameters, which is the topic in Section 5.3.

# **4.4 Summary**

In this chapter, a model for half-cells with hierarchically structured cathodes with porous secondary particles was proposed. As a first step, the well-established half-cell model based on the contributions of Newman and coworkers was recalled in detail. By means of the volume average method, the hierarchically structured half-cell model was proposed as a consistent extension of the classical half-cell model. The mathematical boundary conditions were acquired in detail and physically motivated for both models. Finally, the hierarchically structured half-cell model was qualitatively validated by comparing simulated and measured electrochemical performance based on real-world cathode structures taken from literature. The model successfully predicted the experimentally observed superior electrochemical performance of cathodes with porous secondary particles as compared to monolithic secondary particles. Additionally, the model allowed for the investigation of this observation by discussing the local lithium concentration distribution of the solid phase of the electrode and secondary particle level, respectively. It was shown that nanostructured secondary particles differ from classical ones by having more of a homogeneous concentration distribution for higher C rates. This way, the available active material capacity could be better exploited leading to a higher performance.

# **5 Applications**

In the following, the applicability of the in Sections 3 and 4 presented methods is demonstrated. On the one hand, the resistor network method is employed to derive prediction formulas for the effective resistance of porous secondary particles and for the effective transport properties of sphere packings. In both cases, transport is considered via the surface and the volume of the particles. On the other hand, the hierarchically structured half-cell model is used for parameter studies in order to reveal the influence of different structural and material parameters on the electrochemical performance of cathodes with porous active material particles.

# **5.1 Effective transport properties of hierarchically structured spheres**

In Section 3.1, it was shown that the resistor network approach can be utilized to compute effective transport properties of sphere assemblies. As for classical cell models of the type presented in Section 4.1, for instance, this approach can be used to deliver effective transport properties of the porous electrode structure. That is to say, if the electrode is modeled as sphere packings. Concerning the hierarchically structured cell model from Section 4.2, the active material is modeled as spherical secondary particles built up by smaller spherical primary particles. In this case, in order to model effective transport properties through the active material phase, certain parts of the resistor network method need to be adjusted accordingly.

Consider the exemplary network of porous secondary particles, in Figure 5.1a. In Figure 5.1b, the equivalent node resistor network can be observed, where the secondary particle centers are the nodes and the contacts between the particles are identified as the resistors. In Section 3.1, it was pointed out that the crucial part in the framework of the RN is the computation of the resistances of the individual contact pairs. The most practical way is to have analytical formulas as in Equations (3.3) and (3.11), where the resistances can be evaluated based on the geometry of the overlapping particles.

Therefore, the following section aims at deriving a formula describing the resistance between contact pairs of hierarchically structured spheres by their geometry and inner structure. Moreover, two cases will be considered. In the first case, the transport is through the volume of the primary particles and, in the second case, the transport is via the surfaces of the primary particles.

Figure 5.1: Resistor network method for hierarchically structured spheres. a) Assembly of hierarchically structured spheres. b) Equivalent network with nodes at the centers of the spheres and resistances assigned to the contact pairs.

## **Model**

To begin with, similar to Equation (3.11), the resistance between two contacting hierarchically structured spheres is modeled as a series connection of two resistances. Moreover, those two resistances are computed by utilizing a model where each of the spheres are put between two conducting plates. In Figure 5.2, the replacement model of a hierarchically structured sphere between two conducting plates is sketched. The larger secondary particle comprises of smaller primary particles. The respective primary particle and secondary particle radii are *r* prim and *r*sec, respectively. Furthermore, the secondary particle is positioned between two conducting plates with a gap distance of *h* sec, which results in a contact radius of *r*sec,*<sup>c</sup>* . The internal contact radii *r I*,*J* prim,*c* are between the primary particles with radii *r I* prim and *r J* prim. Clearly, the effective resistance of the model in Figure 5.2 is depends on all the above presented structural properties. A formula representing effective resistances of hierarchically structured spheres must therefore contain all those quantities.

Figure 5.2: Model of hierarchically structured sphere between two conducting plates.

The approach chosen here was to conduct a variety of numerical experiments and deduce structural dependencies. The resistor network method was chosen to perform the numerical experiments. However, before running parametric studies, the applicability of the RN on hierarchically structured secondary particles with inner porosity has to be verified.

# **Verification of the RN for hierarchically structured spheres**

Similar to rectangular domains in Section 3.1, hierarchically structured spheres were generated, the structure was densified and the effective resistance was calculated using the resistor network method. Finally, the results were compared to finite element simulations. Concerning the generation of the assemblies, the RCP was modified to create randomly generated, overlap-free and densely packed collections of spheres inside a spherical domain. The initial structures were then further densified using the numerical sintering algorithm. In this case, the primary particles were fixed in space and their radii were successively increased until a certain densification criterion was reached. This way, the spherical shape of the secondary particle was maintained.

Importing the data, setting up the models and evaluating the effective transport properties for both cases of volume and surface transport followed the descriptions from Sections 3.1.4 and 3.1.2, respectively. In order to model the conducting plates, in case of the RN, sphere caps of a certain height were removed from the secondary sphere, see Figure 5.3a, on two opposing sides in transport direction, see Figure 5.3b. Consequently, all primary particles, with their centers lying inside those virtual caps, were deleted, see the highlighted spheres in Figure 5.3b. Finally, the bottoms of the caps were chosen to be the conducting plates connected to the intersecting primary particles, see Figure 5.3c. This is where the boundary nodes were placed following the description in Section 3.1.

Figure 5.3: Hierarchically structured sphere between two conducting plates.

In the finite element model, the primary particles were simply cut off at a defined position in the transport direction. This can be seen in Figure 5.4b and d. Concerning the FEM models, all nodes on one of the conducting plate surfaces were set to the same temperature, while the potential of the nodes on the opposite surface were dropped by ∆*T* = 1 with respect to this potential. Due to the mathematical equivalence of the transport problems considered in this work, the methods used for heat transport problems solved by finite element methods can be compared equivalently to electric transport problems solved by the resistor network method. Therefore, concerning the RN, an electric potential drop of ∆ϕ = 1 was imposed on the boundary nodes. Analogous to Equations (3.5) and (3.7), the effective resistance of a secondary particle was computed by

$$R^{\rm FEM} = \frac{\Delta T}{\phi\_{\rm q}^{\rm FEM}} \qquad \text{and} \quad R^{\rm RN} = \frac{\Delta \phi}{I\_{\rm eff}}, \tag{5.1}$$

respectively. In the FEM model, the resulting heat flux φ FEM <sup>q</sup> was obtained at one of the surfaces with the applied temperature and, in case of the RN, the effective current *I* eff was calculated using the resistor network solving scheme according to Section 2.3.

Figure 5.4: Verification models for the effective conductivity of hierarchically structured spheres using the resistor network and finite element method. Top: Transport via the surface, i.e. solid-surface. Bottom: Transport via the volume, i.e. solid-volume

### **Test cases**

In the following, 15 scenarios of assemblies were created, according to the method described above. The scenarios were divided into three types of assemblies, see Table 5.1. For all assembly types, the mean radius *r*mean and the mean contact angle θmean were fixed to 1 and 15°, respectively. The only varying parameter was the standard deviation *r*<sup>σ</sup> ranging from 0 to 0.25 such that, with increasing assembly type number, the polydispersity increased. The resulting volume fraction of the solid phase φsolid,mean and pore phase φpore,mean as well as the mean contact radius *r*c,mean is presented in Table 5.1 as well. In addition, the transport is considered either through the volume or the surface of the primary particles. In the latter case, the surface is modeled by shells with a thickness of 0.1 ·*r* prim, i.e. 10% of the primary particle radii.

Table 5.1: Structural parameters of secondary particles for the verification of the transport through the volume and via the surface of solid spheres.


### **Evaluation and discussion**

The results of this investigation can be observed in Figure 5.5. Effective resistances *R*ˆRN eff and *<sup>R</sup>*ˆFEM eff are calculated by Equations (3.8) and (3.17), respectively, and normalized by the bulk or shell conductivity. On the horizontal axis the results provided by the FEM analysis are shown and on the vertical axis the corresponding values by the RN are drawn. The black solid line in this figure represents a perfect match of the two results whereas values above or below indicate an over- or underestimation of the effective transport property by the RN, respectively.

Figure 5.5: Verification of the resistor network method for hierarchically structured spheres using the finite element method. a) Transport through the volume of primary particles. b) Transport via the surface of primary particles.

It can be seen that, while the mean error is around *e* = 6% in the case of surface transport, the mean error is almost zero in the case of volume transport. Therefore, given the considerably low errors, the above observations allow the conclusion that the resistor network approach is suitable to compute the effective transport properties of hierarchically structured spheres between two conducting plates.

### **Parametric studies and prediction formulas**

In order to model the effective resistance, the prediction formula

$$R\_{\rm eff} = \hat{\mathcal{R}}^{\rm fit} \cdot \frac{r\_{\rm prim,m}}{r\_{\rm prim,c,m}} \cdot \frac{r\_{\rm sec}}{r\_{\rm prim,m}} \cdot \frac{r\_{\rm prim,m}}{s\_m} \cdot \frac{1}{r\_{\rm sec}} \cdot \rho\_{\rm bulk,shell} \tag{5.2}$$

is proposed, where *R*ˆfit is a fit function and the rest of the parameters are interpreted as follows. Here, the first model parameter is the ratio *r* prim,*m*/*r* prim,*c*,*m* , which is interpreted as the degree of densification of the secondary particle and relates the mean primary particle radius *r* prim,*m* to the mean contact radius of the primary particles *r* prim,*c*,*m* . Furthermore, *r*sec/*r* prim,*m* describes the tortuosity effect and takes the secondary particle radius *r*sec into account. Finally, *<sup>r</sup>*sec,*c*/*r*sec is the ratio of the contact radius of the conducting plates *r*sec,*<sup>c</sup>* to the secondary particle radius and ρbulk,shell is either the bulk or the shell, resistivity, depending on if volume or surface transport is considered. Additionally, *sm*/*r* prim,*m* relates the shell thickness *s<sup>m</sup>* to the mean primary particle radius, where, in case of volume transport, this parameter is omitted.

The degree of densification *r* prim,*m*/*r* prim,*c*,*m* can theoretically range within the limits of 1 and ∞. In the first case, the contact radii are of the same size as the primary particles, which means that the structure can be regarded as bulk material. The second case represents a structure where the primary particles have no contact to each other, which results in no conducting pathways. The tortuosity effect described by *r*sec/*r* prim,*m* is within the limits of 1 and ∞. Roughly speaking, the closer the primary and secondary particle sizes are, the shorter the conducting paths within the secondary particle. As a consequence, the effective resistance decreases. On the other hand, the effective resistance increases if comparably small primary particles form more tortuous and longer pathways. Finally, the values of *r* prim,*m*/*s<sup>m</sup>* are within the limits of 1 and ∞. The latter of which resembles a mean shell thickness of 0 and no transport paths via the surface. Consequently, the effective resistance tends to infinity as well. In theory, a shell thickness of the size of *s<sup>m</sup>* provides the best transport paths via the surface and, therefore, the lowest effective resistance.

In the following, the verification of the prediction formula in Equation (5.2) is demonstrated by 100 randomly generating hierarchical sphere structures following the tool chain of RCP (modified to account for spherical domains), numerical sintering and RN, see Section 3.1. For each of the scenarios, the combination of the model parameters were chosen at random within the limits shown in Table 5.3.

Table 5.3: Variation of parameters of the hierarchically structured spheres.


Figure 5.6 shows the results of the numerical experiments. In Figure 5.6a, transport through the volume and, in Figure 5.6b, transport via the surface is shown. The effective resistance is plotted over the dimensionless ratio *r*sec,*c*/*r* prim,*m* . It can be observed that, on the one hand, the results can be represented in a dimensionless form using the parameters in Table 5.3. On the other hand, the curves approximately lie on a simple rational 1/*x*-type function, which only depends on the ratio *r*sec,*c*/*r* prim,*m* . This function can be given by

$$
\hat{\mathcal{R}}^{\text{fit}} = \frac{1}{a\_{\text{fit}} \frac{r\_{\text{scc},c}}{r\_{\text{prim},w}}} \,, \tag{5.3}
$$

where *a*fit ≈ 6, in case of volume transport, and *a*fit ≈ 20, in case of surface transport.

Using the prediction formula from Equation (5.2), the behavior of the effective resistance in terms of the structural properties can be investigated. For instance, in both the volume and the surface transport case, the effective resistance tends to increase to rather large numbers if *r*sec,*c*/*r* prim,*m* approaches small values. This means that, given a constant *r*sec,*<sup>c</sup>* , for increasing *r* prim,*m* , the number of primary particles in contact with the conducting plate decreases leading to higher effective resistances. Vice versa, an increasing ratio of *r*sec,*c*/*r* prim,*m* leads to a smaller effective resistance, which means that more primary particles are in contact with the conducting plates.

Figure 5.6: Prediction formulas to estimate the effective resistance of hierarchically structured spheres between two conducting plates. a) Transport is via the volume of the primary particles. b) Transport is via the surface of the primary particles.

## **Conclusion**

Equation (5.2) represents the effective resistance of a hierarchically structured secondary sphere between two plates. The secondary sphere is composed of smaller primary spheres. The formula can be used to model transport either through the volume or via the surface of the primary spheres. Moreover, it can be used to model transport between two contacting hierarchically structured secondary particles. Practically, and similar to Equations (3.3) and (3.11), the fit formula can be used in the framework of the resistor network method to calculate the effective conductivity of an assembly of hierarchically structured spheres.

# **5.2 Effective volume and surface transport properties of sphere assemblies**

The knowledge of effective transport properties in granular materials such as battery electrodes [18, 150] plays an important role for battery modeling. Usually, effective transport properties of granular materials are approximated by the volume fraction of the transport phase under consideration, see [69] for a review. Battery modelers tend to use a Bruggeman type relation [71, 138] to estimate effective transport properties of the porous electrodes [36, 41, 46]. In [92], the resistor network method was used to show that this type of empirical relation works well for the pore phase of sphere assemblies with a polydisperse size-distribution. In order to model effective transport properties of particle networks, however, the Bruggeman relation is not well suited [67].

In order to overcome this drawback, adequate prediction formulas need to be found. Therefore, in the following, the practical use of the solid-volume model from Section 3.1.2 and the solid-surface model from Section 3.1.4 is demonstrated. The models are used to develop prediction formulas to compute effective transport properties of assemblies of overlapping and equal-sized spheres, where the transport can be through the volume or the surface of the spheres.

## **Model parametrization**

A series of 100 scenarios of randomly distributed and monosized sphere assemblies were created in a box-shaped simulation domain using a certain set of model parameters. The model parameters were chosen as follows. A porosity range φpore between 0.1 and 0.4 was targeted representing, on the one hand, a rather dense packing and, on the other hand, one which is near the percolation threshold of monosized spheres [76]. Moreover, the sphere radii *r<sup>p</sup>* were between 0.1 and 0.8 representing an arbitrary but rather large range of particle sizes. This particle size range is comparable to the primary particle sizes possible in [30], which ranged between 0.179 and 0.6 µm. Note, however, an arbitrary length unit lu is chosen, since the investigated transport parameters are dimensionless and independent of the chosen length unit. Finally, in case of surface transport, a constant shell thickness was set to 0.01.

Analogous to Section 3.1.2, the initial sphere packings were generated using the RCP. The structures were then further densified by means of the numerical sintering algorithm. Additionally, the numerical sintering algorithm was slightly modified. In contrast to the algorithm described in Section 3.1.2, the boxshaped simulation domain was isotropically and successively reduced such that the spheres inside the domain were moved upon each other. This was done until a densification criterion, i.e. the targeted porosity, was achieved. In this way, the particle sizes could be preserved.

## **Results and discussion**

The effective transport parameters through the volume and via the surface of the spheres were calculated using the resistor network method as presented in Section 3.1.2 and Section 3.1.4, respectively, and evaluated in Figure 5.7. Here the effective transport parameters are drawn over the porosities and particle radii. In Figure 5.7a, the effective transport parameters of the volume transport case is shown whereas, in Figure 5.7b, the surface transport case can be observed. Recall from the above mentioned sections that the effective transport parameters are defined as the effective transport properties normalized by the transport coefficients of the particle volumes or shells. In the present investigation, without loss of generality, they both set to unity.

Apparently, in Figure 5.7a, the effective volume transport parameters only depend on the porosity of the system. For high porosity values the effective transport parameters are low while for low porosities they increase. Therefore, the numerical results are fitted by the bi-linear surface fit

$$
\hat{k}\_{\rm vol}^{\rm fit} = a\_0 + a\_1 \phi\_{\rm pore} + a\_2 r\_{\rm p} \,, \tag{5.4}
$$

where the parameter values are *a*<sup>0</sup> = 1.0, *a*<sup>1</sup> = −2.130 and *a*<sup>2</sup> = 0.01. Here, *a*<sup>0</sup> is chosen such that, for the special case of a completely dense structure, i.e. if φpore → 0 and *r*<sup>p</sup> → 0, the effective transport parameter becomes 1 resembling bulk material properties.

In Figure 5.7b, the effective transport parameter via the surface is shown. Remarkably, the effective transport parameter is not only dependent on the porosity but also on the particle radius. The best effective transport parameter is achieved for low particle radii and low porosities. Additionally, it can be observed that for a constant porosity, the effective transport parameter can be enhanced by decreasing the particle radius. Finally, the resulting values where fitted by

$$\hat{k}\_{\text{surf}}^{\text{fit}} = \frac{1}{(1 + \phi\_{\text{pore}})^{b\_0}} \cdot \frac{1}{(1 + r\_{\text{p}})^{b\_1}} \,, \tag{5.5}$$

where the fitting parameter values are *b*<sup>0</sup> = 6.717 and *b*<sup>1</sup> = 5.332. Note that the function was chosen such that, for the completely dense structure case, i.e. φpore → 0 and *r*<sup>p</sup> → 0, Equation (5.5) becomes bulk material property.

Figure 5.7: Evaluation of transport parameters by the parametric study. a) Transport through volume. b) Transport via surface. c) Comparison of surface and volume transport over the porosity.

Consider the special case of linearly increasing particle radius and decreasing porosity. A similar case was observed in [151]. There, pellets were investigated made of conductive particles which were pressed together and sintered afterwards. By adjusting the pressure and sinter temperature, different particle radii and the pellet porosities were achieved leading to increasing particle radii while, at the same time, decreasing porosity. In Figures 5.7a and 5.7b, such a path is indicated by ˆ*k* fit,path vol or <sup>ˆ</sup>*<sup>k</sup>* fit,path surf in case of volume or surface transport, respectively. In Figure 5.7c, Equations (5.4) and (5.5) are evaluated along this path. In particular, the curves are plotted over porosity. While ˆ*k* fit,path vol decreases with increasing porosity in the volume transport case, which can be expected, the surface transport curve ˆ*k* fit,path surf rises. Clearly, <sup>ˆ</sup>*<sup>k</sup>* fit,path vol decreases because of reducing overlaps of particles. However, reducing overlaps apply also to ˆ*k* fit,path surf . Therefore, this effect seems to be compensated by decreasing particle sizes because, due to that, the specific surface area increases as well.

Apparently, the specific surface area curve along the given path is similar to the evolution of the transport parameter via the surface. Given that the specific surface area is intrinsically accounted for by the transport angle of the resistance in Equation (3.11), this directly effects the overall effective conductivity.

Figure 5.8: Computation of specific surface area for assemblies of spheres.

The specific surface area is defined as the total free surface area divided by the simulation box volume. The computation of the specific surface area—as used in this work—is sketched in Figure 5.8. For a given assembly of spheres depicted on the left-hand side, the total free surface area is calculated by summing up the surface areas of each individual particle. The magnification on the right-hand side exemplarily highlights a single particle out of the assembly. The individual free surface area is computed by subtracting the surface area of the sphere caps, which are formed by the overlaps with the neighboring particles, from the total surface area of the particle.

In Figure 5.9, the evolution of specific surface area for the present investigation can be observed. In Figure 5.9a, the resulting values are fitted by

$$a\_s^{\rm fit} = \frac{c\_0}{(1 + \phi\_{\rm pore})^{b\_1}} \cdot \frac{c\_2}{(1 + r\_\mathbf{p})^{c\_3}},\tag{5.6}$$

where the fitting parameters are *c*<sup>0</sup> = 3.636, *c*<sup>1</sup> = −1.104, *c*<sup>2</sup> = 5.196 and *c*<sup>3</sup> = 5.397. The evolution of specific surface area along the given path *a* fit,path *s* can be seen in Figure 5.9b.

Figure 5.9: Specific surface area by the parametric study.

### **Conclusion**

The investigation results reveal that Equations (5.4) and (5.5) can be used to model the effective transport parameters via volume and surface of mono-sized sphere assemblies with varying porosity and particle sizes. In case of cell modeling, Equation (5.4) can be used to calculate the transport properties of the active material phase, given that it is modeled by spheres. Additionally, since in Equation (5.5) thin layers of shells are assumed, it can be utilized to model effective transport properties of high-conductive carbon coatings of active material. Finally, the above findings support the observations in [151]. There, an increase in effective ionic conductivity was observed in case of larger porosities. This increase was attributed to the dominant transport via the surfaces of the particles due to an increase in surface area.

# **5.3 Parameter study on hierarchically structured electrode performance**

In Section 4.3, the hierarchically structured half-cell model, as proposed in this work, was validated qualitatively by experiments taken from literature. It was concluded that it can be used to investigate the influence of the electrode's geometry, structure and material on the cell's electrochemical performance. By means of this model, the difference between classical and hierarchically structured cathodes in terms of sensitivity to the diffusion coefficient and electronic conductivity of the active material was investigated in [137]. In the following study, however, the focus is on the hierarchically structured cathodes only and a larger set of parameters is taken into account. The influence of diffusion coefficient and electronic conductivity of the active material as well as secondary particle morphology on the rate capability is analyzed. To this end, the cathode structure n-NMC-F900 from Section 4.3 was used as a reference and the mentioned parameters were varied as follows. Note that, if not stated otherwise, the other parameters concerning geometry, structure, transport and material were taken from Tables 4.1, 4.2, 4.3 and 4.4.

The electronic conductivity was assumed between 1 · 10−<sup>3</sup> , 1 · 10−<sup>4</sup> and 1 · 10−<sup>5</sup> Sm−<sup>1</sup> representing a rather large spectrum of high and low values and corresponding to lithiation states of γ = 0.9, 0.95 and 1.0 in LiγNi1/3Mn1/3Co1/3O2, see Equation (4.124). The diffusion coefficient was set as 1 · 10−13, 1 · 10−<sup>14</sup> and 1 · 10−<sup>15</sup> m<sup>2</sup> s −1 , which are the extreme values found in literature [142–146]. The secondary particle morphology is characterized by the porosity, primary particle radii, specific surface area and effective transport parameter. In this investigation, the primary particle radii and effective transport parameters were varied while keeping the porosity fixed. Note that, indirectly, the specific surface area is influenced by changing the radii since it is calculated as *a* = 3φ *r* [136]. In order to study the influence of the transport via the primary particle network, the effective transport parameter ˆ*k* eon,sec eff is varied between values of 0.25, 0.50 and 0.75, where the first and the middle value is comparable to the lowest and largest parameters which were found for n-NMC-F850 and n-NMC-F900, respectively. An additional effective transport parameter of 0.75 is chosen to see if an improvement of transport path quality would influence performance.

For all 81 combinations of the above described parameters, the rate performance was simulated by using C rates from 1/20 to 10 and the simulation was stopped at a cut-off voltage of 3V. Results of this investigation can be found in Figure 5.10, where some representative cases are selected to be shown. From the top to the bottom row, the electronic conductivity is 1 · 10−<sup>3</sup> , 1 · 10−<sup>4</sup> and 1 · 10−<sup>5</sup> Sm−<sup>1</sup> , respectively. From the left to the right column, the transport coefficient, the primary particle radii and the diffusion coefficient is varied.

In general, it can be observed that from the top to the bottom the rate capability decreases for decreasing electronic conductivity. While a distinct specific capacity drop in case of the largest two conductivity values can be observed after 3*C*, the performance drops after 1*C* for the case of lowest conductivity. Note, however, that the gradient of the drop is steeper in case of the larger two conductivity values. On the other hand, for the second and third column, the variation of particle radii and diffusion coefficient does not seem to influence the shape of the performance curve as the curves within one diagram are almost identical. However, in case of the first column, where the transport coefficient was varied, a clear difference between the curves can be observed for decreasing conductivity of the active material. Especially, in the lower left diagram, it can be seen that for increasing values of ˆ*k* eon,sec eff , the rate performance decreases more visibly.

From the above observations, some conclusions can be drawn. Varying the diffusion coefficient within the range presented here does not seem to influence the performance. While it was found in [137] that for classical cathodes the rate-limiting factor is the diffusion of lithium into the monolithic active material particles, in case of the hierarchically structured cathodes, sensitivity to the diffusion coefficient could not be observed, see also [53]. In addition to that, changing the primary particle radii, in the range presented here, has no visible effect on the performance. This indicates that the primary particle radii in combination with the diffusion coefficient still provide diffusion paths which are short enough to not have an impact. However, electronic conductivity of the active material has a strong effect on rate capability. On the one hand, this can be explained by assuming that lithium diffusion via the electrolyte phase of the secondary particle pores is favored. Here, the diffusion coefficients are orders of magnitudes higher than in the solid phase [146, 152]. On the other hand, low electronic conductivity values lead to a lower performance, which was also observed in [53, 152]. Therefore, electronic conductivity becomes ratelimiting. In addition to this observation, this investigation revealed the secondary particle morphology plays a crucial role when the electronic conductivity reaches lower values. In other words, poor inter-primary particle contact becomes more problematic in case of low electronic conductivity, which was also observed in [152].

Figure 5.10: Parameter study rate performance. Varying electronic conductivity of active material from top to bottom row. Varying transport coefficient, primary particle radii and diffusion coefficient from left to right column.

The effect of the transport coefficient on the electrochemical performance is further investigated with the help of Figure 5.11, where, analogous to Section 4.3, the depth of discharge distribution along the secondary particle radii is shown. The lowest electronic conductivity of 1 · 10−<sup>5</sup> Sm−<sup>1</sup> is chosen at a C rate of 3*C* as, in this case, the differences in performance are visible the most. Three cases of transport coefficients are discussed, where it can be observed that for all cases the value of dodsec is around 1 at the surfaces of the secondary particles. From the surface to the particle center, the depth of discharge decreases slowly at first but at a certain point shows a distinct drop and decreases slowly again until the end. In case of higher values of ˆ*k* eon,sec eff , this distinct drop occurs more towards the surface of the secondary particles as well as results in lower dodsec towards the center. This leads to a lower depth of discharge DOD for the whole cell represented by lower retained specific capacity.

Figure 5.11: Depth of discharge distribution along secondary particle radius direction at 3*C*.

Poor inter-particle contact and low electronic conductivity of the active material can explain the observations in Figure 5.12, which is based on investigations in [137] and where the two cathode structures n-NMC-F850 and n-NCM-F900 are compared to each other. Experiments in Figure 4.8b show that for higher C rates, the two rate performance curves diverge from each other. Notably, however, n-NCM-F900 performs a little better than n-NCM-F850 even though the latter of which shows more preferable structural properties in terms of higher specific surface area of primary particles for electrochemical reactions inside the secondary particles and smaller primary as well as secondary particles leading to shorter diffusion paths. For the simulation, the rate performance for the two structures were calculated using the geometry, structure and material parameters from Tables 4.1, 4.2, 4.3 and 4.4. However, the secondary particles' electronic conductivity was set to the lowest value of 1 · 10−<sup>5</sup> Sm−<sup>1</sup> in both cases. This low value could be justified as a first simple measure to account for additional resistance effects such as contact resistances between primary particles beyond the purely geometrical bottleneck effects considered so far. The results show a clear difference in performance, where n-NMC-F900 achieved better rate capability than n-NMC-F850. By comparing the effective transport parameters from Table 4.3 of the solid phase for both cases, it can be seen that the value of ˆ*k* eon,sec eff = 0.46 in case of n-NMC-F900 is more than twice as large as it is in case of n-NMC-F850, which is 0.20. As was pointed out before, effective transport parameters reflect the quality of transport paths which, ultimately, characterizes the connectivity of primary particles. Adding to the fact that electronic transport is predominantly via the primary particle network of the secondary particles [153], the combination of low electronic conductivity and poor connectivity of primary particles may lead to low electrochemical performance. Therefore, in view of the above observations, the rate-limiting factor is the effective electronic conductivity of the primary particle network inside the porous secondary particles.

Figure 5.12: Simulation results based on real electrode structure and using low electronic conductivity.

## **Conclusion**

A large scale parameter study was conducted by means of the qualitatively validated hierarchically structured half-cell model from Section 4.3. The influence of electronic conductivity, transport and diffusion coefficients of the active material as well as the primary particle radii on the rate performance was discussed. It was found that diffusion coefficients and primary particle sizes have barely visible effect as the diffusion paths might still be short enough to not have an impact on performance. However, it could be shown that low electronic conductivity and large transport coefficient influence the rate capability more significantly. It was concluded that the combination of low electronic conductivity and poor inter-connectivity of the primary particle network inside the secondary particles might become rate-limiting in case of hierarchically structured cathodes.

# **6 Summary and outlook**

In the present work, the resistor network method was extended with regard to the transport through the solid and the pore phase of granular systems represented by sphere packings with polydisperse size-distributions. As for the solid phase, transport through the volume, via the surface or a mix of both were considered. For all cases, appropriate analytically derived formulas from literature describing resistance between two single particles were used or combined accordingly. Finally, these single contact models were embedded into the framework of the resistor network method in order to compute effective transport properties. All the proposed models—single contact as well es effective transport models—were verified using finite element methods.

Regarding the pore phase, a novel method concerning the computation of transport through the pore phase was developed. The key was the so-called Laguerre tessellation, where the pore phase of the system is decomposed into cells with each of them containing a single particle. Consequently, the cell nodes and edges were the basis of equivalent resistor networks. The nodes were identified as pore centers and the edges were the pore throats. As an extension, this model was modified to account for more than one conducting species in the pore phase. Both methods were either verified using the finite element method or validated with the help of experiments taken from literature.

In the application part of this work, it was shown how the resistor network method for volume and surface transport can be used to conduct large scale parameter studies. Due to the resource and time efficiency of this method, a large variety of different structural compositions can be simulated in a short amount of time. In this work, the resulting database allowed for the derivation of prediction formulas for the effective resistance of porous secondary particles and porous electrodes consisting of mono-sized and overlapping spheres.

A model for half-cells with hierarchically structured cathodes with porous secondary particles was proposed. As a preparation step, the well-established half-cell model in accordance to the contributions of Newman and coworkers was recalled in detail first. It was pointed out that the mathematical framework of this model can be derived equivalently by employing volume averaging methods. By means of the same methods, the hierarchically structured halfcell model was proposed as a consistent extension of the classical half-cell model. For both models, the mathematical boundary conditions were acquired in detail and they were physically motivated. The hierarchically structured halfcell model was applied to real-world cathode structures taken from literature. In doing so, the model was validated qualitatively by comparing simulated and measured electrochemical performance as it successfully predicted the experimentally observed superior electrochemical performance of cathodes with porous secondary particles as compared to monolithic secondary particles. Additionally, it allowed for the investigation of this observation by discussing the local lithium concentration distribution of the solid phase along electrode and secondary particle direction. It was shown that nano-structured secondary particles differ from classical ones by having more of a homogeneous concentration distribution for higher C rates. This way, the available active material capacity could be better exploited leading to a higher performance. Finally, the hierarchically structured half-cell model was applied for parameter studies. By varying the diffusion coefficient and electronic conductivity of the active material as well as the secondary particle morphology, it was found that the combination of low conductivity and poor inter-connectivity of the primary particle network inside the secondary particles might become rate-limiting.

As a possible future work, the resource and time efficiency of the resistor network as well as the presented half-cell models can be used for further investigations. A large variety of artificial particle-based structures can be generated and effective transport properties for the solid and pore phase can be computed representing the solid and the electrolyte phase in the cell models. These parameters can be incorporated into the cell models in order to investigate the impact on the electrochemical performance. This way, the electrochemical performance for a large set of parameters regarding the geometry, structure and material properties can be computed efficiently and, ultimately, the best combination can be chosen to design optimized electrodes.

# **Bibliography**


films as anodes for lithium ion batteries," *Journal of Power Sources*, vol. 188, no. 2, pp. 588–591, 2009.


Transport in Single NCM Cathode Active Material Particles for Lithium-Ion Batteries Studied under Well-Defined Contact Conditions," *ACS Energy Letters*, vol. 4, no. 9, pp. 2117–2123, 2019.

[154] W. G. Schlecht, "Calculation of density from x-ray data," *American Mineralogist*, vol. 29, no. 3-4, pp. 108–110, 1944.

# **List of symbols**

#### Constants


### Physics


### Physics cont.


### Chemistry


### Chemistry cont.


### Resistor network method


### Resistor network method cont.


### Cell model


### Cell model cont.


### Cell model cont.


### Cell model cont.


### Sub- and superscripts


### Sub- and superscripts cont.


#### Acronyms


### Acronyms cont.


# **A Mathematical operations**

## **A.1 Divergence operator**

The divergence operator in the cartesian coordinate system is

$$\nabla \cdot \vec{F} = \frac{\partial F\_{\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial F\_{\mathbf{y}}}{\partial \mathbf{y}} + \frac{\partial F\_{\mathbf{z}}}{\partial \mathbf{z}} \quad , \tag{A.1}$$

where *x*, *y*,*z* are the cartesian coordinates and *Fx*,*Fy*,*F<sup>z</sup>* are the respective components of the flux vector.

The divergence operator in the spherical coordinate system is

$$\nabla \cdot \vec{F} = \frac{1}{r^2} \frac{\partial 1}{\partial r} \left( r^2 F\_r \right) + \frac{1}{r \sin(\Theta)} \frac{\partial}{\partial \Theta} \left( \sin(\Theta F\_\Theta) \right) + \frac{1}{r \sin(\Theta)} \frac{\partial}{\partial \Phi} F\_\Phi \quad , \quad (A.2)$$

where *r*,Θ,Φ are the spherical coordinates and *F<sup>r</sup>* ,*F*Θ,*F*<sup>Θ</sup> are the respective components of the flux vector.

# **A.2 Integration**

## **Cartesian domains**

The volume integral of a function *f*(*x*, *y*,*z*) is

$$\iiint\_{V} f(\mathbf{x}, \mathbf{y}, \mathbf{z}) \, \mathrm{d}V = \int\_{0}^{L^{\mathbf{x}}} \int\_{0}^{L^{\mathbf{y}}} \int\_{0}^{L^{\mathbf{z}}} f(\mathbf{x}, \mathbf{y}, \mathbf{z}) \, \mathrm{d}\mathbf{x} \, \mathrm{d}\mathbf{y} \, \mathrm{d}\mathbf{z} \quad , \tag{A.3}$$

where *Lx*,*Ly*,*L<sup>z</sup>* is the edge length of a box. In case of the function being independent of *y*, and *z*, the volume integral is

$$\int\_{0}^{L^{\chi}} \int\_{0}^{L^{\chi}} \int\_{0}^{L^{\chi}} f(\mathbf{x}) \, \mathbf{dx} \, \mathbf{dy} \, \mathbf{dz} = L^{\chi} L^{\varepsilon} \int\_{0}^{L^{\chi}} f(\mathbf{x}) \, \mathbf{dx} \tag{A.4}$$

and the volume average is

$$\frac{L^{\chi}L^{z}\int\_{0}^{L^{x}}f(\mathbf{x})\,\mathrm{d}\mathbf{x}}{V} = \frac{1}{L^{x}}\int\_{0}^{L^{x}}f(\mathbf{x})\,\mathrm{d}\mathbf{x} \quad . \tag{A.5}$$

## **Spherical domains**

The volume integral of a function *f*(*r*,Θ,Φ) is

$$\int\_{0}^{2\pi} \int\_{0}^{\pi} \int\_{0}^{R} f(r, \Theta, \Phi) r^2 \sin(\Theta) \,\mathrm{d}r \,\mathrm{d}\Theta \,\mathrm{d}\Phi \quad . \tag{A.6}$$

In case of spherically symmetric conditions, the angular dependencies are constant such that the concerning volume integrals reduce to

$$4\pi \int\_{r=0}^{\infty} f(r)r^2 \,\mathrm{d}r \quad . \tag{A.7}$$

Furthermore, the volume average of the above equation is

$$\frac{4\pi \int\_{r=0}^{\infty} f(r)r^2 \,\mathrm{d}r}{^{4/3}\pi R^3} = \frac{3}{R^3} \int\_{r=0}^{\infty} f(r)r^2 \,\mathrm{d}r \quad . \tag{A.8}$$

# **B Physics**

## **B.1 Specific capacity from stoichiometry**

By using Faraday's law, which states that the rate of production of a species is proportional to the current, and the total mass produced is proportional to the amount of charge passed multiplied by the equivalent weight of the species [47]:

$$m\_i = -\frac{s\_i M\_i \text{It}}{n \mathcal{F}}\,,\tag{\text{B.1}}$$

where *m<sup>i</sup>* is the mass of species *i* produced by a reaction in which its stoichiometric coefficient is *s<sup>i</sup>* and *n* electrons are transferred, *M<sup>i</sup>* is the molar mass, F is Faraday's constant, and the total amount of charge passed is equal to the current *I* multiplied by time *t*. Rearranging Equation (B.1) brings

$$\underbrace{\frac{It}{m\_i}}\_{\mathcal{Q}\_{i,\text{spec}}} = -\frac{n\mathcal{F}}{s\_i\mathcal{M}\_i} \tag{B.2}$$

which yields the specific capacity of species *i* as *Qi*,spec = *It mi* , where *It* = *Q<sup>i</sup>* is taken as the total capacity of species *i*.

The above mentioned coefficients of species *i* refer to the electrode reaction

$$\rm{s-}M\_{-}^{\;\!\;-} + \rm{s}\_{+}M\_{+}^{\;\!\;+} + \rm{s}\_{0}M\_{0} \rightleftharpoons \rm{n}\mathbf{e}^{-},\tag{\bf B.3}$$

where *M* is the symbol for the chemical formula of species *i*, and *z*−, *z*<sup>+</sup> are signed charge numbers of cation and anion [10].

As for the positive electrode, the reaction in a lithium-ion cell is

$$\rm{M} + \rm{Li}^+ + e^- \rightleftharpoons \rm{LiM},\tag{B.4}$$

where the M stands for some metal oxide in the positive electrode material [10]. Rearranging Equation (B.4) into the form of Equation (B.3) as

$$\text{LiM}-\text{M}-\text{Li}^+ \rightleftharpoons \text{e}^-,\tag{\text{B.5}}$$

leads to the coefficients *z*<sup>−</sup> = *z*<sup>+</sup> = 1, *s*<sup>−</sup> = 0, *s*<sup>+</sup> = −1, *s*<sup>0</sup> = 0 and *n* = 1.

In the following, Equation (B.6) is used to calculate the specific capacity of the positive electrode active material LiNi1/3Mn1/3Co1/3O2, which is used in this work. Therefore, the molar mass *M<sup>i</sup>* ≡ *M*LiNi1/3Mn1/3Co1/3O<sup>2</sup> ≡ *M*NMC is computed as the stoichiometrically weighted sums of each element using their atomic weights. From Table B.1, the sum molar mass is 96.46 gmol−<sup>1</sup> .

Table B.1: Atomic weights of the elements in LiNi1/3Mn1/3Co1/3O2.


By using this molar mass as well as the coefficients from Equation (B.5) and inserting these values into Equation (B.6) yields the specific capacity of NMC as

$$Q\_{\rm NMC,spec} = \frac{\mathcal{F}}{96.46} = 278 \,\text{mA} \,\text{hg}^{-1}.\tag{B.6}$$

## **B.2 Density from X-ray data**

The gravimetric density of a material *i* can be calculated by x-ray data. From [154], the density can be calculated by

$$\varrho\_{l} = \frac{NM\_{l}}{\mathcal{N}\_{A}V\_{l,\text{cell}}} \tag{B.7}$$

where *N* is the number of atoms per unit cell, *M<sup>i</sup>* is the molar mass, N*<sup>A</sup>* is the Avogadro constant and *Vi*,cell is the unit cell volume.

In the following, Equation (B.7) is applied to the positive electrode active material LiNi1/3Mn1/3Co1/3O2, which is used in this work. From Appendix B.1, the molar mass *M<sup>i</sup>* ≡ *M*LiNi1/3Mn1/3Co1/3O<sup>2</sup> ≡ *M*NMC is computed as the stoichiometrically weighted sums of each element using their atomic weights. From Table B.1, the sum molar mass is 96.46 gmol−<sup>1</sup> . Additionally, the values of the number of atoms per unit cell and the unit cell volume are taken from [141], where *<sup>N</sup>* <sup>=</sup> 3 and *<sup>V</sup>i*,cell <sup>=</sup> <sup>100</sup>.769Å<sup>3</sup> . Note that these values refer to a fully lithiated NMC metal oxide, i.e. γ = 1 in LiγNi1/3Mn1/3Co1/3O2. However, densities for other lithiation states are computed accordingly. Finally, the gravimetric density of NMC is

$$\varrho\_{\rm NMC} = \frac{3 \cdot 96.46}{6.022 \, 140.76 \cdot 10^{23} \cdot 100.769} = 4770 \, \text{kg} \, \text{m}^{-3}.\tag{B.8}$$

## **Schriftenreihe des Instituts für Angewandte Materialien**

### **ISSN 2192-9963**


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.



Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.








Band 98 Patricia Haremski **Diffusionseigenschaften von Nickel in einer Festoxid-Brennstoffzelle.**  ISBN 978-3-7315-1124-3 Band 99 Florian Wankmüller **Mehrskalige Charakterisierung der Hochtemperatur-Brennstoffzelle (SOFC).** ISBN 978-3-7315-1142-7 Band 100 Niklas Russner **Modellgestützte Analyse des Stackbetriebs von Festoxidzellen.**  ISBN 978-3-7315-1144-1 Band 101 Theo Fett, Karl G. Schell, Ethel C. Bucharsky, Gabriele Rizzi, Pascal Hettich, Susanne Wagner, Michael J. Hoffmann **Consequences of hydroxyl generation by the silica/water reaction - Part I: Diffusion and Swelling.**  ISBN 978-3-7315-1148-9 Band 102 Theo Fett, Karl G. Schell, Ethel C. Bucharsky, Gabriele Rizzi, Susanne Wagner, Michael J. Hoffmann **Consequences of hydroxyl generation by the silica/water reaction - Part II: Global and local Swelling Part III: Damage and Young's Modulus.**  ISBN 978-3-7315-1159-5 Band 103 Johannes Dornheim **Modellfreies Lernen optimaler zeitdiskreter Regelungsstrategien für Fertigungsprozesse mit endlichem Zeithorizont.**  ISBN 978-3-7315-1158-8 Band 104 Markus Muth **Grundlagenuntersuchungen an intrinsisch gefertigten lasttragenden FVK/Metall-Hybridträgern.**  ISBN 978-3-7315-1161-8 Band 105 Oleg Birkholz **Modeling transport properties and electrochemical performance of hierarchically structured lithium-ion battery cathodes using resistor networks and mathematical half-cell models.**  ISBN 978-3-7315-1172-4

## KARLSRUHER INSTITUT FÜR TECHNOLOGIE (KIT) SCHRIFTENREIHE DES INSTITUTS FÜR ANGEWANDTE MATERIALIEN

Hierarchically structured active materials in lithium-ion battery (LIB) electrodes with porous secondary particles are promising candidates for increasing the gravimetric energy density and rate performance of the cell. However, there are still aspects of this technology that are not fully understood. To gain a deeper understanding of how cathode structure and morphology influence cell performance, the goal of this work is to develop efficient tools to calculate effective transport properties for granular cathode structures that can be imported into cell models to evaluate the electrochemical cell performance of LIBs. To efficiently determine effective transport properties, the resistor network method is extended with respect to transport through the solid and pore phases of granular media and verified using finite element methods. A mathematical model for half-cells with hierarchically structured cathodes is proposed for the determination of electrochemical cell performance and qualitatively validated by experiments from the literature.

ISSN 2192-9963 ISBN 978-3-7315-1172-4 **O. BIRKHOLZ**

Modeling hierarchically structured lithium-ion battery cathodes