Towards a 3D-resolved model of Si/Graphite composite electrodes from manufacturing simulations

Chaoyue Liu, Oier Arcelus, Teo Lombardo, Hassan Oularbi, 1,2 and Alejandro A. Franco Laboratoire de Réactivité et Chimie des Solides (LRCS), CNRS UMR 7314, Université de Picardie Jules Verne, Hub de l’Energie, 15 rue Baudelocque, 80039 Amiens Cedex, France 2 Réseau sur le Stockage Electrochimique de l’Energie (RS2E), Fédération de Recherche CNRS 3459, Hub de l’Energie, 15 rue Baudelocque, 80039 Amiens Cedex, France ALISTORE-European Research Institute, Fédération de Recherche CNRS 3104, Hub de l’Energie, 15 rue Baudelocque, 80039 Amiens Cedex, France Institut Universitaire de France, 103 Boulevard Saint Michel, 75005 Paris, France


Introduction
Lithium ion batteries (LIBs) are already an important part of our everyday life that have enabled the rise of portable electronics. Looking in the short-, middle-term, LIBs also hold tremendous potential to substitute the use of fossil fuels in transportation, with the still incipient Electric Vehicle (EV) market expansion [1][2][3]. They have also found applications in decentralizing energy storage in residential areas, and they are being adopted for partially solving grid-unbalance issues caused by fluctuations in energy generation of renewable power sources [4,5]. However, increasing the energy density of LIBs remains a primary challenge for enabling some of the aforementioned applications. Focusing on the anode side, graphite is the most used material, thanks to its low price and good cycling stability [6]. Nevertheless, the maximum specific capacity that can be theoretically reached with graphite is 372 mAh/g, which limits the range of use cases for LIBs. In these regards, alloy negative electrodes, and in particular silicon (Si) [7], stands out as a promising candidate for next generation LIB anodes. When fully lithiated, Si forms the Li15Si4 alloy, showing a high theoretical specific capacity of 3579 mAh/g, ten times larger than graphite.
It also shows an appropriate voltage (~0.4 V) for being adopted as an anode. However, Si undergoes a ~300% volume expansion upon charging, which induces Si particle pulverization, unstable solid electrolyte interphase (SEI) formation and contact loss between Si and conductive additives because of internal stresses [8,9]. All these issues cause major problems for the performance and lifetime of LIBs using Si as anode. Some of these issues, such as the pulverization of the Si particles, can be solved through nanosizing [10], but much of the other problems that hinder the applicability of these anodes remain. Therefore, alternative strategies to incorporate Si in LIB anodes are currently being developed.
Today, much of the focus is now being given to composite electrodes made of graphite and nanosized Si (n-Si) [11,12], where graphite acts as an electrochemically active and mechanically stable framework to buffer the volume expansion of the n-Si particles. Such approach essentially relies on a compromise between cycling stability, provided by graphite, and higher specific capacities, provided by Si. Composite electrodes with low n-Si content (< 25 wt.%) are of particular interest, since they potentially enable manufacturing processes that are compatible with the current LIBs, while significantly increasing their specific capacity, when compared to graphite only electrodes [13]. For instance, as few as 10 wt.% of Si could, theoretically, contribute to approximately half of the electrode's capacity. Still, similar to electrodes made of Si only, the huge volume of expansion of the n-Si particles results in composite electrodes that, to some extent, also suffer from pretty significant microstructural changes. These changes result in electrode thickness variations [14], loss of electrical particle contact within the microstructure [15,16], continued consumption of electrolyte [17], all of which affects the electrode's electrochemical performance [18]. All those phenomena depend on the exact particles location along the electrode, where the heterogeneous distribution of electrode components in 3D becomes important when trying to understand the coupled mechanical and electrochemical behavior of graphite/Si composite electrodes [19,20].
Experimental 3D probing techniques, such as X-ray tomography [19], or the more recently developed scanning spreading resistance microscopy (SSRC) [20], are critical when trying to extract fundamental insights about the microstructure-performance relationships of battery electrodes in general, and of graphite/Si composite in particular. However, those techniques are too costly to be used systematically, and often require not routinely accessible instrumentations, such as synchrotron beamlines. In addition, those techniques often need to find compromises between the appropriate resolution/phase contrast needed to resolve certain phases, such as the carbon binder domain (CBD) [21,22], and the sample size, which can hinder microstructural representativity [20].
In these regards, computational 3D modelling has been proven to be an efficient tool to better understand battery phenomena at the mesoscale [23][24][25][26][27][28][29], where relationships between the electrode microstructure and performance can be elucidated. Concerning the generation of 3D electrode microstructural models, we can distinguish between two main approaches, each with their own advantages and disadvantages: (i) stochastic generation, where particles are randomly placed or piled on a simulation box until a certain target is met (porosity, particle number, etc.) [30,31] and (ii) physics-based simulations, such as Coarse-Grained Molecular Dynamics (CGMD) [32][33][34] and Discrete Element Method (DEM) [29,32,[35][36][37]. A recently developed third option consists of generating realistic 3D electrode microstructures by generative adversarial networks (GANs) [38,39]. Those models can then be fed to electrochemical models in order to assess their performance as a function of their microstructural properties [24,27,40].
Returning to our graphite/Si composite electrodes, it should be noted that Si is considered to be a challenging material to model, due to its severe volume change and potential hysteresis upon lithiation/delithiation. Many modelling efforts have been done in establishing the connection between the electrochemical and mechanical behavior of Si [41][42][43][44][45][46]. Some of these works decouple both phenomena [41,45], while others empirically relate the Si potential to its inner mechanical stress [42][43][44]. Modeling graphite/Si composite electrodes is even more challenging, and it is currently understudied from a computational perspective. In the work of Sturm et. al. [47], the electrode is treated as a single empirical phase, where properties such as open-circuit voltage (OCV) and maximum solid lithium concentration are considered as effective parameters.
Concerning the mechanical behavior of the composite electrodes, most of the studies are done at the experimental level, focusing mainly on macroscopic observables such as the thickness change of the electrode upon lithiation [14,48,49]. As a result, the influence of Si volume expansion on the composite electrode microstructure still lacks investigation, and this is especially true in modelling, where the 3D location of each electrode component should be resolved explicitly due to the tight link between degradation mechanisms and microstructure.
This work aims at moving a step forward in the modelling of graphite/Si composite electrodes.
Here, we report a new 4D (3D + time) model of a graphite/Si composite electrode, where the graphite and Si electrochemistry and mechanics are coupled. In particular, we show how the electrode microstructure changes upon galvanostatic discharge as a result of the interplay between Si and graphite. The pristine electrode microstructure, including graphite, Si particles and CBD is resolved with CGMD, similarly to what has been previously reported by our group [32,33].

Generation of the microstructural model
The generation of the graphite/Si electrode microstructures is done in several steps. First, a slurry that only contains graphite and CBD particles is simulated by using the experimental composition of a graphite slurry and the corresponding active material (AM) particle size distribution. However, in this case 10% of the total graphite mass is removed, to be later substituted by Si particles. Then, a drying simulation is performed in order to obtain a tentative dry graphite electrode. In this case, the reference composition for the graphite electrode was 94. Here, the graphite particles are considered explicitly, while the additives are considered implicitly in the CBD phase, whose particles intend to represent additive agglomerates. The experimental graphite particle size distribution is shown in Figure S1 in the supplementary information (SI). The graphite particles are approximated as spheres, which is different from the standard flat ellipsoidal graphite, but is nonetheless consistent the experimental shape of the GHDR graphite that we used ( Figure S2). Water was used as a solvent and the slurry solid content (SC) was 38%. Differently from our previous model, in which CBD particles accounted for carbon+binder+solvent at the slurry phase and carbon+binder only at the electrode phase [32], here the CBD particles account for carbon+binder only for both slurry and electrode phases, while the solvent is modelled by using devoted particles, whose total volume is calculated to match the experimental SC. On the one hand, the main idea behind the use of particles to model the solvent is to use highly overlapping particles, thus mimicking a continuum medium. On the other hand, the rationale behind using explicit CBD particles in the slurry is the consideration that carbon black particles and binder form aggregates in the slurry, which further densify upon drying, leading to CBD agglomerates at the electrode phase [50][51][52][53]. In principle, this would imply using smaller CBD particles (mimicking carbon+binder aggregates) at the slurry phase and mimicking their agglomeration along drying. In this work, for simplicity, and given that the objective of this work is not assessing the dynamics of CBD agglomerate formation, but studying the microstructure effects of Si expansion in graphite/Si composite electrodes, the CBD particles are kept constant with a diameter of 1.5 μm and a density of 0.866 g/cm 3 . This corresponds to CBD agglomerates with a nanoporosity of ~50%, as measured by previous FIB-SEM measurements [54]. In this type of CGMD simulations, all constituent particles in the model interact with each other with the Lennard-Jones (LJ) and Granular-Hertzian (GH) force fields (FFs), as introduced by earlier works [32,34].
After the slurry phase is formed, a dry electrode is generated by taking the slurry phase as the input structure and by removing the solvent particles. Then, an additional CGMD simulation is performed in order to equilibrate the dry electrode, by considering 'harder' LJ and GH parameters that mimic the formation of a more solid system, as explained in our previous work [32]. Once the graphite dry electrode is formed, the structure is decorated with the remaining 10 wt.% of Si particles and further equilibrated with CGMD simulations. This last step contains two further approximations: the first one concerning the Si particle size and the second one related to the role of n-Si in affecting the slurry phase and the associated drying.
On the one hand, the size of the Si particles we used (300 nm) is larger than the experimental one, since commercial Si nanoparticles range between 40 and 70 nm in diameter. This assumption is justified by the fact that Si nanoparticles hardly ever disperse perfectly in the dry electrode. Instead, they tend to form aggregates [55]. Additionally, this assumption allows us to build a model with a reasonable number of particles, given our computational resources. Figure S4 shows a preliminary calculation of the number of Si particles to be simulated as a function of the targeted simulation box volume and Si particles size, given the aforementioned formulation. As shown, the number of particles quickly rises when considering smaller Si particles, calling for a compromise between microstructural representativeness and the computational cost (which links to the total number of particles). This compromise is met by selecting a target volume of 30×30×30 μm 3 (even though, the exact final equilibrated volume will depend on the FF parameters) and a Si particle diameter of 300 nm, which results in system that contains ~90000 particles. Even if the Si particles where considered as n-Si aggregates in this work, we must note that a diameter of 300 nm is still a higher limit for the resistance to fracture upon volume expansion [7].
On the other hand, this model assumes that, due to the relatively low amount of added Si particles (10 wt.%), the slurry structure and associated drying (simulated by using graphite, CBD and solvent particles only) are still representative of the ones obtainable by explicitly considering the Si particles as well, which were added only after the drying step. This late addition of Si allows bypassing the full slurry and drying simulations, which would increase the computational cost by an order of magnitude. Of course, the rheological properties of the slurry without Si nanoparticles must not necessarily have to be similar to the ones of a slurry containing a small addition of Si nanoparticles. In fact, evidence suggests that the rheological properties of graphite/Si slurries differ when using μm-sized and nm-sized Si particles [56]. However, Si in dried electrodes is visible by common microscopy techniques such as Scanning Electron Microscopy (SEM), which allows reproducing its location ad hoc without the need of following their particle evolution from the slurry to the dried electrode. This is indeed not the case with CBD particles, thus the need of the full workflow to resolve their position [32].
Lastly, we note that the FF parameters for the CGMD model presented in this work are not fully validated to experimental viscosity curves (for the slurry), or porosity and density measurements (for the dry electrode). Instead, here we aim at qualitatively studying the effect of Si volume expansion on a given graphite/Si microstructure, for which obtaining highly experimentally validated microstructures is not of paramount importance. Nevertheless, using the above presented method to generate microstructures has various advantages when compared to cheaper stochastic microstructure generation. For instance, the model has the potential to be improved by fitting FF parameters to reproduce experimental physical properties. Additionally, DEM methods to simulate electrode calendering could be easily included in the workflow [29,32].
Summing up, Figure 1 shows the overall workflow used to generate the graphite/Si composite microstructures. The structures were equilibrated using a Nose-Hoover barostat and thermostat at 300 K and 1 atm. The simulations were performed using LAMMPS molecular dynamics software [57] using periodic boundary conditions in all directions. The FF parameters used for this work are included in Table S1. All simulations were performed in the TGCC Joliot-Curie (CEA) supercomputer, in 64 core AMD Rome (Epyc) processors. More details on the generation of CGMD electrode structures can be found in Ref [58]. The color and size of each particle is shown in the bottom left. The size of the Si particle is exaggerated for improved visibility. The state of charge map is represented in grayscale (red-blue) for graphite (Si) phases, respectively. The stress field represents the stress of the CBD network.

Meshing of the microstructure
After the graphite/Si microstructural model is generated, the CGMD structure is first voxelated with an in house program and converted to image stacks with voxels sizes of 0.25×0.25×0.25 μm 3 , a little bit smaller than the Si particles. First, it was realized that due to the overlaps between the graphite particles, the meshed structure resulted in a rigid graphite matrix that would not displace the particles as a result of the Si expansion. Therefore, in order to guarantee the free movement of graphite particles, the overlap regions where removed by making sure that at least a 0.5 μm gap existed between them, as shown in the meshed electrode in Figure 1. The AM particles were then dilated in order to compensate for the lost volume due to the cutting.
Due to the small size of the Si particles, and in order to avoid the generation of an excessively dense mesh, which would result in prohibitive simulation times, the Si particles were merged together in a homogenized domain. Additionally, in order to guarantee the connectivity of the CBD network, which is not obvious from the particle-like assumption of the CGMD simulations, the CBD phase was further treated with a 'connect-erode' process, as illustrated in Figure S3. In short, CBD particles are first connected with a cylinder of the same diameter as the particles themselves and the interface between the CBD and the electrolyte phases is eroded in order to recover the original volume. Lastly, a tetrahedral finite element mesh was generated from the voxelated microstructers, with the open-source Iso2mesh software [59].

Continuum Model
The mesh was then imported into COMSOL Multiphysics 5. voltage (OCV) profiles must be considered for charge and discharge. In this work, as single discharge process is analyzed, thus the corresponding Si OCV curve, as well as the graphite OCV curve were considered separately, as shown in Figure S5.
The outputs from the electrochemical model were then coupled to a mechanical model. In the mechanical model, the electrolyte phase was excluded because it was approximated as a free deformation domain, assuming that its effect on the deformation of the electrode was negligible.
Additionally, since the Si phase is approximated as a porous electrode, the inner stresses that are generated can no longer be directly related to the stresses between the explicit Si particles. Then, the equation of motion takes the form in stationary condition as following, where is Cauchy stress tensor: 0 = ∇ • Moreover, graphite, Si and CBD domains are treated as a linear elastic material. Hooke's law is used to describe the relation between the stress and strain tensors: where is the elasticity tensor, is the total strain tensor, and is the inelastic strain tensor.
Here, contains the information about the swelling of the graphite and Si upon lithium intake.
The formula that relates the displacement field and the strain tensor is written as: As already mentioned, both graphite and Si undergo volume expansion with lithium insertion. In fact, even if the volume expansion of graphite (10%) is much smaller than Si (~300%), its contribution to the whole electrode thickness change is non-negligible, since the majority of the AM volume fraction is occupied by graphite [14]. We assumed that the volume of change of both materials follows a linear relationship with the lithium concentration within the particle [60]: where and , are the volume per mole of host material ( = graphite, Si) at the current, and initial states, respectively. is the mole fraction of lithium in the host material, and Ω is the partial molar volume of lithium. The large volume changes of Si also result in large surface area changes, which considerably alter the reaction current densities and the Si overpotential.
Therefore, the Si surface area change was taken into account in the electrochemical model through the following formula: where, ! " is the surface area of Si, with units of m -1 , and $ ," is the initial radius of the Si particles.
Then, the inelastic strain tensor, , is also evaluated as: Additional assumptions for the coupled electrochemical-mechanical model correspond to the fact that the contact interface between the electrode and the current collector is fixed. Also, we note that the mechanical analysis of this work mainly focuses on the electrode microstructure changes. This is to say, even if we applied the mechanical properties of Si to a homogenized Si porous phase, the calculated stresses and strains cannot be directly related to the interactions between the Si particles. Thus, we only focus on the effects of the expansion of the Si phase on the displacements of the rest of electrode components. Additionally, interface between CBD and graphite (or Si) is also fixed. Therefore, no explicit detachment of CBD will occur as a result of active material volume expansion, even though high stress regions will be created.

Results and discussions
We studied the discharge of a half cell model ( Then, graphite starts to dominate, where the slope of the Si SOC curve starts to decrease. This twostage process can be further verified as shown in Figure 2(c), where the current contributions of each AM to the total current intersect in the middle of the discharge, therefore showing how each phase reacts at different discharge depths. There are several aspects that affect this behavior. Figure   S6 shows the average overpotentials calculated for graphite and Si, throughout the discharge process. We notice that the average value of the Si overpotential is maintained at very low values between -1.5 mV and -0.3 mV. On the contrary, graphite experiences a much larger overpotential ranging between -60 mV to -20 mV. In principle, this leads to much lower local current densities for the Si phase, when compared to graphite ( Figure S6). However, due to the fact that the Si particles are considered as nanosized (45 nm), given their extremely high surface area (~90000 μm 2 ), the Si phase actually accommodates rather large reaction currents. On the other side, graphite particles are much larger, and therefore their surface area is much lower (around 8400 μm 2 ), thus resulting in overall lower total currents, where graphite dominates. Looking at this process in more detail, we see that at around 25% of discharge depth, the current contribution curves intersect, and graphite starts dominating the reactions at the second stage. At this transition point, the SOC of Si and graphite is of around 0.5 and 0.1, respectively, which relates to the first plateau in the graphite OCV curve. Afterwards, the graphite OCV drops below 0.2 V and thus starts further reacting and increasing its SOC. Additionally, we can differentiate how each phase experiences the faradaic current, by calculating their respective C rates ( Figure S7). In fact, we can see that Si is experiencing C rates as high as 2.08 at 9% discharge depth. While the highest value for graphite reaches 1.33 at 54% discharge depth. Same for 0.1C scenario, the materials are also experiencing different C rates. This deviation from the C rate that is applied to the composite electrode as a whole, evinces the difficulties when working with electrodes where multiple AM phases coexist, and demonstrates the need of additional considerations to make in the design of efficient composite anodes.
The situation gets more complicated when considering smaller C rates. Comparing to the twostage process that we discussed at 1C, we can see that at 0.1C the SOC curves intersect at around 60% of the discharge depth (Figure 2(b)). At this point, the relative concentration of lithium within graphite becomes larger than that of Si. Looking at the evolution of the current contributions, we can see in Figure 2(c) that the behavior of the electrode is considerably more irregular at 0.1C.
This interweaved profiles show a complex interplay between the AM phases in the composite electrode that cannot be easily interpreted. However, a rather simple realization is that, due to the low applied current, the overpotential that is generated from the Faradaic reactions is omittable. In turn, the electrode potential of both phases will basically follow their respective OCVs. By revising the SOC curves, we can see that their intersection points coincide with those of their respective OCV curves, where each slope changes corresponding to a phase transition of graphite. Therefore, we can see that the potential curve of the composite electrode somehow resembles the shape of the graphite OCV. 1C. This lithium crosstalk between graphite and Si phases is recognized as an additional degradation mechanisms for the composite electrodes, and has recently been observed experimentally [61]. Which further strengthens the validity of our reported model.
Mass transport: Another contributor to the electrode overpotential is the mass transfer limitation within the liquid phase. Figure 3(a) shows the heterogeneous local current densities throughout the electrode, where differences between graphite and Si phases are apparent again.
The heterogeneity of the SOC is also clearly visible from Figure 3(b), where for both materials, the particles on the separator side (top) always show a higher SOC. At the end of the charge, the differences between the SOC of graphite and Si is of 0.81 and 0.37, respectively. Thus, we can clearly say that the overpotential also arises from the lithium concentration gradient in the liquid phase, which is shown in Figure 3(c). Indeed, the concentration difference can be as high as 700 mol/m 3 , which largely affects the lithium exchange current density as included in the Butler-Volmer equation, causing low reaction current densities in low lithium concentration zones. When using lower applied currents, the SOC distribution appears to be more homogeneous along the thickness of the model electrode, where a fully lithiated state is reached for both materials.
Mechanics: As already mentioned, during discharge, both graphite and Si will swell when lithiated, which results in an overall expansion of the electrode (Figure 4(a)).  Considering the 10 wt.% Si content of our model, and no additional degradation mechanisms, such as SEI formation or binder detachment, we think that our results are within a reasonable range. Figure 4(b) shows the evolution of the electrode thickness upon discharge, as we can see, the thickness of the electrodes increases faster when the Si is dominating the reaction. Additionally, small variations exist between 1C and 0.1C, where the latter results in slightly bigger thickness changes, given the higher SOC of the Si phase. Furthermore, we also found the expansion degree of the AM to be inhomogeneous. Figure S8 shows that the volumes ratios of the solid phases at the end of discharge vary across the electrode thickness. It is evident that, at 1C, due to the SOC distribution caused by mass transfer limitations, leads to solid volume ratios that are larger near the separator side. detachment from AM particles due to cohesion failure, which leads to a degradation of the structural integrity of the electrode, and a loss of electrical conductivity due to a damaged conductive percolation network. In order to give a closer look at these effects, we also analyzed the averaged pressure and Von Mises stress profiles along the thickness of the electrode. As can be seen in Figure 5(b) there is a clear correlation between the pressure felt by the CBD and its Von Mises stress, which is an indicator of mechanical failure of the CBD phase. Along the whole electrode thickness the CBD feels a negative pressure where, due to the expansion of the graphite and Si phases, and the electrode thickness change, it is being elongated. This stretching effect is particularly strong in the middle of the electrode, as can also be visually seen from the higher concentration negative pressure (blue) regions in Figure 5(a). Interestingly, we see that the minima of the pressure profiles coincide with the maximum points of Von Mises stress, where again the maximum points of potential CBD failure are located in the middle of the electrode. Additionally, the Von Mises stress also shows particularly high values near the CC and Separator sides, even if pressure does not, indicating more potential degradation effects on these interfaces.

Conclusions
In this work, we have presented the first steps towards a 3D modelling approach to capture manufacturing-performance relationships of graphite/Si composite electrodes. We presented a workflow for the generation of the graphite/Si electrodes using Coarse-Grained Molecular Dynamics (CGMD), where the location of the Si particles was resolved by an ad hoc decoration approach, followed by a relaxation run. Then, these microstructures were used as inputs for simulating the electrochemical performance of the graphite/Si composite electrodes by coupling electrochemistry and mechanics, given that Si suffers a ~300% volume expansion upon lithiation. we show that the model captures some of the basic observables for microstructural evolution upon Si and graphite volume expansion, such as the evolution of the electrode thickness. We also assessed the effects of stress and pressure on the electrode's CBD phase. The model captures potential failure points for the CBD phase, by looking at regions where Von Mises stresses are high. This high stress zones are correlated with negative pressure zones, indicating the stretching of the CBD phase. We argue that this potential failure points may cause detachment of the CBD phase from the active material particles, thus damaging the conductive percolation network. We also observe potential failure points at the CBD and CC (Separator) interfaces.
Summing up, all the aforementioned points have been a subject of discussion in many experimental works in the literature. We believe and we hope that the present model will help obtaining a more precise, and spatially resolved, understanding on the complex interplay between the electrochemical and mechanical behavior of composite graphite/Si anodes, eventually guiding the design of optimal electrode microstructures, potentially unlocking next generation of lithium ion batteries.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper