Modeling the Influence of Particle Shape on Mechanical Compression and Effective Transport Properties in Granular Lithium‐Ion Battery Electrodes

Lithium-ion batteries (LIBs) are presently one of the most favorable energy storage technologies. Due to high energy densities, good performance, and durability, the applications range from the automotive sector to consumer electronics. Especially for hybrid and electric vehicles, new active materials (AMs) are needed for high-capacity and high-power batteries. Regarding the performance of LIBs, the effective transport properties are of great importance. The microstructure of LIBs is typically granular. The positive electrode, for instance, consists of AM, binder, and conductive agent. About 90% of the total mass of the electrode is the AM, which is lithiated and delithiated during charge and discharge cycles. The remaining 10% are shared between binder and conductive agent, which often is carbon black (CB). Liquid electrolyte surrounds the porous solid skeleton, consisting of AM, CB, and binder. The interplay between all electrode components is quite complex. Themicrostructure as a result of the processing of a battery strongly influences the effective transport properties and, thus, the performance of a cell. Significant changes in the electrode microstructure may occur during manufacturing as a result of mechanical compaction, such as the calendering step. It has been proved that calendering plays a fundamental role, as the process significantly increases the number of particle-to-particle and particle-to-current collector contacts, which improves the electrical and thermal conductivities. Furthermore, the thickness is reduced while the volume of the AM is constant, which results in an increase in the volumetric energy density. At the same time, however, the porosity of the structure shrinks, which results in an increase in the electrode tortuosity and, thus, leads to a lower ionic conductivity of the electrode. To efficiently enhance the cell performance, it is important to find the optimal level of compaction. For this reason, understanding the underlying physics of granular media is of strong interest. As experimental measurements of 3D granular structures are complicated and difficult to analyze, computational techniques to test and to predict the behavior of granular materials have been frequently applied. In particular, the discrete element method (DEM) is a favorable modeling technique, as it can provide the relation between the macroscopic and the microscopic behavior of granular (particulate) assemblies. Furthermore, it allows to model nonspherical particle shapes, which are more realistic, concerning the natural grain shape of granular media. It has been widely recognized that the particle shape strongly influences the dynamics of granular systems, which, in turn, have an impact on physical properties. V. Becker, O. Birkholz, Dr. Y. Gan, Prof. M. Kamlah Institute for Applied Materials Karlsruhe Institute of Technology 76344 Eggenstein-Leopoldshafen, Germany E-mail: verena.becker@kit.edu

DOI: 10.1002/ente.202000886 The calendering step during manufacturing of lithium-ion batteries is an essential process in the production, as it significantly influences the microstructure of electrodes and, therefore, the performance of the battery. Within this context, this article investigates the influence of particle shapes on the micromechanical responses during calendering and, in turn, their impact on the effective transport properties of battery electrodes. The electrodes are modeled using discrete elements. For this reason, a novel algorithm for the generation of random stressfree particle assemblies consisting of superellipsoids is presented. The effective conductivities of solid and pore phase are calculated with a resistor network approach. In this context, a new analytical formula for calculating the individual resistance between two mechanically deformed ellipsoidal particles is presented. Furthermore, a geometrical approach is chosen for the pore phase for calculating individual resistances of pore throats in superellipsoidal particle assemblies. With the theoretical fundamentals, the effective conductivities of solid and pore phases of uniaxially compressed ellipsoidal particle assemblies are investigated. The mechanical response and its influence on the evolution of the effective conductivities are discussed. The deeper insight into the interplay between the calendering process and electrode microstructure can be a helpful information regarding a specific electrode design.
Several publications investigated the relationship between microstructure and effective transport properties in granular media making use of DEM. Ott et al. present a method based on discrete elements to model solid oxide fuel cells (SOFCs) and LIBs. [8] In their approach, the system consists of a binary mixture representing the AM and the CB phase. A fit formula for the estimation of the resistance between two mechanically deformed spherical particles is presented. Particle assemblies are densified dependent on the different manufacturing processes either by mechanical compression (LIB) or by geometrical adaption of the particle size to represent the sintering in SOFC. Furthermore, the volume change during operation in LIBs caused by intercalation and deintercalation of lithium into and out of the AM is modeled as a geometric increase or decrease in the particle radius. The authors present a resistor network (RN) approach for the estimation of the effective transport properties, which can be applied to uncompressed or compressed, and unlithiated or lithiated particle assemblies. For LIBs, they show that the conductivity generally increases during intercalation. A positive influence of calendering on the percolation threshold of particles in an assembly is found, as well. However, particles are assumed to be perfectly spherical, only.
Birkholz et al. present a similar approach to investigate the effective transport through the solid and the pore phase of LIBs. [12] In their work, they consider sintered structures for which the sintering process is simulated as in the previous study. [8] The pore phase is assumed to be shared by CB and electrolyte. Due to this assumption, not only a method to determine the effective solid phase conductivity is presented but also an approach to calculate the effective transport property of the pore phase. The pore phase is disassembled into pores and pore throats using a Voronoi tessellation. The individual resistance of a pore throat is then calculated based on a geometrical approach. The resistance formulas for the solid and the pore phase are used together with a RN method to determine the effective transport properties of the single phase. [12] Concerning the simulation of the mechanics of LIBs with the help of DEM, Giménez et al. investigate the influence of the calendering step on the mechanical properties of LIB electrodes. [13] They propose a DEM approach, which comprises the reproduction of real electrode structures as random spherical particle assemblies and their mechanical compression to simulate the calendering step during manufacturing. In this context, a contact model is presented, which combines the well-known elastic Hertzian model and a bond model to capture the viscous behavior of the binder. In the previous study, [14] they add an elastoplastic contact law to their modeling approach, which follows the theory of Thornton and Ning, [15] to better capture the mechanical properties of a lithium-nickel-manganese-cobalt oxide (NMC) AM. The calibration and validation are shown by comparison with nanoindentation tests on single NMC particle and nanoindentation measurements on whole electrode structures. In the previous study, [16] Giménez et al. consider the effective transport behavior in LIBs. They establish key correlations between the microstructure of an electrode with its specific properties and the effective ionic and electronic transport properties of an electrode by comparing simulations and real experiments. For the electronic conductivity, the parameters in the key relation are the solid fraction of the representative volume element (RVE), the coordination number meaning the number of contacts per AM particle, the fabric tensor component based on the interparticle contacts in the direction of compression, and the ratio of the contact radius to the particle radius and the conductivity of a single NMC particle. In case of the ionic conductivity, the key relation includes the specific free surface area, the ratio of the coordination number in direction of compression to the direction perpendicular to it, and the bulk conductivity of the electrolyte. In addition, both formulas include proportionality factors to fit them to real experimental results. Stershic et al. follow the idea to relate the electronic conductivity to structural quantities, too. [17] In their paper, they propose a fabric tensor formalism based on the interparticle contacts to describe the structure of a LIB and its evolution in case of mechanical loading. The fabric tensor analysis is applied to experimental data sets obtained by X-ray tomography and to DEM simulations of electrode microstructures. For the latter, the open source DEM code LIGGGHTS is used to create and to compress particle assemblies consisting of spherical particles with a size distribution. While the authors do show that fabric tensors capture the evolution of the interparticle contact distribution in real electrodes well and may, for this reason, be a good measure for the electronic transport within an electrode, the findings from DEM do not represent the experimental trend very well. Therefore, Stershic et al. draw the conclusion that a spherical representation of particles in DEM simulation is not sufficient to describe the mechanics of real electrode structures. [17] While the authors in the aforementioned publications make use of microscopic quantities such as the contact radius between particles in DEM or even on macroscopic properties of electrodes such as the fabric tensor calculated based on the particle contacts, Kespe et al. apply a spatially resolved 3D electrochemical continuum model of a LIB half-cell. In contrast to the previously mentioned publications, their method does not allow to simulate the influence of the compression. [18] For the calculation of the influence of the electrode microstructure on the cell performance, the approach comes along with an enormous computational effort. For this reason, small parts of LIB electrodes with a maximum of %100 particles may be considered, only, which may not be perfectly representative for the behavior of a whole electrode. This does not apply to the approaches by Ott et al., Birkholz et al., Giménez et al., or Stershic et al., where more general quantities with less detailed spatial information are used. In these cases, a large amount of several thousand particles can easily be analyzed.
Even though the previously mentioned works helped to give an insight into the calendering step and the calculation of effective transport properties, they almost exclusively concentrate on spherical particle shapes. In real electrode structures, however, particle shapes may strongly deviate from a sphere. [19] The transfer of findings from spherical to nonspherical particle systems is subject to doubts. Thus, the influence of particle shapes on the micromechanical responses during calendering and the impact on effective transport properties of battery electrodes are investigated in this work. We present a method to generate superellipsoidal particle assemblies at first. The particle assemblies represent RVEs inside of a real electrode. The calendering of an electrode is simulated as mechanical compression of the assemblies. Extending the approach given in the previous study, [12] a new method for calculating the effective conductivity of the solid particle phase and the pore space in (super-) ellipsoidal particle assemblies is shown, as well.
This article begins with the presentation of the methodological background including the particle shape description and all relevant contact relations. A new algorithm to create initial structures of superellipsoidal particle assemblies is derived by extending the random close packing (RCP) algorithm for spheres. [20] Next, the RN method for the calculation of the effective conductivity of solid and pore phases of cathode structures is reviewed. The formula used for the calculation of the resistance between two sintered spheres in the previous study [12] is taken and adapted to calculate the resistance between mechanically compressed ellipsoidal particles. A formula for the determination of the resistance through a pore throat in a superellipsoidal assembly is derived, generalizing the geometric approach in the previous study. [12] These formulas are verified on the basis of finite-element (FE) simulations. In addition, the validity is checked with the help of theories presented in the literature. Finally, random assemblies of spheroidal particles with different aspect ratios are created and compressed, to discuss the mechanical response and the resulting effective conductivity in view of the influence of particle shape on the physical properties and the impact on the performance of a cell.

Methodology for Calculation of Mechanical
Response and Transport Behavior

Definition of a Superellipsoid
Usually, spherical particles are used to simulate the behavior of granular media by DEM due to simple contact laws and, as well, a simple geometry-based contact point definition. However, in reality, individual grains are not perfectly spherical. A more general particle shape description, which not only covers spheres, but also ellipsoids, cylinders, or boxes, is the superquadric shape. By increasing the shape parameter from one in the case of spherical particles, namely, the radius, to five for superquadrics, which are three half-axes a, b, c and two blockiness parameter n, m, a lot more shapes can be easily defined; see Figure 1.
The surface of superquadric particles is described by the implicit equation of Barr, which is in local coordinates [21] The variables a, b, and c are the half-axes along the local main coordinates x, y, and z of the superquadric, respectively. The blockiness parameters n and m describe the sharpness of edges of the shape. If n ¼ m ¼ 2.0, an ellipsoid is described. If further a ¼ b ¼ c, the special case of a sphere follows. Convex and concave particle shapes may be generated. Sometimes, the use of the parametric representation of a superquadric is beneficial, where a specific point x is defined as For the trigonometric arguments in Equation (2) applies To avoid the presence of multiple contacts, the particle shapes in this work are restricted to smooth and convex forms, so-called superellipsoidal shapes, which means n, m ∈ ½2.0, ∞Þ. If the blockiness parameters in the exponents are n ≫ 2.0 and m ¼ 2.0, a cylindrical particle is considered. Furthermore, for n ≫ 2.0 and m ≫ 2.0, box-like shapes result. The edginess of the particle shape increases with increasing n and m beyond 2.0.

DEM for Mechanical Behavior of Granular Materials
In the so-called soft-sphere DEM, particles are individually identified based on their shape, the mass, and the physical and mechanical properties. [22] A contact between two particles is characterized by a virtual overlap with respect to the undeformed shape, which represents the local deformations in the contact zone. Detection of the contact interactions between particles and tracking the motion of each individual particle determines the behavior of the whole assembly. An interaction depends on the present contact model, which depends on the material properties of a particle. The particle motion is described by Newton's second law of motion, as well as the Euler differential equations. In this work, DEM is used to perform uniaxial compression tests (UCTs), which represent the calendering process during manufacturing of LIB electrodes. In a UCT, one dimension is reduced while keeping the dimensions perpendicular to the direction of reduction constant. A UCT is simulated for a www.advancedsciencenews.com www.entechnol.de particle assembly in a box domain, which characterizes an RVE of a battery electrode, by applying a prescribed deformation to one surface of the box in a series of load steps, which are separated by time increments Δt. Artificial damping is used to gradually remove all kinetic energy from the system, as quasi-static conditions are of interest. [23,24] The load is applied first for each increment, whereas for the relaxation process, the equations of motion are solved numerically by an explicit time stepping scheme to find an equilibrium state. [25] 2.

Mechanical Contact Model
The contact force in normal direction for elastic material behavior is well described by the widely known Hertzian model. [26,27] The Hertzian theory applies only under the condition that two contacting elastic solids have continuous surfaces, and the deformation induced is significantly smaller than the dimensions of each solid. Thus, the problem remains of localized nature, and the deformed regions in the contact zone can be approximated as quadratic curved surfaces, regardless of the overall shape of the solids. In the general case, when loaded in normal direction, two smooth convex solids will deform under compression and touch each other over an elliptical region, having the two semiaxes a c and b c ; see Figure 2. Furthermore, we define the equivalent radius c c ¼ ffiffiffiffiffiffiffiffi ffi a c b c p as the contact radius of a circle of the same size as the area of the elliptical contact region. The contact properties are mainly influenced by the local geometries in the vicinity of the contact point. These are well described by the principal The Hertzian normal force follows as The variable f 2 is one of two correction factors f 1 and f 2 to capture the eccentricity of the contact area. The correction factors are equal to unity for a circular contact. Although their expressions are quite complex, as described in the previous study, [27] approximate equations are often sufficient in terms of precision for practical purposes. Thus, the correction factor functions can be approximated according to the previous study [32] as The variable δ n describes the size of the virtual normal overlap, which depends on the contact point definition and which is described in Section 2.2.2.
In tangential direction, we make use of the simplified approach by Di Renzo and Di Maio. [33,34] The tangential force is taken as minimum between shear force and friction with the direction of the tangential displacement t. In this equation, F t;s is the tangential shear force, whose incremental change in each step is defined as with the effective shear modulus the correction factor to account for the ellipticity of the contact area and the incremental tangential overlap www.advancedsciencenews.com www.entechnol.de Finally, the total contact force, acting on a particle, can be calculated according to

Mathematical Contact Point Definition
The contact point definition of nonspherical shapes is, in contrast to spherical particle shapes, where the contact point is geometrically defined along the connection between the two sphere's centers, not straightforward. [11] As the superquadric equation defines the particle as a function in local frame, the contact detection algorithm chosen here follows the approach by Houlsby. [35] This algorithm makes use of the special definition of the superquadric equation, where the particle's surface is described by f ðxÞ ¼ 0 according to Equation (1). If f ðxÞ < 0, the point x lies inside of the particle. Furthermore, f ðxÞ > 0 describes points, which lie outside of the particle. Any surface must be strictly convex, which is automatically satisfied if the blockiness parameter of a superquadric is n ≥ 2.0, m ≥ 2.0, to assure to eliminate the possibility of multiple contact points between particles. Another necessary condition is that for any point x, not only the function itself must be evaluable but also the first and second partial derivatives. The surface of particles I and J can be expressed in a common global coordinate system X ¼ ðX , Y, ZÞ T in the form F I ðXÞ ¼ 0 and F J ðXÞ ¼ 0. Considering two superellipsoidal particles with strictly convex shapes and their shape functions F I ðXÞ and F J ðXÞ in global coordinates, the contact point is defined as point midway to both and at the same time closest to both, which leads to the optimization problem minimize½F I ðXÞ þ F J ðXÞ, subjected to F I ðXÞ ¼ F J ðXÞ This optimization problem can be solved by applying Lagrange multipliers, leading to a system of four equations with four unknowns, which can be solved by Newton's method. Details concerning the calculation of the contact point between two particles are given in the previous studies. [24,35] Having the contact point X c found, the normal overlap direction is defined as The contact point and the contact normal define the contact line; see Figure 3. The normal overlap δ n follows as distance between the two intersection points X I c and X J c of the contact line with the particle's surfaces. To find the intersection points, the following nonlinear algebraic equations (19) have to be solved with respect to the two scalars ξ I and ξ J . These equations can easily be solved in local frame by Newton's iteration, which is described in the previous study. [24] Having the two intersection points found, the normal overlap follows as δ n ¼ jX I c À X J c j.

Equations of Motions for Superellipsoids
The translational motion of each particle obeys Newton's second law of motion, which is where m and X C are the mass and the position of the particle center, respectively. F c is the total force, acting on a particle according to Equation (16). Due to contact forces, particles will translationally move. At the same time, contact forces may cause torques on particles, which lead to reorientation. The rotational motion is described by the Euler equations. [36] Referred to the local body-fixed frame, they are defined as where ω I is the angular velocity, i I is the tensor of inertia, and t I is the total torque acting on particle I in local frame. The representation in local coordinates is beneficial in this case, because the tensor of inertia in local frame contains nonzero entries in the diagonal, only. The solution of the equations of motion is done with a time stepping scheme. The so-called Velocity-Verlet algorithm can be used to integrate the translational movements, whereas a direct multiplication method serves to solve the rotational motion of a particle. Details can be found in the previous studies. [37][38][39][40] The transformation between global and local coordinate system can be performed with a quaternion multiplication, which is beneficial compared to the use of rotation matrixes in terms of memory and computation time. An explanation for the transformation with the help of quaternions is given in the previous studies. [39,41]

RCP for Generation of Initial Overlap-Free Particle Assemblies
For the simulation of the micromechanical behavior of LIB electrodes, we want to investigate an RVE. In the easiest case, spherical particle assemblies represent a periodic part of the real electrode structure. For the generation of spherical initial structures, the so-called RCP algorithm, first proposed by Jodrey and Tory, is a widely known algorithm. [20] As the generation of the RVE is an essential step, we present a novel RCP, which is able to produce overlap-free, dense, and random superellipsoidal particle assemblies. Monosized assemblies or assemblies having a particle size distribution may be created with the adapted algorithm, as well as structures consisting of different particle types, e.g., cylinders and boxes.
The idea of the modified RCP is still to eliminate the overlap δ I,J between two contacting superellipsoidal particles I and J by iterative size reduction and separation. In the initial state, N particles of k different types with fixed shapes are created in a virtual box of size L 3 . Fixed shape means that the aspect ratios α k ¼ a k =b k and β k ¼ a k= c k and the blockiness parameters n k and m k are constant during the whole procedure. Each position x ¼ ½x, y, z T of an individual particle is chosen at an arbitrary position within the box 0 ≤ x < L. Furthermore, each orientation, described by the Euler angles ϕ ¼ ½α, β, φ T , is random and within the range 0 ≤ ϕ < 360 . The orientation follows the x-convention of Euler, where the first rotation α takes place around the z-axis, the second rotation β around the new x-axis, and the third rotation φ around the new z-axis. [11,30] The first crucial step is the definition of the so-called outer and inner shape parameters s out ¼ ½a out , b out , c out T and s in ¼ ½a in , b in , c in T , as schematically shown in Figure 4. These parameters are defined as follows: the outer shape parameters characterize the overlapping state, whereas the inner shape parameters represent the overlap-free situation, where particles are just in touch in a single point but do not overlap. Decisive for the inner shape parameters is the worst overlap between two particles. As the outer shape parameters are always larger than or equal to the inner shape parameters, an outer packing factor η out can be calculated, defined as the ratio of the sum of all particles having size s out divided by the box volume V box ¼ L 3 . Calculating this ratio with the size of the inner shape parameters s in provides the so-called inner packing factor η in During the algorithm, the worst overlap has to be identified and to be removed by iterative size reduction and separation in each iteration step. Finally, when η out ¼ η in , which is the case for s out ¼ s in for all particles, the assembly is free of overlaps, and the configuration free of contact forces is reached, which can serve as a start configuration for DEM simulations.
To get the very first size of the outer shape parameters at the beginning of the RCP process, the half-axis a 1 of particles of type 1 is chosen to be the reference axis a R for all other half-axes of all particle types defining the ratios and for particles of type k 6 ¼ 1. Initially, the outer packing factor is chosen to be η out ¼ 1.0, which means that the sum of the volume of all particles is the same as the box volume. This leads to the condition where V k i is the volume of particle i of type k, and N k is the number of particles of type k. B is the so-called beta function. [42] Reformulation with the use of Equation (24) and (25) leads to Finally, solving Equation (27) gives the initial size of the leading axis a R by which all other outer shape parameters can be calculated according to Equation (24) and (25). Initially, many overlaps are in the system due to the initial packing factor of η out ¼ 1.0. The worst one has to be identified to get the values of the inner shape parameters. If the particles have different sizes, the overlap has to be related to the mean www.advancedsciencenews.com www.entechnol.de particle radius of both particles to calculate how much percentage the overlap occupies relative to the particle size. Thus, the absolute largest overlap in the whole system does not have to be the proportionally worst one. Keeping this in mind, the overlap of two particles in contact with each other is not just simply taken as the largest one. Instead, the overlap is weighted with the mean bounding sphere radius r I,J H of the two contact partners. This gives the related weighted overlapδ I,J where r I,J H is the mean bounding sphere radius, meaning the mean of the minimum radii r I H and r J H of the surrounding spheres that do touch the particles only in single points, but do not intersect with them. The calculation of the bounding sphere of a superellipsoid is given in the previous study. [24] In each iteration step, all contacts are searched, and the weighted overlaps are calculated. These are stored in a list among their wickedness starting with the worst one. For all contacts and belonging overlaps, is applies The inner shape parameters are calculated as a function of the worst overlap by searching for the size of the shape parameters for which the contact point X c of the worst contact lies exactly on the surface of both contacting particles, which is the case when f ðX c Þ ¼ 0. This leads to the condition for the inner shape parameters s in Introduction of the factor ρ gives Next, the two worst overlapping particles are moved away from each other dependent on their amount δ k , k ¼ I, J of the total overlap δ I,J n ¼ δ I þ δ J along the normal direction of the contact n I,J ; see Equation (18) and Figure 3. Subsequently, the outer shape parameters of all particles are contracted according to where i is the number of iterations, and Δη i is the difference between outer and inner packing factors calculated by Equation (22) and (23) based on the current values of s out and s in . ⌊•⌋ represents the Gaussian brackets. Applying the Gaussian brackets to a real number x, returns the largest integer, which is smaller or equal to x. The contraction of the particle's outer shape parameters gives the new values for the outer shape parameters in the next iteration step. The parameter c r in Equation (33) is the so-called contraction rate and an input parameter, which controls how fast the outer radius is contracted during each iteration step. As reported earlier, for example the previous study, [37] the resulting packing factor can be roughly controlled by c r . To check if this is still valid for the modified superellipsoidal RCP, spherical, ellipsoidal (spheroids), cylindrical, and box-like assemblies with different aspect ratios α r ¼ a=b ¼ β r ¼ a=c are created with four different values of the contraction rate c r . The exact input parameters are given in Table 1. Each assembly consists of 500 monosized particles, and three random cases are generated for each geometry data set.
As a first result, it is shown in Figure 5a,c that the trend of an increasing final packing factor with decreasing contraction rate is still valid for ellipsoidal, cylindrical, and box-like assemblies. Nevertheless, due to the more complex shapes, nonspherical structures reach lower packing factors for the same values of the contraction rate compared with spheres, which is shown in Figure 5b. Here, the final packing factor for the contraction rates c r ¼ ½0.01, 0.1, 1.0 is plotted over the aspect ratio α r ¼ β r for spheroidal shapes. It is evident that for an aspect ratio of 1.0, which represents a sphere, the final packing factor reaches the highest values. For box-like and cylindrical shapes, this trend is the same. Cylindrical and box-like structures with the same length of all half-axis a ¼ b ¼ c reach the highest packing factors, whereas the packing factor with an increase in the aspect ratio of the rotational symmetric structures decreases. In addition, the maximum packing factor reached for spheroids with a ¼ b ¼ c at a fixed contraction rate is higher than for cylinders with a ¼ b ¼ c. The packing factor for cylinders with a ¼ b ¼ c is, in turn, higher than for box-like particles with a ¼ b ¼ c. These observations occur, because the orientation of each individual particle is chosen randomly and cannot change during the algorithm. Furthermore, the higher the deviation from the spherical shape is, the more complex the shape and the mutual interference of particles. Due to the fixed orientation, more iterations and a slower contraction are necessary with increasing aspect ratio and higher edginess of the particle shape, to reach as high packing factors as in assemblies consisting of spherical particles.
On the other hand, we want to assure random orientation of particles to have no texture in the initial structure prior to calendering. The adapted RCP that fulfills this condition is shown in Figure 6, where the probability of the shortest half-axis for an oblate particle assembly (left) and the probability of the largest half-axis for a prolate particle assembly (right) are plotted based on the orientation of the considered half-axis for all particles in the assembly. This can be done using a fabric tensor analysis. Fabric tensors are a numerical approximation of a real directional data set in a multidimensional space. [43] There are three kinds of fabric tensors, where the first is the weighted moment tensor (called moment tensor) N, the second a least-square-based tensor (called fabric tensor) F, and the third an orthogonal decomposition of the second (called deviatoric tensor) D. Fabric tensors, which represent the directional distribution of unit vectorsñ, are symmetric referred to the origin f ðñÞ ¼ f ðÀñÞ. They sum up to one on a unit surface ∫ Γ f ðñÞdΓ ¼ 1, where Γ is a unit sphere in 3D. Based on the true directional distribution of the directional data f ∨ ðñÞ, a kth order approximation f ðkÞ ðñÞ can build up, where the entries of the fabric tensor serve as coefficients for a polynomial expansion. [17] Making use of spherical coordinates in 3D, this approximation function can be illustrated. To analyze the orientational distribution, we make use of the moment tensor of the second rank, which can be calculated as where k is the total number of unit vectors. As a unit sphere represents the state of perfectly random distribution, the almost perfectly spherical shapes in Figure 6 show that the individual orientations of particles in the assemblies are distributed randomly.

Resistor Network Approach
In this work, we make use of the RN method to calculate the effective conductivity, for both the solid and the pore phase of a cathode structure inside a LIB. The general idea behind RN, as used here, is the so-called node potential method. [44] This method is explained in-depth in the previous study [12] and, thus, will be briefly summarized here as background information for the extension to ellipsoidal and superellipsoidal particle structures, only. The mathematical formulation for the node potential method can be explained with the help of the exemplary network in Figure 7. In this graph, any nodes I and J are labeled with N I and N J , whereas the current between them is I I,J . Considering the boundary nodes of the network, N 0 and N 1 , between which all other nodes lie, an unknown effective current I eff results due to the potential drop between these boundary nodes. At each single node in the network, Kirchhoff 's current law represents the conservation of charge. The variable n neigh is the number of neighboring nodes of node N I . Furthermore,    Ohm's law is valid, where R I,J is the resistance of the resistor between two nodes N I and N J , and U I,J is the voltage drop between them, which can be expressed by the potential difference φ I À φ J . Combining both formulas results in which gives one equation for each node. Formulation of this equation for each node results in a system of linear algebraic equations (SLAE), which has to be solved to get the values of the unknown node potentials φ and the unknown current I eff . To this end, an arbitrary potential difference φ 0 À φ 1 ¼ U 0,1 between the two boundaries is chosen. Having the SLAE solved and, thus, knowing all potentials and the effective current, the effective resistance of the system follows as Finally, the dimensions of the considered domain have to be considered to get the effective conductivity. In this work, rectangular domains are considered with the cross-sectional area A dom and the domain length L dom . Thus, the effective conductivity follows as For details concerning the solution of the SLAE, we refer to Birkholz et al. [12]

Effective Transport Property of Solid Phase
We approximate the granular structures of LIB cathodes by discrete elements, where for each element, overall properties are known, only. Concerning the geometry, each discrete element is clearly defined by the position of the particle center, the orientation, and the particle shape. Now, regarding the effective conductivity through particles of a certain phase, we present, in this section, a formula for the resistance due to the bottleneck effect along the contact area between two particles pressed against each other.
First of all, to make use of the RN method for the solid phase, conducting pathways have to be searched in the particle assembly. Such percolated paths are symbolized in Figure 8. To find all percolated paths, all clustering particles, which connect two opposite boundaries of the structure with each other, have to be identified. Subsequently, these paths are converted into a RN, where we take the center of each particle I as a node N I , to which we assign the potential φ I , whereas a resistor R I,J is attributed to the edge between any nodes I and J in mechanical contact. The resistance between two individual particles can be determined dependent on the mechanical deformation of the contact zone and the arising mechanical contact area between them. At particles in touch with the boundary, additional nodes are modeled to impose the boundary conditions there. This procedure results in a RN, which is needed to calculate the effective solid conductivity according to the scheme presented in Section 2.4.1.
The calculation of a single resistance between two individual particles is based on the mechanical contact area, which arises with increasing compression. The formation of the contact area, in turn, depends on the material properties, e.g., elastic or possibly plastic material behavior. [45] In the elastic case, for general shapes, the contact area will be elliptical, as explained in Section 2.2.1. In this context, the size of the contact area calculated by DEM is defined by Equation (4)- (7).
For the geometrical resistance between two sintered spheres, the authors of the previous study [12] make use of the analytical resistance formula presented in the previous study [46] Here, κ I and κ J are the bulk conductivities of the conducting materials, and r c is the contact radius of the two spheres. This solution can be used in conjunction with the equivalent radius of the mechanical elliptical contact area according to Section 2.2.1, which gives r c ¼ c c , to describe the elliptical situation. Thus, the solid resistance formula between two mechanically compressed ellipsoids follows as

Verification of Resistor Network Method for Solid Phase
To verify the solid resistance formula, taking advantage of the mathematical equivalence between stationary charge and thermal transport, [12] FE analyses (FEAs) are taken as a reference solution. The FEA provides the spatially resolved solution of the stationary boundary value problem. If the mesh refinement www.advancedsciencenews.com www.entechnol.de is sufficient, the solution can be considered to be exact. The FEAs are calculated with the software Abaqus. [47] In the first step of the FEA, two elastic ellipsoidal bodies, as shown in Figure 9, are compressed displacement driven in x-direction. The contacting particles are either spheroids or triaxial ellipsoids with aspect ratios in the range of α r ¼ 0.5..2.0 and β r ¼ 0.5..2.0, where α r ¼ a=b and β r ¼ a=c. The half-axis a points in the polar direction, whereas b, c are the half-axes in the equatorial plane. For prolate spheroids, it applies a > b ¼ c, and thus, α r > 1 and β r ¼ α r , whereas for oblate spheroids, the aspect ratio is α r < 1 and β r ¼ α r , and thus, a < b ¼ c. A sphere is defined through α r ¼ β r ¼ 1 and a ¼ b ¼ c. Triaxial ellipsoids are characterized by α r 6 ¼ 1.0 6 ¼ β r . The considered spheroids are volume consistent, and thus, for different aspect ratios, the half-axes sizes vary, only, whereas the mean particle radiusr ¼ ffiffiffiffiffiffiffi abc p is the same. For spheroids, four different contact configurations are considered; see Figure 10. Three additional contact situations with triaxial ellipsoids of different particle sizes are investigated, as well. The contact configurations are chosen, such that symmetry can be exploited. The geometry data and the material properties are listed in Table 2.
In the second step, the deformed ellipsoidal geometries are imported, and a temperature gradient is imposed on the xy-symmetry planes of both ellipsoids, which is shown in Figure 11. The value of the temperature gradient is set to one for simplicity. Thus, a flux through the mechanical contact area emerges between both particles. To determine the resulting resistance R I,J between them, the resulting total heat flux Q FEA at the xy-symmetry plane of one of the two particles is evaluated, to finally obtain the resistance as This value is then compared with the result according to the analytical formula in Equation (41), which is plotted in Figure 12.
In the left graph, the verification for the contact between two spheres, two oblates, touching as in case four of Figure 10, and two prolates, touching as in case five of Figure 10 is shown. As the calculated curves and points of all cases according to Table 2 lie close to each other, these three cases are shown exemplarily here, only. The conductivity in all cases is κ I ¼ κ J ¼ 1.0. In the right part of Figure 12, the comparison between FEA and analytical formula for the three triaxial cases is presented. In these cases, the ratio of the conductivities of both contact partners I and J varies between κ I ¼ κ J ¼ 1.0 (hollow markers) and κ I =κ J ¼ 1=10 000 (filled markers). The average mean error of all considered cases between FEA and analytical formulae is below 5%. The good agreement between FEA and analytical formulae approves the here chosen approach for the calculation of the effective conductivity between two solid, mechanically deformed ellipsoidal bodies.   From the computational point of view, a FEA for a FE model of two contacting eighths ellipsoids with around 300 000 nodes for the mechanical compression and the heat transfer analysis calculated with eight central processing units (CPUs) on the workstation Intel(R) Xeon(R) CPU E5-1620 v4 @ 3.50 GHz with 32 GB ram specifications needed about 23 GB ram and around 86 400 s, while the resistance calculated by Equation (41) is immediately there.

Effective Conductivity of Pore Phase
For the application of the RN method to the pore phase, we follow the approach developed by Birkholz et al. [12] and extend it to superellipsoidal particle shapes. The idea is to discretize the pore phase by a Voronoi tessellation. In general, the Voronoi tessellation assigns a spatial cell to an isolated point N I in a domain. Each cell contains all points, whose distance is less or equal to a neighboring cell node N J . Thus, a Voronoi cell has, in 3D, the shape of a convex polyhedron. For the pore phase, each particle center is chosen to be the center of a Voronoi cell, which is consequently build up around the related particle. In this way, the whole domain can be discretized into cells. For this discretization process, the open source library Voroþþ is used. [48] Figure 13 shows a Voronoi tessellation. The part of a Voronoi cell, which is not covered by the particle, is part of the pore space. Each Voronoi vertex is the center of a pore, whereas the corresponding edges represent pore throats. Hence, network nodes N I are set to the pore centers, where a potential φ I is assigned, whereas a resistance R I,J is related to an edge connecting to pore center N J . Furthermore, all boundary cells have to be identified to apply the boundary conditions there. This procedure results in a RN according to Section 2.4.1.
The crucial part is the calculation of the individual resistances of pore throats, where we use a geometrical approach similar to the previous study. [12] Each pore throat element is subdivided into m subpore elements along the considered edge; see Figure 13. Furthermore, each subpore element is divided into n sufficiently small increments ΔL, where for each segment, the resistance can be calculated according to Here, ρ is the resistivity, and A mean ¼ 0.5ðA n þ A nþ1 Þ is the mean cross-sectional area for the nth increment; see Figure 14.
The cross-sectional area cannot simply be determined based on the radius and the boundary nodes of the cross section, as it can be done for spherical particles. Instead, all intersection points    www.advancedsciencenews.com www.entechnol.de of one cross-sectional element with the particle's surface have to be determined iteratively. Having in mind that for any point on the surface of a superellipsoidal particle f ðxÞ ¼ 0 applies, the condition, which has to be solved, written in global frame, is In Equation (44), X k describes all points of the cross section, which intersect with the particle's surface. n k are normal vectors along the "particle-free" edge of the considered cross-sectional area and, as shown in the right part of Figure 14, which point from the edge to the particle center. Equation (44) can easily be solved in local frame by Newton' iteration with respect to ξ k . If enough increments k are chosen, the intersection line between particle and cross-sectional area can be determined precisely. When all intersection points are found, the respective cut sections can be calculated using the Boost library. [49] Having all cut sections found, the resulting resistance of one pore sub-element follows as series connection of the n incremental resistances The length ΔL of the n increments is successively reduced during the calculation process, until the subresistance value of the next increment does not differ from the previously calculated one by more than a certain threshold, which is chosen here to be 5%. Finally, the total pore resistance follows as parallel connection of all m subpore elements along one edge as

Verification of Resistor Network Method for Pore Phase
To verify the abovementioned formula for the pore throat resistance of superellipsoidal particle systems, we again make use of FEAs. Due to the absence of a fully analytical formula, the FEAs are assumed to be the exact solution. Three different types of particle assemblies are considered. First of all, ellipsoidal systems with 200 particles are created by the RCP presented in Section 2.3. Assemblies with different size distributions with the standard deviations of σ ¼ 0.0, 0.05, 0.1, 0.2 are considered. The geometry data is listed in Table 2.
In a second step, the overlap-free initial systems are artificially densified, by increasing the particles half-axis about 1.0, 2.5, 5.0% relative to their size, neglecting any mechanics. This leads to overlaps in the structure, as they would occur in mechanically compressed discrete element structures. For the verification, we fall back to an artificial densification due to the computational limitation of the FEA. Namely, the mechanical compression of 200 particles, spatially resolved with sufficient fine meshing using Abaqus, is computationally not possible. In addition, as particle deformation due to mechanical compression is limited to local regions, only, and virtual overlaps are merely of a few percent, the overlapped state can be assumed to be sufficiently accurate to approximate the compressed deformed state of a particle assembly. For the simulation of the conductivity with FEs, 200 particles are approximately the maximum number of particles for which a solution can be gained in terms of computation time and memory. On the other hand, going back to Section 2.4.3, when calculating the resistance for transport through two particles in contact, it turned out that such a simplification for that case is not admissible. Evaluating the resistance in Equation (41) based on a contact area obtained by growing the particles self-similarly gives no reasonable approximation for the resistance based on the Hertzian contact area.
For the FEA, the geometry data of the superellipsoidal structures generated with RCP and artificially densified are imported into a box-shaped domain, where all parts of particles outside of the domain boundaries are cut off. The complement of the particle phase is the pore phase of the structure in the box domain. As we are interested in transport through opposite faces, all nodes of the top surface are set to the same temperature. The temperature assigned to all nodes of the opposite surface, the bottom surface, is dropped by a temperature gradient of ΔT ¼ 1.0, for simplicity. The other four box faces are assumed to be insulated.
Finally, the resulting heat flux Q FEA has to be evaluated at one side of the box faces, where the temperature boundary condition is applied, to get the effective conductivity of the FEA. Knowing the domain length L dom , as well as the cross-sectional area A dom , the effective conductivity follows as For the calculation of the effective conductivity with the RN method, a potential drop of Δφ ¼ 1.0 is imposed on the opposing boundary faces, whereas no boundary condition is applied at the nodes of the four other box faces. According to Equation (39), the box dimensions have to be considered to calculate the effective conductivity. It yields (48) Note that the used bulk transport properties of the materials are set to a unit value, so that all results, finite element method (FEM) and RN, refer to unit bulk properties. Multiplying these normalized values with the corresponding real material bulk property gives the effective transport property of the considered granular material.
The pore phase conductivity of ellipsoidal particle assemblies according to Table 3 is presented in Figure 15. The black solid line represents a perfect match, whereas the blue dashed line demonstrates the mean deviation between FEM and RN. Each case represents a geometry data set according to Table 3 for which four different calculations have been performed. In general, for all 48 simulations, the RN method overestimates the FEM solution by about 0.37%. The absolute error between both methods is 2.22%, which leads to the conclusion that the FEM and RN are in very good agreement.
To verify the pore phase conductivity of superellipsoidal particle assemblies according to Table 4, the comparison of FEM and RN for cylindrical and box-like structures is presented in Figure 15. Again, the black solid line represents a perfect match, whereas the blue dashed line demonstrates the mean deviation between FEM and RN. Here, according to the 44 simulations, the RN method underestimates the FEM by about 0.42%. The absolute error between both methods is 1.81%, which is a good match again.
Due to the overall good agreement between FEM and analytical prediction, the RN method is successfully verified for the calculation of the effective pore phase conductivity in nonspherical superellipsoidal particle assemblies ( Figure 16).

Discussion of Derived Resistance Formulas for Solid and Pore Phases
For the pore space of an assembly of particles, a verification could be achieved by FEAs. However, other than a two-particle contact, FEA of mechanical compression of a whole particle assembly is not feasible. Thus, so far, the resistance for transport through two particles in mechanical contact could be verified, only. To check the validity of the effective conductivities calculated with the RN method, the results are compared with well-known theories in the literature. An upper and lower bound can be calculated according to the definition of Hashin and Shtrikman (HSh). [50] These well-known bounds are defined as κ eff hs1 ¼ 2κ solid þ κ pore À 2ðκ solid À κ pore Þϕ pore 2κ solid þ κ pore À ðκ solid À κ pore Þϕ pore κ solid , κ eff hs2 ¼ 2κ pore þ κ solid À 2ðκ pore À κ solid Þϕ solid 2κ pore þ κ solid À ðκ pore À κ solid Þϕ solid κ pore (49) Furthermore, the result of h ð3ϕ pore À 1Þκ pore ð3f1 À ϕ pore g À 1Þ þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3ϕ pore À 1Þκ pore þ ð3f1 À ϕ pore g À 1Þκ 2 the effective medium theory (EMT) divides the area between the two boundaries by HSh into an upper and a lower part. The upper part belongs to materials, which have an inner porosity such as the structure of a sponge, which applies to the pore phase of a LIB electrode. On the contrary, structures having an external porosity, such as the solid phase of battery cathodes, can be found in the lower part. [51] The effective conductivity for the solid phase, as well as the pore phase, of all particle assemblies according to Table 3 and 4, is plotted in Figure 17 together with the theoretical bounds. In the left graph, the comparison of the effective solid conductivity, provided by the RN method, is shown. In this case, the pore phase does not contribute to the transport, and thus, it is κ pore ¼ 0. As κ solid > κ pore , the lower bound is described by the second equation of Equation (49). The solid phase should be located, as explained before, between the lower bound of Hashin-Shtrikman and the bound calculated based on EMT, Table 3. Geometry data of ellipsoidal particle assemblies.  Figure 15. Comparison of FEA results and RN results for the effective conductivity of the pore phase for polydisperse ellipsoidal particle assemblies. Mean deviationε ¼ þ0.37%. Geometry data according to Table 3. Table 4. Geometry data of cylindrical and box-like superellipsoidal particle assemblies.  Figure 16. Comparison of FEA results and RN results for the effective conductivity of the pore phase for polydisperse cylindrical (filled points) and box-like (not filled points) particle assemblies. Mean deviation ε ¼ À0.42%. Geometry data according to Table 4.
www.advancedsciencenews.com www.entechnol.de which applies. Furthermore, the pore phase represents the case of a sponge-like structure, wherefore its effective conductivity values are expected to be found between the EMT and the upper HSh bound. Now, we assume that the solid phase does not contribute to the transport process through the pore phase, which means that κ solid ¼ 0. The interchange of the upper and lower bound in case of κ solid < κ pore must be considered. Thus, the upper bound is described by the second equation of Equation (49) for the pore phase, again. The comparison in the middle and right graphs of Figure 17 supports the expectation that the results of the RN method lie between EMT and upper HSh bound. Consequently, from this point of view, the RN method provides reasonable results. The solid phase conductivity decreases with increasing porosity, whereas the pore phase conductivity increases, which is in accordance with the theories. Finally, we conclude that for the considered particle assemblies, particle shapes, and aspect ratios, the RN method provides good results, for both the pore and the solid phase of granular structures.

DEM Simulation of Representative NMC-Cathode Volume Element
In this study, we investigate the influence of the particle shape on the mechanics and on the effective transport properties of battery electrode microstructures. Therefore, RVEs consisting of 1800 monosized ellipsoidal particles are created with the modified RCP presented in Section 2.3. The ellipsoids have either oblate or prolate shapes with the aspect ratios of α r ¼ β r ¼ 0.667, 0.769, 0.714, 0.833, 0.909, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, respectively. For each aspect ratio, five random cases are created. The RVEs have periodic boundaries and represent a small part inside the electrode. The initial packing density of all systems is around ϕ ini % 55.0 AE 0.03%, which is visualized in Figure 18.
Here, each point represents the mean value of five cases per aspect ratio with the standard deviation shown as error bars.
As mechanical DEM is very sensitive to the initial packing factor, it is important to have an almost constant initial packing factor to be able to compare the different assemblies and their properties with each other. First, we analyze the initial configurations in terms of their effective transport properties. In a second analysis, the UCTs are performed to simulate the calendering step during manufacturing of a battery electrode. Therefore, the particle assemblies are subjected to a uniform uniaxial strain field in z-direction, while keeping the dimensions perpendicular to the direction of reduction constant. The prescribed compression is applied to the upper box surface in z-direction in a series of load steps, which are separated by time increments Δt. Furthermore, a strain-rate tensor ε : ij is specified to control the deformation of the periodic cell. [52] Dependent on this strain-rate tensor and the position x i of a particle in the cell, the centers of all particles are moved, as if they would be points in a continuum by The total reduction is 20% of the box length in the compressed direction. For each converged step, the average stressσ ij in the system, which can be obtained bȳ is analyzed. [53] The variables, f n and f t are the magnitudes of the normal and tangential contact forces applied from particle J on I, δ c denotes the distance of the particle's centers, and n i , t i are the normal and tangential unit vectors, respectively. [52] The orientation change of the particles is characterized by a fabric tensor analysis to see if any texture forms in the assemblies as a result of compression. As indicated in Section 2.3, fabric tensors are a numerical approximation of a real directional data set in a multidimensional space. Here, we make use of the weighted moment tensor N, see Section 2.3 and the previous study, [17] to analyze the orientation of particles, and of the fabric tensor F, see the previous study, [17] to characterize the contact normal of the arising contacts with increasing compression. The fabric tensors F of the second and fourth ranks, respectively, can be calculated as Figure 17. Comparison of RN results with theoretical boundaries of effective conductivity for porous media. Effective conductivity of solid phase of ellipsoidal particle assemblies (left) and of pore phase of ellipsoidal particle assemblies and of superellipsoidal particle assemblies (right). Geometry data according to Table 3 and 4. Figure 18. Initial packing factor ϕ ini of particle assemblies with 1800 monosized spheroidal particles versus aspect ratio α r ¼ β r .
www.advancedsciencenews.com www.entechnol.de and may be visualized by a rose diagram or a polar plot that indicates the likelihood of distribution as a radius length for the given direction. In the last step of the analysis, the effective conductivities for the solid and the pore phase of each assembly are calculated for different compression states according to Section 2.4.1.
The particles are assumed to be elastic and to have the properties of lithium-NMC according to the previous study. [16] Thus, Poisson's ratio of ν ¼ 0.25, Young's modulus of E ¼ 142 000 MPa, and a mean radius of 5 μm are chosen.
At this point, it has to be remarked that a unit value of the bulk transport property for the material of the considered phases is assumed for all results shown here. These normalized values have to be multiplied by the corresponding bulk transport property of the material of interest to obtain the real effective transport properties of a granular system.
Starting with the analysis of the initial configurations, the effective conductivities of the solid and the pore phase are of interest. Due to the definition of the RCP, there are no overlaps in the system yet, and thus, no mechanical contact areas are present, which is synonymous with no transport through the solid phase. Nevertheless, transport through the pore phase is possible and can be analyzed with the formulas presented in Section 2.4.4. In Figure 19, the effective conductivity is plotted over the aspect ratios of the monosized ellipsoids. As the initial orientation of the particles is random, see Figure 6, the effective conductivities in the x-, y-, and z-directions are in the same range, and thus, the average effective conductivityκ pore eff ¼ 1=3ðκ pore eff ,x þ κ pore eff ,y þ κ pore eff ,z Þ is analyzed. Each point in Figure 19 represents the meanκ pore eff ¼ 1=5 P 5 n¼1κ pore eff ,n of the average effective conductivity of all five cases per aspect ratio with the standard deviation shown as error bars.
It is obvious that the effective conductivity for the different aspect ratios strongly differs. Taking the sphere solution α r ¼ β r ¼ 1.0 as a reference, there is a slight increase in the effective conductivity for prolate shapes, beforeκ pore eff linearly drops for aspect ratios larger than α r ¼ β r ¼ 1.4. The effective conductivity reaches its maximum value for oblate shapes for α r ¼ β r ¼ 0.615. In general, higher effective conductivities compared with the sphere solution are found for oblate spheroids and the considered aspect ratios. To understand these observations, the pore phase structure is analyzed. We recall that for the calculation of transport through the pore phase, the latter is disassembled into Voronoi cells. Each Voronoi cell contains several subpore elements, which belong to a certain pore throat. To get the resistances of the different pore throats, the mean cross-sectional areas A mean of each subpore element are of interest; see Section 2.4.4 and Figure 14. As a result of this, the pore structure can be characterized by the sizes of the mean crosssectional areas. In this context, it is important to keep in mind that a certain increment length ΔL belongs to each mean crosssectional area, and that the sizes of the increment lengths may not be the same. Thus, we sum up the sizes of all p mean crosssectional areas of all cases P p A mean;p and weight the result by the sum of the total increment length P p ΔL p . This gives the size of the overall average mean cross-sectional areaĀ mean ¼ P p A mean;p = P p ΔL p in the pore space of all considered initial cases and aspect ratios. Then, for each case and aspect ratio, the distribution of the mean cross-sectional areas is discretized into 100 segments and weighted by the size of the overall average mean cross-sectional area. The relative frequency for each segment cannot be simply gained by summing up the amount of cross-sectional areas per segment. As the increment length of each cross-sectional area may not be the same, the size of the increment length must be considered. Thus, the sizes of all belonging increment lengths of the cross-sectional areas in one segment are summed up and weighted by the mean increment length to get the relative frequency of one segment. This is visualized in Figure 20 and 21. Again, each curve represents the average of the five considered cases per aspect ratio.
Spheres have a clear maximum at a certain size of the mean cross-sectional area. For oblate ellipsoids, Figure 20, this maximum decreases and shifts toward almost infinitesimal small areas, so that for the case of α r ¼ β r ¼ 0.5, a continuously drop in the relative frequency with increasing mean cross-sectional Figure 19. Effective conductivityκ pore eff of pore phase of particle assemblies with 1800 monosized spheroidal particles versus aspect ratio α r ¼ β r . Figure 20. Relative frequency pðA mean =Ā mean Þ of weighted average mean cross-sectional areas A mean =Ā mean of pore phase of particle assemblies with 1800 monosized oblate ellipsoids.
www.advancedsciencenews.com www.entechnol.de area size is found. For prolate ellipsoids, Figure 21, two maxima evolve with increasing aspect ratio, where the first and higher maximum lies again at the position of almost infinitesimal small areas, whereas the second maximum forms at a certain segment size of the mean cross-sectional areas; see Figure 21. Comparing this observation with the results for the effective conductivities for the pore phase in Figure 19 leads to the conclusion that a pore structure with one mid-size peak, in the optimum case close to the one obtained here for an aspect ratio of α r ¼ β r ¼ 0.615, seems to be beneficial. On the contrary, the formation of two different maxima results in a strong deterioration of the effective conductivity.
In the next step, we analyze the compression tests of the prolate and oblate particle assemblies. The stress-strain relations for prolate (left) and oblate (right) ellipsoids are shown in Figure 22. Each curve represents the mean of the five considered cases per aspect ratio. It is evident that the stress response at the same level of compaction, for both prolate and oblate shapes, decreases with increasing aspect ratio. This can be explained with observations made for the theoretically densest packing of spheres and ellipsoids. The theoretical densest packing of spheres is the facecentered cubic (fcc) packing with a maximum packing fraction of ϕ fcc;sph ¼ 0.7405. [54] In contrast to spheres, ellipsoids have three additional degrees of freedom and, thus, an orientation in space. Dependent on the aspect ratios and the particle orientations, higher packing densities than ϕ fcc;sph are possible. Two pronounced maxima of ϕ ell ¼ 0.770732 can be observed for an aspect ratio of α r ¼ β r ¼ 1= ffiffi ffi 3 p (oblate ellipsoids) and α r ¼ β r ¼ ffiffi ffi 3 p (prolate ellipsoids), when the number of in-plane neighbors reaches six. This maximum packing density decreases for both types of spheroids towards ϕ fcc;sph with four in-plane neighbors, when the spheroidal shapes get closer to the shape of a sphere; see the previous study [54] and Figure 24. In addition, a constant maximum packing density has been found for ellipsoids with α r ≥ ffiffi ffi 3 p and α r ≤ 1= ffiffi ffi 3 p . [54] As the initial packing density of all considered cases is nearly the same, ϕ ini % 55.0 AE 0.03%, spheroids show a softer stress response due to their higher theoretical packing density. Furthermore, as the theoretical packing factor increases with increasing aspect ratio, the stress response gets softer, the larger the aspect ratio of a spheroid.
To check if any texture forms due to reorientation, we analyze the orientation of the spheroids with a moment tensor, see Section 2.3, for the initial state and at maximum compression. In Figure 23, the moment tensor of the shorter half-axis in case of oblate ellipsoids and the longer half-axis in case of prolate ellipsoids is plotted before (top) and at maximum (bottom) compression. Starting with a random initial orientation, which is represented by the spheres in the top line of Figure 19, it can be seen that for oblate ellipsoids, the z-direction gets more favorable, which means that the shorter half-axis tries to align along the z-axis due to compression in z-direction. For prolate shapes, the z-direction gets less favorable, which is synonymous with a rotation of the longer half-axis toward the xy-plane. From a mechanical point of view, this behavior is understandable. If we compare the resulting texture with the theoretically favorable texture according to Figure 24 and the previous study, [54] one can see that prolate ellipsoids move closer to their theoretically densest configuration due to reorientation, in which the longest halfaxis in one layer is aligned along one direction, whereas this direction is shifted by a certain angle in the next layer. Figure 21. Relative frequency pðA mean =Ā mean Þ of weighted average mean cross-sectional areas A mean =Ā mean of pore phase of particle assemblies with 1800 monosized prolate ellipsoids.   Nevertheless, the axis of rotational symmetry being here the longest axis is always found in the xy-plane, which is shown in Figure 24. This observation underlines the softer stress response of prolate shapes. In the densest configuration of oblate ellipsoids, the shortest half-axis of one layer is again aligned along a certain direction, whereas the favorable direction is rotated about a certain angle in the next layer. Again, the axis of rotational symmetry, i.e., the shortest axis, in both layers is found in the xy-plane. If we compare this favorable configuration with the observed texture, we see that the reorientation of oblate spheroids, where the shortest half-axis aligns along one axis, the z-axis in this case, is not in accordance with the ideal densest packing structure, where the favorable direction of the shortest axis changes from layer to layer. This may explain why the stress response of oblate ellipsoids, in particular for the aspect ratios α r ¼ β r ¼ 0.769, 0.714, 0.667, is not as soft as for their prolate equivalents α r ¼ β r ¼ 1.3, 1.4, 1.5.
In the next step, it is of interest weather and how the arising texture during compression affects the effective transport properties of the assemblies. Therefore, the effective conductivities of the pore and the solid phase are plotted over the compression in Figure 25 for prolate and in Figure 26 for oblate spheroids.
First of all, as it is expected, the effective conductivity of the pore phase reduces during compression, as the amount of pore phase decreases. Although a texture is formed in the assemblies, we do not observe big differences for the effective conductivity of the pore phase in the x-, y-, and z-directions, and thus, the average mean effective conductivityκ pore eff is considered here. However, the effective pore phase conductivity decreases faster for an increasing aspect ratio with further compression. This effect is more pronounced for prolate shapes. The analysis of the pore phase structure in the way it has been done for Figure 20 and 21 is visualized in Figure 27. To understand the larger drop of the effective conductivity of the pore phase with increasing aspect ratio, we compare the results in Figure 27 with the initial states shown in Figure 20 and 21. Starting with the sphere situation, no big differences between the initial state and at maximum compression can be found. This fact applies to the aspect ratios α r ¼ 0.769 and α r ¼ 1.3, as well, for which the decrease in the effective conductivity of the pore phase during compression happens to the same extend than for spheres; see left graphs of Figure 25 and 26. For the aspect ratio α r ¼ 1.5, the shape of the pore structure curve gets closer to the initial situation of α r ¼ 2.0, which with its two peaks was found to be disadvantageous. Big changes occur for α r ¼ 0.667, as well, where the pore space structure changes toward a shape similar to the initial state of α r ¼ 0.5. This curve progression was found to be not as bad as the shape with two maxima, but worse than one maximum at a certain mean cross-sectional area size larger than zero. Thus, the smaller drop for α r ¼ 0.667 compared with α r ¼ 1.5 is in accordance with the observed changes of the effective conductivities for the different pore structures in Figure 19.
Concerning the effective conductivity of the solid phase, it is visible that a certain threshold has to be passed before enough contacts in the assemblies have been formed and the effective  . Theoretically densest laminate crystal packing ϕ max of ellipsoids versus aspect ratio α r ¼ β r . The fcc packing of spheres ϕ fcc;sph ¼ 0.7405, when α r ¼ β r ¼ 1.0. Maximum density ϕ ell ¼ 0.770732 of ellipsoids for ffiffi ffi 3 p (oblate ellipsoids) and α r ¼ β r ≥ ffiffi ffi 3 p (prolate ellipsoids). Graph based on the previous study. [54]  www.advancedsciencenews.com www.entechnol.de conductivity increases. An interesting observation regarding the effective conductivity of the solid phase is that the latter increases faster with increasing aspect ratio. This effect is found for both prolate and oblate ellipsoids. While, for example, the stress at maximum compression for α r ¼ 1.3 is only one-fifth of the stress reached by spheres, the effective conductivity of the solid phase of α r ¼ 1.3 is about three-fourth of the value of spheres. To understand this trend, the mechanical contact areas and the related contact forces, as well as the coordination number, which is defined as the average number of contacts of a particle, are analyzed. The result of the aspect ratios α r ¼ 0.769, 0.667, 1.0, 1.3 is shown in Figure 28, where the ratio of the total sum of the mechanical contact areas of the ellipsoidal systems to the sphere assemblies is plotted in the left graph, the ratio of the total sum of the normal contact forces of the ellipsoidal systems to the sphere assemblies is shown on the right-hand side, and the coordination number is illustrated in the bottom graph. It can be seen that the increase in the ratio of the contact areas with increasing compression is steeper and larger than the increase in the ratio of the contact forces. This means that there can be larger contact areas in ellipsoidal systems compared with sphere systems for the same contact force. One could ask how this is possible. Let us take a step back and look at the texture that forms under compression. The shorter half-axis of oblate ellipsoids aligns in the z-direction, whereas the longer half-axis of prolate spheroids rotates into the xy-plane. As the assemblies are compressed in the z-direction, we expect more contacts along this direction. To this end, we analyze the contact normal with a fabric tensor analysis. Figure 29 exemplarily shows the fabric tensor of rank four F ð4Þ for α r ¼ 0.667, 0.769, 1.3, 1.5 at maximum compression. From this graph, it gets clear that most of the contact normal are aligned along the z-axis, because the z-direction is more favorable.
Together with the observations made concerning the texture of the systems, we can conclude that the oblate ellipsoids touch preferentially at poles, whereas for prolate spheroids, contact at the equator is favorable. The relation between contact force and related contact area for idealized configurations of two particles is visualized in Figure 30. In this graph, the ratio of the spheroid relation to the sphere relation is shown for oblate ellipsoids touching at their poles on the left and for prolate ellipsoids touching at their equators aligned or rotated about 90 on the right-hand side. The ratio is evaluated for different overlaps represented by the different dots, whereas the sphere contact force and the contact area at maximum overlap are taken as reference values. The blue curve with circular points represents the sphere solution. As for the same force ratio, the ratio of the contact area for the contact situations spheroids lies always right of the sphere relation, it is obvious that spheroids form a larger contact area than spheres for the same contact force. This trend is more pronounced with increasing aspect ratio and more noticeable in the case of a pole contact between two oblate ellipsoids. Based on the relation between contact force and contact area for idealized configurations of two particles together with the observed texture that occurs as a result of compression, the steeper increase in the effective conductivity of the solid phase for spheroids with Figure 27. Relative frequency pðA mean =Ā mean Þ of weighted average mean cross-sectional areas A mean =Ā mean of pore phase at maximum compression ε z ¼ 20%.  www.advancedsciencenews.com www.entechnol.de increasing aspect ratio in Figure 25 and 26 can be explained. The larger coordination number, which is present in spheroidal particle assemblies and was found in other studies before, see the previous study, [55] has an additional positive impact to this trend. As was demonstrated in this section, significant changes may be induced in the electrode microstructure due to mechanical impact, such as the calendering step. The mechanics themselves are strongly influenced by the particle shape. As the mechanics, in turn, impact on the effective transport properties, big differences in the preservable effective transport properties between spheroids and spheres can be observed. Given the fact that the redox reaction, which takes place during loading and unloading of a battery, is an interplay between ionic and electronic transport, an optimal level of compaction has to be found to efficiently enhance the cell performance and to ensure sufficient conductivity characteristics for both phases. In this context, the influence of the particle shape on the mechanics, as well as on the effective transport properties, as has been demonstrated here, should be kept in mind for an ideal battery electrode design.

