Spatially resolved lithium‐ion battery simulations of the influence of lithium‐nickel‐manganese‐cobalt‐oxide particle roughness on the electrochemical performance

This publication deals with the influence of the roughnesses of single lithium‐nickel‐manganese‐cobalt‐oxide‐based active material particles on the electrochemical behaviour. A roughness model is used to create particles with a defined roughness depth. Subsequently, comparative calculations are performed using a 3D spatially resolved electrochemical model under different operating conditions. The results show that for small C‐rates below 1C, the usage of the widespread assumption of smooth and spherical active material particles in models is justified. For higher C‐rates, the particle's shape and roughness are not negligible. Small, rough particles show the best performance characteristics for high rate applications. For operation in the diffusion‐limited state (low diffusion coefficient and high C‐rate), the smallest absolute diffusion length represents a better correlation in terms of utilizable capacity compared to the volume‐specific active surface area.


| INTRODUCTION
Lithium-ion batteries (LIBs) are used as energy storage devices in many portable devices such as smartphones, laptops and power tools. They are also being used more and more in the field of electric mobility, such as in electric bicycles, electric cars and local public transport, and as stationary energy storage devices to buffer the fluctuating energy conversion in renewable energy systems. [1][2][3] Thanks to their high energy density and high rate capability, as well as their increasing cyclical and calendrical aging resistance, LIBs are still outperforming other battery technologies. 4 Nevertheless, there is a demand for ever more powerful mobile energy storage devices. For this reason, the further development of LIBs with the aim of increasing their power and energy density is a focus of research.
The positive electrode of a LIB consists of a heterogeneous particulate system, whereby the active material (AM) absorbs lithium into or releases lithium from its crystal lattice structure by intercalation or deintercalation during the discharge or charging process respectively. 5 Commercial LIBs consist of stacked or wound positive electrode layers, which amount 30 to 45 μm of thickness for high-performance cells 6 and 70 to 320 μm for highenergy cells. 7 Although the thickness of these layers is very small, the electrode plane in total covers an area of several cm 2 to m 2 . Nevertheless, the performancedetermining transport and intercalation mechanisms in the micro-and nanometer range take place at the particle and sub-particle scale. However, these processes are only experimentally accessible with considerable metrological effort. At this point, spatially resolved numerical simulations at the particle level are particularly suitable for expanding and further investigating the understanding of the processes taking place. 8 Starting at Garcia et al, 9 the area of completely spatially resolved microstructure simulation was developed, which is based on continuum theory. Due to the insufficient computing capacity at that time, the 2D-resolved simulation was carried out first, using circular AM particles. 9 Subsequently, Garcia and Chiang 10 extended this model to three dimensions with smooth spherical particles. Goldin et al 11 developed a three-dimensional particle-resolving model for the parameterization of reduced-order models. They used regularly arranged cut, smooth and spherical active particles. Less et al 12 published an experimentally parameterized, spatially resolved model consisting of smooth cubic AM particles and derived a correlation between the electrode morphology and the macroscopic battery performance. Wiedemann et al 8 performed 3D simulations based on reconstructed electrodes and showed that significant local fluctuations in the electrode structure occur, which cannot be represented with idealized microstructures, such as spherical AM particles. Mai et al 13 investigated the influence of microstructure and discharge rate on AM loading and mechanical stresses using smooth, spherical and randomly shaped particles. Hutzenlaub et al 14 modelled a conductive carbon black-binder domain in addition to the electrolyte and AM domain and performed simulations with reconstructed electrode structures to obtain a more realistic transport of charge and species. Roberts et al 15,16 numerically investigated the electrochemical, mechanical and thermal behaviour of LIBs using reconstructed, real electrode structures. The focus was mainly on mechanical degradation. Kespe et al [17][18][19] developed a 3D model taking into account periodic electrode structures consisting of smooth and spherical particles, whereby different particle sizes, the electrode structure, the carbon black distribution and blended cathodes were used to derive recommendations for optimized manufacturing processes. Mistry et al 20 investigated the complex interaction of the electrochemically coupled transport processes using spherical AM particles and the effect of the carbon black morphology and electrode porosity on battery performance. Rucci et al 21 pursued a simulative approach to correlate the correcting variables in electrode manufacturing with the battery performance, also using smooth and spherical particles.
From the literature review it can be seen, that up to now 3D simulations have mainly dealt with smooth spherical or idealized AM particles and their collectives or reconstructed real electrode structures have been used. The latter represent only a small section of the entire electrode plane, which has a negative effect on the representativeness of the results for the entire electrode. Reconstructions of real electrode structures, such as those published by Wiedemann et al, 8 Hutzenlaub et al 14 and Roberts et al, 15,16 show that the morphology of the active particles deviates significantly from the ideal smooth spherical shape. Furthermore, the investigations of Chung et al 22 showed that the surface roughness in real lithium-manganese-oxide electrodes is 2.5 times higher than that of ideally smooth and spherical particles, with the advantage of high available power. This publication continues the research on this topic by providing a systematic simulation based investigation of the influence of the roughness of AM particles on the electrochemical performance of batteries.

| Electrode particle generation
To investigate the influence of roughness and the volumespecific active surface area on the utilizable capacity, a computer-aided generation of particles with a defined roughness depth is applied. For this purpose, the roughness model proposed by Rumpf et al 23 is used. Those investigated the most important mechanisms of adhesion between particles or between a particle and a plate respectively. Part of this study was to determine the influence of surface roughness on the adhesive forces. They modelled the roughness as a small sphere on a larger sphere or plate.
For the intended investigations in this paper, lithiumnickel-manganese-cobalt-oxide (NMC) is used. NMC is a commercially available and very frequently used AM on the market. Dreizler et al 24 characterized commercially available NMC. The almost spherical AM particles of densely packed primary crystals had an internal porosity of less than 4%. Thus, the assumption of rough, pore-free particles is considered justified. Furthermore, experimental investigations show that NMC particles typically have primary particles in the range between 0.2 μm and 3.6 μm and a volume-related mean diameter d 50.3 of the secondary particles between 8 μm and 12 μm. [24][25][26] The primary particle size thus also determines the roughness of the particles.
For the generation of the differently rough particles, a purely geometric approach is used. The roughness radii R R decrease in relation to the Feret radius R F starting from R R = R F /4 by half to R R = R F /64. A smooth, spherical particle with Feret radius R F and roughness radius R R = 0 μm serves as reference. Figure 1 shows the definition of the rough particles using the roughness model proposed by Rumpf et al 23 The radius of the inner sphere R I , which also represents the smallest absolute diffusion length of the particle and on which the centres of the roughness spheres lie, is calculated by the difference of the Feret radius R F and the roughness radius R R (1).
The calculation of an approximate equal distribution of roughness spheres on this inner sphere is done using the Fibonacci-Lattice method. The angle θ i of each roughness sphere i describes the azimuth angle of the spherical coordinates with values between 0 and 2π and is calculated using (2). The golden ratio 1 + ffiffi ffi 5 p À Á =2 (middle expression) serves as the basis for the Fibonacci numbers, which is why this approach is called the Fibonacci-Lattice method. 27 The angle φ i describes the pole angle of the spherical coordinates and can have values between 0 (north pole) and π (south pole). It is calculated from (3), via the consecutive roughness sphere index i {1, …, n} ∈ Ν and the total number of spheres n. 27 The sphere count n results from geometrical considerations and can be found in Table 1 with further geometrical properties of the regarded AM particles. For simplification, the particles are named according to their roughness as an example instead of R R = R F /4 with R/4. The smooth sphere is always designated with R. Table 1 shows that the geometry named R/16 has the largest active surface area A. However, the volume V of the respective particles decreases linearly with increasing roughness, starting with the smooth sphere. On the other hand, the volume-specific surface area S V , increases non-linearly with increasing roughness of the particles considered here. Figure 2 shows the six AM particles with a Feret diameter of 10 μm, which are examined for their performance in this study. The smooth sphere serves as a reference.

| Electrochemical model
In this publication, the electrochemical model published by Kespe et al [17][18][19] is used. It consists of two nonoverlapping domains and considers the charge and mass transport in a spatially resolved positive electrode. The anode is assumed to consist of metallic lithium. The behaviour of the anode and the cathode current collector are modelled by appropriate boundary conditions.
The electron conduction in the rough NMC particles is described by the charge conservation equation in (4). Here, κ S represents the electrical conductivity and ϕ S the electrical potential.
The lithium transport in the AM particles is modelled via solid diffusion (5). c S represents the lithium concentration and D S the solid diffusion coefficient. In this article, only isotropic transport properties are used.
The electrolyte, assumed to be liquid, is modelled using the theory of concentrated solutions, taking electroneutrality into account (6). 28,29 Here, c E is the lithium concentration, D E the diffusion coefficient and i In equation (7), ϕ E stands for the electrical potential of the electrolyte and R for the universal gas constant.
T describes the absolute temperature and κ E the ionic conductivity. The electrochemical kinetic relationship of the Butler-Volmer type 28,30,31 in (8) is used for coupling the solid electrode and the electrolyte. The exchange current density i ! BV depends on the material-specific kinetic constant k BV , the concentrations at the particle surface of electrolyte and AM Γ S, E , the maximum concentration of lithium in NMC c S,max and the reference concentration c E,ref in the electrolyte. The anodic or cathodic apparent transfer coefficient is indicated by α a or α c respectively, 28 and the overpotential by η.

T A B L E 1
Properties of the computer-generated particles with a Feret radius of 5 μm

Name of particle
The overpotential η is calculated according to (9) from the difference between the electrical potentials of the solid ϕ S as well as the electrolyte ϕ E on the particles surface and the material-specific, concentration-dependent equilibrium potential U eq (c S ).
The continuity condition in (10) states that the ionic current density i The electrochemical reactions at the surface of the lithium metal anode are assumed to be fast, which is why the electrical potential there is set to 0 V. To satisfy the charge conservation equation over the entire electrode, the absolute ionic current at the anode is equal to the electric current at the current collector. The cathode current collector Γ CC is impermeable to lithium ions, correspondingly the ionic current density there in the normal direction to the arrester is i The material-specific model parameters of NMC, the electrolyte and the Butler-Volmer kinetics used in this study are summarized in Table 2.
For the course of the equilibrium potential U eq of NMC, the correlation of Stewart et al. and Kindermann et al. is used 30,31 (11).
The dependence of the ionic conductivity κ E of the electrolyte on its concentration c E is implemented by means of the correlation proposed by Less et al 12 (12).
To consider the concentration dependence of the diffusion coefficient D E (13), a correlation published in the same paper by Less et al. 12 is used.

Parameters of Butler-Volmer kinetics
Cathodic apparent transfer coefficient α c 0.5 -35 Anodic apparent transfer coefficient α a 0.5 -35 Butler-Volmer reaction rate constant k BV 2.895Á10 −7 A m mol −1 31 An important battery parameter for the evaluation of the discharge curves is the depth of discharge (DOD). In the spatially resolved model, this is given for the electrode domain Ω S in (14). 19 In order to determine the performance of a simulated battery half-cell, a further parameter was published by Kespe et al, 19 which is called "utilizable capacity" (UC), see (15).
In this context the terminal voltage also referred to as cut off voltage (COV), which is COV = 3.2 V for the investigations shown here.
The C-rate refers to the charge or discharge rate of a battery half-cell. 36 It is a quantity that makes different types of batteries and materials comparable with one and another. Furthermore, the C-rate is defined as the current that is necessary to fully charge or discharge the half-cell times the number of C-rate within 1 hour. In (16), i BV is the mean intercalation current density, A represents the active particle surface area and V the particle volume. Figure 3 shows a diagonally cut simulation result of the R/16 geometry at 4C and DOD of 0.49. From this, it can be seen that the cases are configured as half-cells with an active particle as positive electrode. The electrolyte domain is coupled via periodic boundary conditions in the plane. The assumption of isothermal operation at T = 298 K applies.
The model equations were implemented using the finite volume method in the open source simulation platform OpenFOAM of version 6. The simulations were conducted on the supercomputer ForHLR, which is funded by the Ministry of Science, Research and Art Baden-Wuerttemberg and the Federal Ministry of Education and Research. Figure 4 shows the course of the half-cell potentials of particle R/16 with Feret radius R F = 5 μm at different discharge rates from 0.1C to 10C over the DOD. It can be seen that as the C-rate increases, the DOD decreases until the COV at 3.2 V is reached. This has a direct effect on the UC, see (15). The reason for that lies in the intercalation current density, which increases with the C-rate, leading to a higher overpotential η on the particle surface according to (8). The latter can be seen from the vertical difference between the equilibrium potential and the half-cell potential in Figure 4, which rises with increasing C-rate.