Conclusion
The major goal of this study was to investigate the influence of the particle shape on mechanics and on effective transport properties in granular LIB electrodes. To this end, all necessary equations for the implementation of superellipsoidal shapes into a mechanical DEM code were presented at the beginning of this work. This especially included the definition of superellipsoids and the formulation of the contact forces as well as the contact point definition. In addition, a microstructure generator based on the RCP algorithm for the generation of initial particle systems with superellipsoidal shapes was introduced. It was shown that random, highly packed but overlap-free initial states for superellipsoidal particles systems could be generated with the modified algorithm. An analytical formula for the calculation of the resistance between two mechanically compressed ellipsoids was derived in the theoretical part of this article, as well, whereas a geometry-based method was presented to determine the resistance of an individual pore throat in a superellipsoidal particle system. Based on these individual resistances, it was shown how a RN can be build up for the solid or the pore phase of a particle assembly to determine the effective transport properties of such structures. The derived formulas were verified by FEA and proved to deliver realistic values for the fast estimation of the effective conductivity of granular media by the comparison with well-known theories in the literature.
Based on the presented theoretical formulas, the influence of the ellipsoidal shape on the mechanical behavior and the effective transport properties under compression, such as in the calendering step during manufacturing of a LIB electrode, was investigated. Therefore, different random monosized ellipsoidal particle structures were generated with the newly RCP and uniaxially compressed. In the scope of this contribution, different mechanical and structural quantities during compression were examined to better understand the response of the particle system to the external impact. For that purpose, the stress evolution and the effective transport properties of the solid and pore phases with increasing compression were studied. Significant differences in the evolution of these properties were observed. In general, spheroidal particle systems showed a less stiff stress response at the same level of compaction. While a faster decrease in the pore phase conductivity with increasing compression and particle aspect ratio was observed, the effective conductivity of the solid phase showed a reverse trend. Here, the effective conductivity increased stronger with increasing aspect ratio. The different response behavior of spheres and ellipsoids could be explained with the help of structural investigations, such as the analysis of the particle's orientation, the pore throat structure, or theoretical bounds such as the maximum possible packing density in dependence of the particle's aspect ratio. Based on this work, deeper insight into the electrode behavior during calendering with the presented simulation techniques becomes possible, which, with regard to the requirement of specific cell properties, could be helpful to choose an appropriate electrode design.