| RESULTS AND DISCUSSION
F I G U R E 3 Exemplary representation of the particle R/16 at 4C and DOD = 0.49 using the numerical setup of the 3D structures employed in this simulation study. The shown structure is cut at the diagonal. The colouring of the AM particle represents the relative lithium concentration related to the maximum concentration of NMC (legend left), the differently coloured arrows represent the ionic current density within the electrolyte (legend right) In Figure 5, the respective UC of the different particles is plotted over the C-rate. It can be seen that with increasing C-rate the particles differ more and more in their achievable UC. The reason for this is that at low Crates the half-cells are close to thermodynamic equilibrium. For this reason, the particle morphology initially has only a small influence. At high C-rates, the importance of particle morphology takes into effect, which can be seen to a large extent in the results at the discharge rate of 10C. Here, the smooth sphere shows a relatively seen performance loss of 43.2% compared to the best AM particle R/4. One reason for this is that the intercalation current strength behaves proportionally to the active volume and thus different mean Butler-Volmer exchange current densities prevail at the particle surfaces. Relatively speaking, this current density in the smooth particle is 154.2% of that of the R/4 particle. The differences in UC grow with an increase in the C-rate from 0.7% at 0.1C to 6.2% at 1C and 24.9% at a discharge rate of 4C. This result allows the conclusion that the usage of smooth spherical AM particles in simulations of high-energy applications with low charge and discharge rates below 1C is justified. At higher discharge rates the simplified consideration leads to a significant underestimation of the UC. Figure 6 shows the concentration profiles of the investigated AM particles at an identical discharge state of DOD = 0.49 and a discharge rate of 4C. Here, the geometry R/4 has a 33.2% higher available capacity than the smooth one "R" at the same operating conditions. The better volumetric accessibility of lithium in the active particles is also evident from the smaller dimensions of the dark blue colouration inside the particles. The smaller volume of the particles with greater roughness due to the space between the roughness spheres and the lower intercalation current density on their active surface is also evident from the cut particles, as the particle surface is not as deep red in colour as in the particles with lower roughness.
In order to obtain reliable results, a scaling of the six initial geometries took place, resulting in 22 additional AM particles. In each case, fractions with 4 μm and 6 μm Feret radius, as well as particles with the same active volume as the smooth reference particle with 5 μm radius and with the same active surface as the smooth reference particle are generated. The relative particle roughness remains identical for the scaling of the particles, so that for the recognizability of the particle type, the naming remains the same despite different particle sizes. However, the different simulation cases are labelled and colour-coded in the diagrams.
In Figure 7, the usable capacity is plotted over the roughness radius R R of the corresponding AM particles. This figure shows that there are linear correlations between the roughness and the UC of each particle section with the same Feret radius (4 μm, 5 μm, and 6 μm), which extends from the smooth sphere (rectangle) to the particle with roughness R/4 (downward pointing triangle).The measures of determination of the linear regressions R 2 of these three case studies range from 99.78% to 99.88%. However, taking into account all single particle simulations, no direct dependence of the available capacity on the particle roughness can be identified. By enlarging the particles with the same structure, the roughness is increased by the factor x, the volume and thus also the overall intercalation current rate increases proportionally x 3 and the surface area, which correlates F I G U R E 4 Galvanostatic discharge curves at different discharge rates exemplarily for the particle R/16 with Feret radius R F = 5 μm. The course of the equilibrium potential of NMC published by Stewart et al 30 is shown as a dashed line F I G U R E 5 Utilizable capacities of the differently rough active material particles with Feret-Radius R F = 5 μm plotted over the applied C-rates with the intercalation current density, increases by x 2 . Thus the absolute diffusion length increases and the usable capacity deteriorates with increasing particle size, which is reflected in the vertical difference of the usable capacity of the respective particle type in Figure 7 (eg, rectangle, star, circle, etc.) from green R F = 4 μm (top) to blue R F = 6 μm (bottom).
In Figure 8, the UC is plotted over the radius of the inner sphere. From the diagram a linear relationship can be seen, which shows that for high discharge rates, such as the present 4C, the smallest absolute diffusion length F I G U R E 6 Concentration profiles of the considered AM particles at identical DOD of 0.49 F I G U R E 7 Representation of the usable capacity at a galavanostatic discharge rate of 4C applied over the roughness radius of the different AM particle geometries. The different colours depict the different simulation cases. At the top are the particles with constant Feret radius R F = 4 μm (green). Below that are the particles with the same active surface as the smooth sphere with R = 5 μm (grey). Then the particles with R F = 5 μm (black), the particles with the same volume as the smooth sphere with R = 5 μm (white) and finally the particles with R F = 6 μm (blue) follow F I G U R E 8 UC plotted over the inner radius R I of the particles. Here the colours show the different simulation cases. Beginning at the top left the particles with constant Feret radius R F = 4 μm (green) are depicted. Then the particles with the same active surface as the smooth sphere with R = 5 μm (grey), the particles with R F = 5 μm (black), the particles with the same volume as the smooth sphere with R = 5 μm (white) and the particles with R F = 6 μm (blue) follow to the particle centre is the relevant morphological parameter for the volume accessibility of lithium. This is an indication that the half-cells operate in the diffusionlimited state, which is in good agreement with the results of Du et al. 37 Thus, at least for NMC single particles, the smallest absolute diffusion length turns out to be a set screw for performance in high-power applications. It is to be expected that the herein examined AM particles with their roundish shape will exhibit, when used in representative, periodic half-cell electrode structures, qualitatively similar behavior like the single particles. Whether the latter is true and these deductions are also valid for active particles of different shapes is subject of ongoing research.

| CONCLUSIONS AND OUTLOOK
The results of the spatially resolved single particle simulations show that the roughness of AM particles has an influence on the electrochemical performance. A Crate study shows that for the simulation of high-energy applications with charge and discharge rates below 1 C the representation of electrode structures using smooth, spherical AM particles is justified. At a C-rate of 1C, the deviation between the particle with the greatest roughness depth and the smooth particle with a Feret diameter of 10 μm is less than 6.2%. At higher C-rates, the simplification of AM particles to smooth spheres leads to a significant underestimation of the utilizable capacity.
Looking at high performance applications, rough and small particles show the best performance. However, for operation in the diffusion-limited state, it is not the roughness depth or the volume-specific active surface that correlates with the available capacity, but rather the smallest absolute diffusion length towards the particle center. Due to small dimensions it is possible to lithiate a large part of the particle volume despite limiting conditions. Purely diffusion-limited operation prevails for the common particle sizes of about 10 μm and high discharge rates of 4C when using NMC according to Du et al. 37 Furthermore, due to the generation of the results using the material-specific properties of NMC, the findings cannot be directly transferred to other AMs. The different material data of diverse AMs, such as their diffusion coefficients and reaction rate constants, cause diffusion and kinetics limitations of different degrees throughout their cycling.
To what extent the results of the simulations shown are transferable to representative NMC electrode structures is the subject of ongoing research.