# Discrete modeling of the calendering process for positive electrodes of Li-ion batteries

Max Sonzogni, Jean-Mathieu Vanson, Yvan Reynier, Sébastien Martinet, Katerina Ioannidou, Farhang Radjai

► To cite this version:

Max Sonzogni, Jean-Mathieu Vanson, Yvan Reynier, Sébastien Martinet, Katerina Ioannidou, et al. Discrete modeling of the calendering process for positive electrodes of Li-ion batteries. Powder Technology, 2024, 442, pp.119891. 10.1016/j.powtec.2024.119891. hal-04722526

HAL Id: hal-04722526

https://hal.science/hal-04722526v1

Submitted on 5 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Discrete modeling of the calendering process for positive electrodes of Li-ion batteries

Max Sonzogni<sup>a,b</sup>, Jean-Mathieu Vanson<sup>b,*</sup>, Yvan Reynier<sup>c</sup>, Sébastien Martinet<sup>c</sup>, Katerina Ioannidou<sup>a</sup>, Farhang Radjai<sup>a,*</sup>

## HIGHLIGHTS

* We introduce a new discrete element model for calendering Li-ion battery electrode.

* This model involves a contact law accounting for both active and CBD materials.

* Calendering speed affects the elongation of the electrode and its internal anisotropy.

* Calendering leads to vertical tensile contacts and horizontal compressive contacts.

* Electronic and ionic conductivities are controlled by connectivity and porosity.

## GRAPHICAL ABSTRACT

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Keywords:**

Li-ion battery

Calendering

Discrete Element Method

Cohesion

Contact law

Fast Fourier Transform

## ABSTRACT

The electronic properties of Li-ion batteries crucially depend on the microstructure of their electrodes. One step of the manufacturing process, called ‘calendering’, consists in compressing the electrodes between two counterrotating cylinders to increase their density. Through a new simulation model, we investigate the effect of calendering on microstructural and electronic properties of the electrodes. Our model takes into account the real geometry of the rotating cylinder and a new contact law involving the elastic behavior of the active material and the cohesive-plastic behavior of the binder layer. Our results align well with experimental data for porosity, final thickness, and elongation. We show that the bonding structure induced by calendering involves mostly vertical tensile contacts and horizontal compressive contacts. This unexpected observation highlights the importance of shear deformation induced by rolling and thickness reduction. Using the FFT method, we also investigate the ionic and electronic conductivities of the numerically calendered electrodes.

<footnote><sup>a</sup> LMGC, CNRS, University of Montpellier, Montpellier, F-34090, France</footnote>

<footnote><sup>b</sup> CEA, DES, IRESNE, DEC, Saint-Paul-lez-Durance, F-13108, France</footnote>

<footnote><sup>c</sup> University Grenoble Alples, CEA, Liten, DEHT, Grenoble, F-38000, France</footnote>

<footnote>* Corresponding authors.
E-mail addresses: jean-mathieu.vanson@cea.fr (J.-M. Vanson), franck.radjai@umontpellier.fr (F. Radjai).</footnote>
# 1. Introduction

Li-ion batteries have proven their worth for many years in several domains and play a key role for various technologies [1–4]. The manufacturing process of their electrodes involves several stages designed to create a controlled material capable to store and deliver large amounts of energy. The conventional Li-ion cell manufacturing process involves a slurry preparation, coating, drying and calendering steps for their electrodes [5]. The calendering step concerns the compaction of the electrode between two rotating cylinders to reduce its thickness. This compaction of the electrode increases its energy density and enhances its storage capacity. However, this process also decreases the ionic conductivity of the electrode and thus the ability of the battery to store and deliver energy quickly [6,7]. The reduction of porosity increases the tortuosity and leads to the clogging of the pore network. The link between calendering parameters and the final microstructure and performance of batteries is therefore a crucial feature of manufacturing Li-ion batteries that remains, however, poorly understood [8].

In this paper, we are interested in the link between the final properties of the positive electrode and the calendering parameters. Since the Li-ion electrode has a granular microstructure, we employ the well known Discrete Element Method (DEM) to compute the mechanical behavior of the material during the calendering process [9]. This method is based on explicit representation of the active material particles and their interactions and has been recently employed by several authors to simulate the manufacturing process of electrodes [10–14]. These DEM-based models use different representations of the electrode within the algorithmic constraints imposed by the framework of DEM.

A natural DEM approach proposed by some authors consists in introducing two types of particles, namely the NMC active material and CBD (Carbon Binder Domain) [15]. The CBD is usually composed of sub micron carbon conductive additive particles and PVDF binder. This approach implies, however, a high computational cost as it requires hundreds of small CBD particles for each active material particle. Computationally affordable models are based on the representation of active material as rigid particles while the CBD is taken into account through contact interactions between particles. For example, the simplest DEM models of the electrode use a Hertzian-bond contact model to account for the elastic behavior of the active material particles and an adhesion force accounting for the action of the binding material [16]. This approach considers a solid bond formed by the binder during the drying step between the active material particles. Some models have improved this representation by adding a plastic component to the contact law, which allows one to simulate the damage of active material particles [10,11,17,18]. Some other models assume that the binder induces a constant attraction force between the active material particles and the particles follow an elasto-plastic behavior [19–21]. The use of irreversible cohesive force by some authors in this context implies irreversible loss of cohesion, which contradicts the chemo-physical nature of CBD.

Another aspect of DEM-based calendering models concerns the simulation of the external forces applied to the system. The loading applied on the material during the calendering process is usually represented by a triaxial compression test with a confining pressure applied on the lateral walls and vertical compression at constant speed. This simplified representation is based on the fact that the diameter of the calendering roll is much larger than the length of its contact with the electrode. The contact zone between the calendering roll and the electrode is therefore assumed to be flat, neglecting the curvature of the roll and thus, also, the resulting shear stress. Several studies have also been reported in which the compression is modeled with periodic boundary conditions along horizontal directions and vertical compression between two platens [12,19,20]. Since the horizontal cell dimensions are fixed, such a configuration is equivalent to uniaxial compression and cannot account for the free spreading of the material and its lateral extension during calendering. The shear stresses induced by friction

with the cylinders are also neglected like in triaxial compression and the contact surface is assumed to be flat. However, during the industrial process, the motion of the electrode between two counter-rotating rollers is driven by friction forces mobilized between the rollers and the surface of the electrode. Hence, shear stresses induced by friction and rolling are intrinsic to the process and might not be neglected. The effect of these shear stresses is probably less crucial for understanding the process at low calendering speeds. The models using simplified boundary conditions can therefore be relevant for low calendering speeds (< 1 m/min) usually employed in laboratory setups while the speeds used in industry are usually close to 10 m/min. The calendering speed affects horizontal elongation [22,23], and, beyond a certain level of elongation, the collector may break during winding [24].

In this work, we present a new DEM-based model for the simulation of the calendering process of Li-ion battery electrodes. The CBD material is represented by its adhesive-plastic behavior at the interface between NMC particles and most parameters are calibrated from experiments. We also model the calendering process by introducing a cylinder that drives the electrode by its rotation. We carry out a detailed parametric investigation by varying thickness reduction ratio and calendering speed and analyze their effect on porosity reduction, elongation, microstructure, and electronic properties of the calendered zone. Our results will be compared to previous DEM simulations reported in the literature when possible and, as we shall see, the full modeling of the calendering process and contact interactions turns out to be crucial for the simulated behaviors.

In the following, we first introduce in Section 2 the numerical method with focus on the interaction laws, calendering process, system calibration, boundary conditions, and FFT method for the computation of electronic conductivities. In Section 3, we analyze in detail the time evolution of system-scale and particle-scale variables and the influence of system parameters on the final state (porosity, elongation, microstructure, and conductivities). We conclude in Section 4 with the most salient findings of this work and its potential extensions.

# 2. Materials and methods

## 2.1. Discrete element method

To simulate the behavior of the granular electrode, we use the Discrete Element Method (DEM) [9]. Each particle is considered as a rigid element with its own shape, size, mechanical properties, and degrees of freedom. The equations of dynamics for the translational and rotational degrees of freedom are integrated in time for all particles according to a time-stepping scheme. The forces and moments resulting from the interactions between particles are computed by means of contact force laws. The total contact force $f$ between two particles is the sum of a normal force $f_n$, where $n$ is the contact normal, and a tangential force $f_t$:

$$f = f_n n + f_t. \qquad (1)$$

The normal force law is assumed to be a function of the overlap $\delta_n$ while the tangential force is defined as a function of the cumulative elastic tangential displacement $\delta_t$ at the contact point. The position and velocity of each element are then incrementally updated through the equations of motion. We used an in-house code called Rockable [25, 26]. Classical contact detection procedures and a velocity-Verlet time-stepping scheme are used in this code [27].

## 2.2. Contact force laws

Numerical models for electrodes either explicitly represent both active material (AM) and Carbon-Binder Domain (CBD) particles or only the particles of the active material and the action of the binding material is modeled through a contact force law with its plastic and cohesive components [28]. We consider the latter approach with
<figure><img src="image_5.png" /><figcaption>Fig. 1. Normal force $f_n$ as a function of overlap $\delta_n$ and schematic representation of the binder-covered particle model of the electrode. The colored zones represent the contact zones between binder layers and between particles, respectively.</figcaption></figure>

carbon-binder phase described implicitly as a layer surrounding the particles of active material. The AM particles define a polydisperse assembly of spheres [29] and the CBD is a layer of thickness $R_p$ coating AM particles, as schematized in Fig. 1. The thickness of the binder layer is fixed from the ratio of the mass $m_{CBD}$ of the binding material to the mass $m_{AM}$ of active material. The thickness of the surrounding CBD layer is kept constant during the calendering process. This yields two equations. The first equation is based on the porosity of the electrode $\epsilon$ computed as $1 - (V_{AM} + V_{CBD})/V_{tot}$, leading to the following expression:

$$ \epsilon = \frac{4\pi}{3V_{tot}} \sum_{k=1}^{N_p} (R_k + R_p)^3. \qquad (2) $$

The second equation is based on the mass distribution between active material and CBD:

$$ \frac{m_{AM}}{m_{CBD}} = \frac{\rho_{AM}}{\rho_{CBD}} \frac{\sum_{k=1}^{N_p} R_k^3}{\sum_{k=1}^{N_p} (R_k + R_p)^3 - R_k^3}. \qquad (3) $$

Based on this representation of binder-covered particles, there are three types of contacts between two particles depending on the overlap $\delta_n$: binder-binder, binder-particle, and particle-particle. We assume an elasto-plastic behavior with adhesion for binder-binder and binder-particle contacts and an elasto-plastic behavior for particle-particle contacts. Fig. 1 shows the evolution of normal force $f_n$ as a function of the overlap $\delta_n$. Upon the initial contact at point A between two particles, the interaction between the two binder layers is governed by a linear elastic law:

$$ f_n^{AB}(\delta_n) = k_1 \delta_n, \qquad (4) $$

where $k_1$ is the reduced stiffness of binder-binder and binder-particle stiffnesses. Note that overlaps are counted as positive. At point B the AM cores of the particles touch each other ($\delta_n = R_p$). The AM particles

being more rigid, the elastic interaction is governed by a larger stiffness $k_2$ when $\delta_n > R_p$ and we have

$$ f_n^{BC}(\delta_n) = k_1 R_p + k_2 (\delta_n - R_p). \qquad (5) $$

The active material particles used in Li-ion battery electrodes have relatively low yield strain and undergo damage and cracking under low compression [10]. We account for the resulting plastic behavior by introducing a yield strain $\delta_e$ (point C). Beyond this point, the normal force evolves following a strain-hardening power-law behavior:

$$ f_n^{CD}(\delta_n) = k_1 R_p + k_2 (\delta_e - R_p) + k_2 (\delta_n^e \delta_e^{1-\zeta} - \delta_e), \qquad (6) $$

where the exponent $\zeta$ controls the intensity of plastic hardening.

Upon unloading from any point, the normal force follows a linear elastic path, starting from the largest overlap $\delta_n^{max}$ reached during loading. Since the deformation of the CBD is irreversible, the elastic stiffness of the unloading path is that of the active material ($k_2$). The evolution of the normal force is then given by:

$$ f_n^{load}(\delta_n) = f_n^{max} + k_2 (\delta_n - \delta_n^{max}), \qquad (7) $$

where $f_n^{max}$ is the normal force reached at $\delta_n^{max}$ along the loading path. The adhesion of the binder is modeled by allowing the normal force to take negative values down to a pull-off force $f_p$ depending on both binder properties and the maximum overlap reached. The cohesive bond breaks if this tensile path is continued.

Since adhesion depends on the amount of binding material between two particles, we assume that the absolute value of $f_p$ increases linearly with the area of the contact surface. Geometrically, the contact area can be approximated by $2\pi\delta_n^{max}R^*$ as far as $\delta_n^{max} < R_p$, where $R^* = (1/R_1 + 1/R_2)^{-1}$ is the reduced radius. Hence, in this range $f_p$ increases linearly in absolute value with $\delta_n^{max}$. The largest contact area between the binding layers occurs for $\delta_n^{max} = R_p$ and we therefore assume that the $f_p$ is constant for $\delta_n^{max} > R_p$. This is shown by a dashed line in the tensile force domain of Fig. 1 and can be written as

$$ \begin{aligned} f_p(\delta_n^{max}) &= -2\pi R^* \sigma_a \min(\delta_n^{max}, R_p) \\ &= -k_p \min(\delta_n^{max}, R_p), \end{aligned} \qquad (8) $$

where $\sigma_a$ is the yield strength of the binder. The highest bonding force $f_p(R_p) = -2\pi R^* R_p \sigma_a$ is the adhesion force given by the DMT theory [30]:

$$ f_a = -2\pi\gamma R^*, \qquad (9) $$

where $\gamma = \sigma_a R_p$ is the surface energy. This value is reached if the binder layer is fully crushed, i.e. when the contact surface area between the binder phases of the particles reaches its maximum value and the loading is tensile.

For the tangential force, we use a linear elastic law combined with a Coulomb dry friction criterion:

$$ f_t = \begin{cases} -k_t \delta_t & \text{if } \|f_t\| \le \mu f_n^e \\ -\mu f_n^e / \|\delta_t\| & \text{otherwise} \end{cases}, \qquad (10) $$

where $\mu$ is the friction coefficient, $k_t$ is the tangential stiffness and $f_n^e = f_n - f_p(\delta_n^{max})$. Note that when contact slip occurs, the value of $\delta_t$, which represents only the cumulative elastic tangential displacement keeps its last value reached before slip. If the contact becomes non-slipping again, $\delta_n$ will be cumulated from its last value, i.e. the value reached before the slip begins. In other words, the elastic energy of the tangential spring is conserved.

Throughout this work we set $\mu = 0.4$, which is a common value of friction coefficient, and $k_t$ is computed following the Hertz-Mindlin theory for frictional contacts [31]:

$$ \frac{k_t}{k_2} = \frac{2(1 - v^*)}{2 - v^*}, \qquad (11) $$
where $v^* = (1/v_1 + 1/v_2)^{-1}$ is the equivalent Poisson ratio. Compared to the classical Coulomb criterion $\|f_f\| \le \mu f_n$ for cohesionless contacts, here the Coulomb cone is shifted to account for the adhesion force present in the total normal force [31,32]. This means that only the repulsive part $f_n^e$ of the normal force comes into play for the tangential force computation.

In addition to the explicit plastic-adhesive behavior induced by the binding material, we introduce a normal damping force $f_n^v$ to the elastic and adhesion forces:

$$f_n^v = -2\alpha_n \sqrt{k_1 m} \delta_n, \quad (12)$$

where $m = (1/m_1 + 1/m_2)^{-1}$ is the reduced mass between two particles. The parameter $\alpha_n$ is related to normal restitution coefficient $\epsilon_n$ in the absence of adhesion and plastic behavior as follows:

$$\alpha_n = - \frac{\ln(\epsilon_n)}{\sqrt{\ln^2(\epsilon_n) + \pi^2}}. \quad (13)$$

In all simulations, we set $\epsilon_n = 0.45$, corresponding to $\alpha_n \approx 0.25$. We set the cohesionless tangential damping to zero, so that the energy dissipation is controlled by adhesion, explicit plastic behavior of the binding material, and implicit normal damping.

As compared with other models used for the simulation of the calendering process, the model presented in this section has the advantage of both accounting for the binding material and being computationally efficient. The spatial distribution of CBD between active material particles is in reality inhomogeneous. The assumption of a uniform coating of the active material particles by the binding material provides an alternative mean-field approach that accounts for its mechanical behavior and binding effect without resorting to the computationally expensive approach based on its explicit representation as a continuum matrix or a discrete material composed of small particles as in Ref. [15]. Furthermore, the thickness of the binding layer directly accounts for the volume of the binding material so that the binding volume effects are incorporated into the model while contact laws based only on the representation of the active material particles ignore such effects. In our model, the plastic behavior and adhesion are combined into a single law with the adhesion force depending on the history of the contact. This contact law is similar to that previously developed for plastic particles [33], but generalized here to a coated particle. Finally, it is also worth mentioning that all the parameters introduced in this model have physically clear experimental counter-parts and can therefore be calibrated by direct comparison with experiments.

## 2.3. Simulation of calendering

Our DEM model and its parameter values were based on the electrodes that were manufactured and characterized by means of an in-house setup performed in CEA-LITEN Grenoble. The cathode structures were made of NMC-811 (LiNi$_{0.8}$Mn$_{0.1}$Co$_{0.1}$O$_{2}$) – referenced thereafter as active material (AM) with density $\rho_{AM} = 4.80 \text{ g/cm}^3$, carbon black (CB) and polyvinylidene fluoride (PVDF), with mass proportions 96-2-2 of the three components. The last two components form the conductive CBD matrix of density $\rho_{CBD} = 2.00 \text{ g/cm}^3$. The polydispersity of the NMC particles had a log-normal distribution, with $d_{10} = 6.1 \text{ µm}$, $d_{50} = 12.1 \text{ µm}$ and $d_{90} = 21.8 \text{ µm}$. The slurry was coated onto an aluminium foil of thickness $h_{al} = 20 \text{ µm}$, resulting in electrodes of initial thickness $h_t = 2h_i + h_{al} = 194 \text{ µm}$. The initial porosity of the electrodes was $\epsilon_0 = 44.5\%$. The electrodes were calendered with different speeds and intercylinder gaps using a rolling press with cylinders of diameter $D_c = 20 \text{ cm}$. Disk samples of diameter 14 mm were extracted from the calendered electrodes and the thickness $h_f$ of each side was measured. By measuring the mass of the cylindrical samples, their porosity $\epsilon$ was computed from the volumes, densities, and proportions of the components. The measured values of porosity, thickness, and elongation are given in Table 1.

<table><caption>**Table 1**
Calendering parameters used to produce the experimental electrodes described in Section 2.3 and resulting structural electrode properties.</caption><thead><tr><th>Target thickness (µm)</th><th>Calendering speed (m/min)</th><th>Measured thickness (µm)</th><th>Porosity (%)</th><th>Elongation (%)</th></tr></thead><tbody><tr><td>Initial</td><td>-</td><td>87</td><td>44.5</td><td>0</td></tr><tr><td rowspan="3">65</td><td>1</td><td>66.5</td><td>27.5</td><td>1.11</td></tr><tr><td>1.5</td><td>68</td><td>28.8</td><td>0.99</td></tr><tr><td>2</td><td>68</td><td>29.2</td><td>0.87</td></tr><tr><td rowspan="3">62.5</td><td>1</td><td>65</td><td>26.2</td><td>1.48</td></tr><tr><td>1.5</td><td>65</td><td>25.4</td><td>1.36</td></tr><tr><td>2</td><td>64</td><td>25.0</td><td>1.23</td></tr><tr><td rowspan="3">60</td><td>1</td><td>62</td><td>23.4</td><td>1.72</td></tr><tr><td>1.5</td><td>64</td><td>24.2</td><td>1.60</td></tr><tr><td>2</td><td>63.5</td><td>24.2</td><td>1.48</td></tr></tbody></table>

<table><caption>**Table 2**
Values of the simulation parameters for particle-particle, particle-collector, and particle-cylinder contacts.</caption><thead><tr><th></th><th>Sph - Sph</th><th>Sph - Col</th><th>Sph - Cyl</th></tr></thead><tbody><tr><td>k<sub>1</sub> (kN/m)</td><td>2.5</td><td>4.5</td><td>4.7</td></tr><tr><td>k<sub>2</sub> (kN/m)</td><td>130</td><td>93</td><td>165</td></tr><tr><td>k<sub>p</sub> (kN/m)</td><td>4</td><td>27.5</td><td>0</td></tr><tr><td>γ (mJ/m<sup>2</sup>)</td><td>34</td><td>240</td><td>0</td></tr><tr><td>ζ</td><td>0.15</td><td>0.275</td><td>1</td></tr><tr><td>k<sub>t</sub> (kN/m)</td><td>117</td><td>77</td><td>138</td></tr></tbody></table>

For DEM simulations, we consider spherical particles with the same particle size distribution as in experiments. The value of $R_p$ was computed introducing the experimental data into Eqs. (2) and (3) and set to $R_p = 0.34 \text{ µm}$. The yield strain value $\delta_e$ was taken from Gimenez et al. [10] and set to $R_p + d_{50}/100 = 0.40 \text{ µm}$. The current collector was considered rigid and undeformable. This assumption holds only for relatively low/medium calendering degrees since high levels of compression tend to warp the current collector [18]. The calendering rolls are represented explicitly in the simulations. Since our experimental data were obtained using a stainless steel cylinder of radius $R_c = 10 \text{ cm}$, the length of our samples was set to $l = 1 \text{ cm}$. The width was set to 200 µm with periodic boundary conditions. With an initial thickness $h_i = 87 \text{ µm}$, we built a numerical electrode composed of 24 425 spherical particles. The values of the contact law parameters are given in Table 2. The values of the hardening parameter $\zeta$ were calibrated through a sensitivity analysis over the porosity and thickness of the sample. Since performing this preliminary study on the complete calendering model would require a relatively high computational cost, we used a uniaxial compression/decompression test, similar to those found in the literature [11–14]. The value of $\zeta$ was varied between 0.05 and 0.25 for the contacts between spheres and between 0.15 and 0.35 for the contacts between spheres and the current collector. The value of $\zeta$ for the contacts with the calendering roll was fixed to 1 since no plastic deformation of the cylinder is assumed. Three levels of compaction (thickness reduction) were studied and compared with the experimental data. Fig. 3 shows the results of porosity and final thickness for different values of $\zeta$ used for the contacts between spheres. We see that both porosity and final height increase with $\zeta$, since this parameter controls the hardening of the material and hence its maximal normal force under compression. The dotted lines in each plot correspond to the targeted experimental values. The values of $\zeta$ which lead to results close to these targeted values are located between 0.12 and 0.17. We therefore chose $\zeta = 0.15$ for the contacts between spheres, and 0.275 for the contacts with the current collector.

The simulations were performed in 3 steps. First, the particles were deposited on the current collector under gravity. The relaxed sample is shown in Fig. 2. Then, the cylinder is lowered slowly until it reaches the
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Macro view of the calendering simulation sample before lowering the cylinder (left) and snapshots of a portion of the electrode during the passage of the calendering roll (right). The current collector and calendering roll are colored in gray and black respectively, while the particles are color-coded relatively to average particle stress in the snapshots. The snapshots correspond to the initial and final states and crossover points between the four stages of evolution of porosity and pressure identified in Fig. 4(a).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Evolution of (a) porosity and (b) final height as a function of hardening parameter $ζ$ for contacts between particles in a uniaxial compression test. Each color corresponds to a compaction level and the dotted lines to the reference experimental values.</figcaption></figure>

specified position. Finally, the cylinder starts to rotate at an increasing rate until the target rotation speed is reached. The motion of the sample is thus driven by friction with the calendering roll. The translation of the current collector on Z and X axes are disabled so that the gap $h_g$

between the cylinder and current collector is constant. Fig. 2 shows close-up snapshots of the electrode at several time steps of simulation.

## 2.4. Computation of ionic and electronic conductivities

To compute the ionic and electronic conductivites of the samples, we employ the Fast Fourier Transform technique [34]. The FFT is a homogenization technique which is used to compute the effective or equivalent property of a multi-phase medium. It is extensively employed in mechanics, thermal transfer or diffusion. The multi-phase geometry of the sample is described explicitly on a cartesian meshing grid. Each voxel of the grid has a value corresponding to its local property. The diffusion gradient is considered as homogeneous inside the volume. The method assumes periodic boundary conditions and the effective properties of the matrix (mesh) is computed considering a steady state regime. Compared with finite-element method, the FFT is able to handle more degrees of freedom and may be easily accelerated thanks to parallelization (openMP here). Moreover, the cartesian mesh allows a very simple definition and easy to process geometry.

FFT computation needs a discretization of the bed defined on a cartesian mesh (voxelisation). It basically defines the property/value of a voxel thanks to its phase location. Most voxels contains only one phase where the property to assign is clearly defined. For the voxels containing two phases (interfaces) the property is not clearly defined. We employ the procedure developed in [35]. The ill-defined voxels need a special treatment. For a two-phase microstructure, the most simple strategy consist in assigning a defined property on each voxel (either phase 1 or phase 2). It gives a lower and an upper bound depending on the phase chosen for the ill-defined voxels. However, depending on the ratio between phase properties these bounds might be far from each other and thus far from the real value.

An alternative and “mean” procedure consists in assigning the property of the center of the voxel to the entire voxel. If the center of the voxel is in phase 1, we affect the property of phase 1, and if it is in phase 2, we affect the property of phase 2. This leads to a “mean estimate” which is an approximation. In fact, it does not rely on the volume proportion of phases in the voxels. However, we assume the error is averaged by the broad variety of ill-defined voxels in all configurations. The results from [35] confirms that this procedure provides a reasonable estimate. In our case, we deal with 3 phases. The rules are the same as that for 2 phases. For ill-defined voxels, we affect the properties of the center of the voxel. The effective property of the voxellized microstructure is then computed using the FFT.
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Time evolution of porosity (a), pressure (b), coordination numbers (c) and deviatoric fabric tensor components (d) for imposed thickness reduction $C_r = 0.28$ and calendering speed $v_{cal} = 2$ m/min. The labels of four successive stages of the evolution of porosity are marked. The coordination numbers $Z^+$ and $Z^-$ for compressive and tensile contacts, respectively, are also shown in (c). The deviatoric fabric tensor components in (d) are calculated with (dotted lines) and without (plain lines) accounting for contacts with the cylinder.</figcaption></figure>

# 3. Results and discussion

## 3.1. Calendering steps

We performed 28 simulations with different values of the gap $h_g$ between the cylinder and the current collector and the calendering speed $v_{cal} = R_c \omega_{cal}$, where $\omega_{cal}$ is the rotation speed of the roll. The values are summarized in Table 3. We investigate here the deformation of the electrode and the evolution of its microstructure in the ‘calendered zone’ defined as the portion of the electrode corresponding to the contact area between the cylinder and the material and which is fully compressed once it has moved a distance equal to its own length $L_c$ due to the rotation of the cylinder. Calendering is controlled by the thickness reduction ratio $C_r$ defined by

$$C_r = \frac{h_i - h_g}{h_i}, \qquad (14)$$

where $h_i$ is the initial thickness of the electrode.

Fig. 4 shows the time evolution of porosity, average pressure, coordination number and deviatoric fabric tensor components in the calendered zone for $C_r = 0.28$ and speed $v_{cal} = 2$ m/min. The initial porosity is 0.44. We distinguish four different steps:

1. Initial state: The cylinder is away from the calendered zone and the pressure is zero. The calendered zone is stable and fully relaxed under the action of its own weight and internal cohesive forces.

2. Compression: The average pressure increases in the calendered zone and porosity declines.

<table><caption>**Table 3**
Calendering gap $h_g$ and calendering speed $v_{cal}$ used in the DEM simulations.</caption><tr><td>$h_g$ (μm)</td><td>7 values ∈ [75, 50]</td></tr><tr><td>$v_{cal}$ (m/min)</td><td>{0.5, 1, 2, 5}</td></tr></table>

3. Relaxation: The cylinder leaves the calendered zone, which freely relaxes with a small elastic rebound.

4. Final state: A new stable state is reached with zero internal pressure and a porosity of 0.27.

In the following, we focus on the influence of calendering parameters on this process and final relaxed state reached by the calendered zone after calendering. Note that each data point represents the average value over 5 calendered zones selected from the electrode during a single simulation run, and error bars are their standard deviation.

## 3.2. Thickness, porosity, and elongation

Because of elastic rebound, the final thickness reduction ratio $C_f$ is not exactly equal to the imposed value of $C_r$. Fig. 5(a) displays the evolution of $C_f$ as a function of $C_r$ for different values of calendering speed $v_{cal}$. The simulations are in good agreement with experimental data, and we observe that $C_f$ is independent of $v_{cal}$. The relationship between $C_r$ and $C_f$ is linear with a shift of ≈ 2.5% between the target and measured values due to elastic rebound. This linear evolution and a constant shift from the perfectly plastic deformation $C_f = C_r$, independently of the value of $C_r$, is obviously a consequence
<figure><img src="image_9.png" /><figcaption>Fig. 5. System parameters as a function of thickness reduction ratio $C_r$ for different values of the calendering speed in simulations (circles) and experiments (triangle). Error bars represent standard deviation over several calendered zones from a single simulation run. (a) Effective thickness reduction $C_f$ of the electrode. The dashed line represents the ideal thickness reduction without elastic rebound ($C_f = C_r$). (b) Elongation of the electrode. (c) Porosity $\epsilon$. The red dashed line represents the fitting function given in Eq. (15). (d) Maximum stress induced by calendering. The stress is directly calculated on the calendering roll.</figcaption></figure>

of linear plastic behavior of the contacts defined in Eq. (7) and the vanishing of the pressure exerted on the electrode during calendering. The compression induces a high level of self-equilibrated forces since thickness reduction ratio is high. This comes from the fact that adhesion force increases with plastic surface area between particles. The elastic rebound is therefore mainly due to contacts between the cylinder and the particles. The maximum elastic displacement after the vanishing of pressure is $\approx f_a/k_1$ at the particle-cylinder contacts, leading to a deformation of $f_a/(k_1d_{50})$. With the values of parameters used in the simulations ($f_a = 1360$ µN, $k_1 = 4.7$ kN/m, and $d_{50} \approx 12.1$ µm), we find an elastic deformation of the order of 2.4%, which is consistent with what we measure both in simulations and experiments.

The reduction of the electrode thickness has two origins: (1) reduction of porosity and (2) elongation. Fig. 5(c) shows the final porosity $\epsilon$ as a function of $C_r$ for different values of calendering speed. As expected, $\epsilon$ declines as $C_r$ increases but does not depend on $v_{cal}$. The porosity has to tend to a constant low value for high compression level. The numerical data is well fitted by a power-law function:

$$ \epsilon = A (B + C_r)^{-\theta} + D, \qquad (15) $$

with $A \approx 0.17$, $B \approx 0.57$, $D \approx 0.06$, and $\theta \approx 1.5$. We have $\epsilon (C_r = 0) \approx 0.45$, which is close to the initial porosity $\epsilon_0 = 0.44$ of our numerical electrode. The predicted lowest value of porosity by this function is $\approx 0.15$ for $C_r = 1$. The lowest value of $h_g$ cannot be below one particle diameter, implying that the largest value of $C_r$ is $(h_i - d_{50})/h_i \approx 0.89$. For this value the porosity is around 0.16 which is well aligned with the fitting function.

Fig. 5(c) shows also that the numerical values of porosity are close to the experimental data. However, the experimental value of porosity

tends to decline faster with increasing $C_r$ than the porosity simulated. This discrepancy reflects the rather crude assumption in simulations that the binding material does not deform and remains attached to the particles as a plastic layer. In reality, the binder is unevenly distributed in the pore space and deforms with compression [29]. Furthermore, at high pressure levels, the active particles may also deform and break, allowing the material to reach even lower levels of porosity.

Since the electrode is compressed vertically and sheared horizontally, it is expected to expand, specially along its longitudinal y direction. Indeed, since in our simulations we applied periodic boundary conditions along the x direction, the electrode cannot expand along this direction. The longitudinal elongation $\Delta L/L$ of the electrode is shown in Fig. 5(b) as a function of $C_r$ for different values of calendering speed. The elongation is nearly proportional to $C_r$, and decreases with increasing calendering speed. This effect is unexpected since both final thickness and porosity does not depend on the calendering speed. As shown on Fig. 5(b), the experimental data show a similar correlation between elongation and calendering speed. It is noteworthy that none of the reported numerical calendering simulations available in the literature has been able to predict this dependence on calendering speed.

In our model, we have an explicit representation of the calendering cylinder, which sets the electrode in motion via the action of friction force on top of the electrode. The observed decrease of elongation at higher speeds might come from the shearing of the electrode as a result of the mobilization of friction force. If a relative slip of the cylinder with respect to the electrode occurs at higher speeds, the electrode is sheared less and the resulting elongation is smaller. Also, high calendering speeds result in shear stresses strong enough to either rip the
<figure><img src="image_10.png" /><figcaption>Fig. 6. Numerical microstructural parameters as a function of thickness reduction ratio C, for different values of the calendering speed (a) Final values of coordination number Z and (b) compressive and tensile coordination numbers Z+ and Z-. (c) Diagonal elements of deviatoric fabric tensor after calendering for calendering speed vcal = 2 m/min. (d) Vertical fabric component Fzz. The dashed line in (c) and (d) represents the isotropic state.</figcaption></figure>

granular mixture from the collector or directly warp it. Cathode current collectors are usually thin aluminium foils which cannot endure large deformations and tend to break during the assembly step [24]. This observation clearly shows the importance of shear stresses developed in the electrode during calendering and the need for the real geometry of the manufacturing process in numerical simulations.

The pressure acting on the calendered zone is a consequence of constriction imposed by the gap between the cylinder and collector. Fig. 5(d) shows the maximum stress σmax calculated on the calendering roll as a function of C. We see that σmax increases almost linearly with C, in exception to the largest values of C, and it is independent of calendering speed.

Fig. 7 shows the final porosity ε as a function of σmax. Considering the initial porosity at zero stress, the relationship between porosity and maximum vertical stress is not perfectly linear. It might be fitted to a power-law function when forced to pass by the initial porosity at zero stress:

$$ \epsilon = G (\sigma_{max} + H)^{-\beta} + I, \quad (16) $$

with G ≈ 9.5, H ≈ 140 MPa, I ≈ −0.26, and β ≈ 0.53. We have ε (σmax = 0) ≈ 0.43, which is the initial porosity of the electrode. This evolution following an inverse power law has already been observed with previous simpler calendering models [10,36,37]. This fitting function implies, however, that ε vanishes at a finite stress while we have seen that the fitting function (15) predicts a finite porosity. If we take 0.16 as the lowest porosity, as suggested by Eq. (15), the maximum pressure needed to reach this porosity is 212 MPa, above which the fitting form of Eq. (16) is unphysical. Further simulations are necessary to check the validity of the fitting forms proposed here for higher values of C.

### 3.3. Evolution of microstructure

DEM simulations provide access to microstructural variables such as the contact and force networks. We are interested in the parameters that control the electric conductivity across the granular microstructure. The lowest-order parameter is the coordination number Z, defined as the average number of contact neighbors per particle:

$$ Z = 2 \frac{N_c}{N_p - N_0}, \quad (17) $$

where N_c is the total number of contacts, N_p is the total number of particles, and N_0 is the total number of floating particles that have one or no contact neighbors. In a cohesive granular medium, contacts can either be compressive or tensile. Compressive contacts are those which have a positive normal force whereas tensile contacts have a negative normal force. By restricting the contact neighbors only to compressive or tensile contacts, we also define the coordination numbers Z+ for compressive contacts and Z− for tensile contacts with Z = Z+ + Z−. The proportion of tensile contacts is a descriptor of the stress state. When the external confining pressure is high compared to adhesion forces acting between particles, Z+ is larger than Z−.

Fig. 4(c) shows the time evolution of Z, Z+, and Z−. The coordination number increases as the cylinder approaches the calendered zone, reaches a maximum value slightly above 6, and then declines as the cylinder leaves the calendered zone. The final value of Z is much higher than its initial value. The compressive coordination number Z+ follows a similar evolution. The tensile coordination number declines and remains nearly constant during the rolling of the cylinder before relaxing to a higher value in the final state. It is noteworthy that before calendering, Z+ and Z− are nearly equal with a slightly lower value of
<figure><img src="image_11.png" /><figcaption>Fig. 7. Final porosity ε after calendering as a function of maximum stress σmax for different calendering speeds. The dashed line is a power-law fitting function given by Eq. (16) forced to pass by the initial porosity at zero stress.</figcaption></figure>

$Z^{-}$, such that $Z^{-}/Z \approx 0.45$. In the final state, the proportion of tensile contacts is $Z^{-}/Z \approx 0.43$, which is only slightly above its initial value.

Fig. 6(a) and (b) respectively shows the values of Z, and $Z^{+}$ and $Z^{-}$ after calendering as a function of thickness reduction ratio $C_r$ for different values of calendering speed. The latter has a small effect on average but it is not significant within the statistical precision of the data. The three coordination numbers increase linearly with $C_r$. Z increases from 4.1 before calendering to 6.4 for $C_r = 0.43$. Interestingly, the ratio $Z^{-}/Z$, i.e. the proportion of tensile contacts in the system is nearly constant ($\approx 0.43$) and independent of $C_r$. This shows that the stress state after relaxation and under gravity is the same and independent of $C_r$. In other words, the compression of the electrode due to calendering is large enough to drive the microstructure to a state which is independent of the initial state. As we shall see below, this final state reflects the anisotropic structure of the calendered material.

As compared to Z, which represents the average connectivity of the particles, the anisotropy of the contact network is a higher-order descriptor of granular microstructure. It is conveniently described by the fabric tensor F defined as

$$F_{ij} = \frac{1}{N_c} \sum_{k=1}^{N_c} n_i n_j. \qquad (18)$$

By definition, we have $tr(F) = 1$ so that its deviatoric part is given by

$$F'_{ij} = F_{ij} - \frac{1}{3} \delta_{ij}, \qquad (19)$$

where $\delta$ is the Kronecker delta. The deviatoric fabric tensor F' quantifies the relative deviations of the proportions of the contacts in each direction from the perfect isotropic state, in which the contact orientations are random and uniformly distributed in all space directions. By definition, we have $tr(F') = 0$. Hence, with respect to an isotropic distribution of contact orientations, a positive value of a component in a given direction reflects an excess of contacts whereas a negative value means a lack of contacts in that direction.

In the calendered zone, we consider two different tensors by either including or excluding the contacts with the cylinder. Fig. 4(d) displays the time evolution of diagonal components of F' calculated in the calendered zone for both of these tensors. The components $F'_{xx}$ and $F'_{yy}$ are almost equal. The electrode is initially in an anisotropic state with higher value of $F'_{zz}$ compared to $F'_{xx}$ and $F'_{yy}$ as a consequence of gravitational deposition used to build the electrode. As the roll advances in the calendered zone, $F'_{zz}$ increases due to the new contacts created between the electrode and the roll but it declines if these contacts are not included. At the same time, since deviatoric fabric

tensor is traceless, the two planar components $F'_{xx}$ and $F'_{yy}$ decrease when the contacts with the roll are included and increase otherwise. As the roll leaves the calendered zone, $F'_{zz}$ decreases and $F'_{xx}$ and $F'_{yy}$ increase. As a result, all components are lower in absolute value in the final state. In other words, calendering reduces the initial anisotropic structure of the electrode.

The decrease of $F'_{zz}$ during calendering is counterintuitive as vertical compression is expected to induce new contacts along the vertical direction thereby increasing $F'_{zz}$. This is what occurs when a granular material is subjected to triaxial compression. In contrast, in the calendering process the variation of fabric components suggests that new contacts are gained in the horizontal direction and lost in the vertical direction. In fact, the horizontal motion of the cylinder with an imposed thickness reduction, the horizontal force exerted by the cylinder on the top calendered layer, and mobilization of friction forces at the interface between cylinder and electrode induce a complex shear deformation that controls the gain and loss of contacts. The evolution of fabric tensor reflects this deformation and its evolution during loading and unloading, showing that the calendering process cannot be reduced to simple compression.

Fig. 6(c) shows the fabric components after calendering as a function of $C_r$ for $v_{cal} = 2$ m/min. We see that, independently of $C_r$, the two horizontal components are always nearly equal. All components decrease almost linearly in absolute value with $C_r$ and tend to zero at $C_r = 0.43$. A similar evolution was observed in other studies [10, 20]. Fig. 6(d) displays $F'_{zz}$ after calendering as a function of $C_r$ for different values of the calendering speed. In all cases, the anisotropy declines with increasing $C_r$ but we observe a slight dependence on the calendering speed. Increasing the speed leads to less reduction of $F'_{zz}$. Consistently with the effect of calendering speed on the elongation of the electrode, we may attribute this effect to the overall shear deformation of the material. Less elongation implies lower shear deformation although the vertical compression is the same for a given value of $C_r$, and lower shear deformation leads to less evolution of the fabric.

The stress and fabric states of the electrode before and after calendering can be visualized through the force network as displayed in Fig. 8. In the initial state, the contact forces are induced by gravity and the contacts are oriented around $\pm 45^{\circ}$ with respect to the vertical axis due to the initial gravitational deposition. After calendering we observe a large number of contacts oriented along the horizontal and vertical directions. Many vertical contacts are tensile whereas compressive contacts occur predominantly along the horizontal direction. This organization of the force network is consistent with the orientation of the fabric tensor. It is also important to note that both compressive and tensile forces are much larger than in the initial state. The larger values of forces is a consequence of higher mobilization of tensile forces by the action of calendering and self-balanced by compressive forces of the same order of magnitude.

### 3.4. Electric properties

The effective electronic and ionic conductivities of our numerically calendered electrodes are computed by means of the FFT method. We extracted 5 samples of size $200 \times 200 \times 50 \ \mu m^3$ from our simulations and used a homemade code to voxelize them. The voxel size was set to $(0.5R_p)^3 \ \mu m^3$ in order to represent efficiently the binder layer. For the spherical particles of size considered here and high contrast between conductivities ($> 100$), the error attributed to the voxel size (discretization) is expected to be below 5% regarding the sensitivity analysis made in Ref. [3]. The error increases rapidly with voxel size. Since FFT computations require the sample to be periodic in all directions, we did not include the current collector. The type of material attributed to each voxel is chosen from the DEM sample. If a voxel is on the edge of a particle, it is identified as CBD phase and if the voxel is inside a particle, it is considered as active material. In all other cases, the voxel is in the electrolyte phase.
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Force chains inside the electrode before (a) and after (b) calendering. Line thickness is proportional to force magnitude (×100 in the initial state).</figcaption></figure>

<table><caption>**Table 4**
Conductivity values used in the FFT computations.</caption><thead><tr><th></th><th>NMC</th><th>CBD</th><th>Electrolyte</th></tr></thead><tbody><tr><td>Electronic conductivity (mS/cm)</td><td>10<sup>-2</sup> [38]</td><td>100 [39]</td><td>10<sup>-5a</sup></td></tr><tr><td>Ionic conductivity (mS/cm)</td><td>10<sup>-4a</sup></td><td>10<sup>-2</sup> [40]</td><td>10 [41]</td></tr></tbody></table>

Each voxel has its own bulk conductivity corresponding to its phase. The values that we employed are summarized in Table 4. Some conductivity values – namely the electronic conductivity of the electrolyte and the ionic conductivity of the active material – are vanishingly small. However, too high contrasts between the highest and lowest conductivities in FFT computations lead to very high computational costs and slow convergence rates. We did a parametric study on one reference case with values in the range [10⁻⁷, 10⁻³] mS/cm in order to select conductivity values which reduce most the computational cost while keeping the effective conductivities of the electrode close to reference values.

Electronic conductivity $\lambda_{el}$ represents the flow of electronic current across the active material particles to the current collector. Fig. 9(a) shows $\lambda_{el}^{FFT}$ the electronic conductivity computed through FFT as a function of thickness reduction ratio $C_r$ for $v_{cal} = 1$ m/min. We observe a clear correlation between electronic conductivity and $C_r$. $\lambda_{el}$ increases as a function of $C_r$ along the three directions and seems to tend to a plateau value, although a clear plateau is not reached. The gain in electronic conductivity is nearly the same (≈ 30%) in the three directions. Nevertheless, the electronic conductivity along vertical direction is slightly above those in the other directions. Fig. 9(b) shows $\lambda_{el}^{FFT}$ as a function of $C_r$ along the vertical axis (Z) for different values of calendering speed. The calendering speed does not seem to affect vertical electronic conductivity while, as shown in Fig. 6(d), $F'_{zz}$ slightly depends on the calendering speed. In fact, the variation of $F'_{zz}$ as a function of speed is too small to affect significantly the electric conductivity.

The electronic conductivity depends on both contact network and electronic conductivity at each contact. The latter is a function of the contact area and varies therefore with normal force and the plastic deformation of the contact [11]. Dimensional analysis implies that $\lambda_{el}^{FFT}$ is

proportional to the conductivity $\lambda_{el}^{NMC}$ of NMC particles. Furthermore, the number density of contacts $n_c = Z\phi/2V_p$, where $V_p$ is the average volume of one particle and $\phi = 1 - \epsilon$ is the packing fraction, and their orientations through the fabric tensor control the amount of electric current and thus the conductivity of the network. Since $F_{ii} = 1/3 + F'_{ii} \approx 1/3$, we may neglect the effect of fabric anisotropy. Hence,

$$ \lambda_{el} \propto \phi Z \lambda_{el}^{NMC}, \qquad (20) $$

In Fig. 9(b) we have plotted $\lambda_{el}$ as a function of $C_r$ from this equation. We see that the prediction is quite good for all values of $C_r$ with a proportionality factor $\approx 0.27$ independently of $C_r$. This relation shows that the influence of $C_r$ is due to the variations of $\epsilon$ and $Z$.

Ionic conductivity $\lambda_{ion}$ reflects the diffusion of ions in the pore space between particles under the influence of chemical potential gradient. Fig. 9(c) shows $\lambda_{ion}^{FFT}$ calculated by FFT in three directions for $v_{cal} = 1$ m/min as a function of $C_r$. $\lambda_{ion}$ declines with increasing $C_r$ in all directions. Fig. 9(d) displays $\lambda_{ion}^{FFT}$ in the vertical direction for different values of the calendering speed. We see that $\lambda_{ion}^{FFT}$ is independent of calendering speed. This result is expected since the ionic conductivity mainly depends on the conductivity of the electrolyte and the porosity $\epsilon$ [21], which was observed to be independent of the calendering speed. We have

$$ \lambda_{ion}^{FFT} \propto \epsilon \lambda_{ion}^{electrolyte}. \qquad (21) $$

The fitting curve in Fig. 9(d) was based on this equation with proportionality factor $\approx 0.51$. We see that the fit correctly adjusts the data points.

## 4. Conclusion

We developed a DEM-based model for the simulation of the calendering process of Li-ion battery electrodes. The CBD material was taken into account as a thin layer coating NMC particles and governed by a plastic-adhesive behavior. Most parameters were calibrated from real experimental results. The process was modeled by introducing a cylinder that drives the electrode by its rotation via friction force mobilization at its interface with electrode. As a result, the calendered zone undergoes a complex deformation combining shear and compression. A parametric investigation was performed by simulations for a range of values of thickness reduction ratio and calendering speed. The effect of calendering and its parameters was analyzed in terms of porosity, elongation, microstructural parameters, and electronic properties of the electrode.

We found that most results on porosity and elongation of the electrode were in well aligned with experimental data. The elongation of the electrode and some other properties such as the vertical fabric component were shown to be slightly dependent on the calendering speed. This dependence was not previously observed in simulations of the literature. It shows the important role of explicit representation of the calendering roll and accounting for the calendering speed to correctly estimate the elongation of the electrode.

We showed that, as a result of contact plastic behavior and particle rearrangements, the elastic rebound is small, as in experiments. The porosity decreases as a power law with increasing thickness reduction ratio and seems to tend to a limit value representing the lowest reachable porosity. The highest stress reached in the calendered zone was found to increase almost linearly with thickness reduction ratio. As a result, the total coordination number and its tensile and compressive component increase linearly with thickness reduction ratio, resulting in much higher self-balanced tensile and compressive force chains. Interestingly, the proportion of tensile contacts is constant and independent of thickness reduction ratio, suggesting that the self-balanced structure induced by calendering is similar for all levels of thickness reduction. An important result shows that this structure involves mostly vertical tensile contacts and horizontal compressive contacts which is in radical contrast with the expectation that vertical compression tends to induce

<footnote>a Values are tuned to ease FFT computations while keeping consistent results as explained in text.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Electronic conductivity $\lambda_{el}^{FFT}$ (and ionic conductivity $\lambda_{ion}^{FFT}$ respectively) computed through FFT as a function of thickness reduction ratio $C_r$ (a,c) in all directions for $v_{cal} = 1$ m/min, and (b,d) along the vertical direction for different values of the calendering speed. The red dashed curves in (b) and (d) correspond to the functions defined in Eq. (20) and (21) respectively.</figcaption></figure>

compressive contacts along the vertical direction. This counterintuitive observation was attributed to shear deformation induced by rolling and thickness reduction.

The electronic properties of our numerically calendered electrodes were computed by means of the FFT technique. The effective electronic conductivity increases with thickness reduction ratio and is independent of calendering speed. The ionic conductivity decreases with thickness reduction and is also independent of calendaring speed. The electronic conductivity was shown to mainly depend on packing fraction and coordination number while ionic conductivity depends only on porosity.

With industrial applications going towards larger calendering rolls (>50 cm) and faster production lines (up to 100 m/min) in order to produce electrodes faster, this work can be extended to investigate the effects of larger calendering speeds and larger thickness reduction ratios. A more detailed analysis is necessary to quantify the deformation field in the calendered zone and the degrees of local slip and shear at the electrode-cylinder interface as a function of calendering speed. In this work, the rolling resistance at the contact point due to surface roughness was neglected since initial porosities close to their experimental values were obtained without introducing rolling resistance. It is known that the frictional strength of granular materials is an increasing function of rolling resistance [42,43]. We therefore expect that a major effect of rolling resistance would be the increase of stresses induced by calendering. It should be interesting to investigate numerically the effect of this parameter, which controls the stress but also the mobility of particles during calendering. The flexibility of the current collector should also be investigated and taken into account and the spatial repartition of the binder layer.

Nomenclature

<table><tr><td><i>f</i></td><td>Total contact force (N)</td></tr><tr><td><i>n</i></td><td>Contact normal</td></tr><tr><td><i>f</i><sub>n</sub></td><td>Normal force (N)</td></tr><tr><td><i>f</i><sub>t</sub></td><td>Tangential force (N)</td></tr><tr><td>δ<sub>n</sub></td><td>Normal overlap between two particles (μm)</td></tr><tr><td>δ<sub>t</sub></td><td>Cumulative elastic tangential displacement between two particles (μm)</td></tr><tr><td><i>R</i><sub>p</sub></td><td>Coating layer thickness (μm)</td></tr><tr><td><i>k</i><sub>1</sub></td><td>Reduced stiffness of binder-binder phase (kN/m)</td></tr><tr><td><i>k</i><sub>2</sub></td><td>Reduced stiffness of AM-AM phase (kN/m)</td></tr><tr><td>δ<sub>e</sub></td><td>Yield strain of the active material (μm)</td></tr><tr><td>δ<sub>n</sub><sup>max</sup></td><td>Largest overlap reached before unloading (μm)</td></tr><tr><td>ζ</td><td>Hardening exponent</td></tr><tr><td><i>f</i><sub>a</sub></td><td>Pull-off force (N)</td></tr><tr><td>σ<sub>a</sub></td><td>Yield strength of the binder (MPa)</td></tr><tr><td>γ</td><td>Surface energy of the binder (mJ/m<sup>2</sup>)</td></tr><tr><td><i>R</i></td><td>Radius of a particle (μm)</td></tr><tr><td><i>k</i><sub>t</sub></td><td>Tangential stiffness (kN/m)</td></tr><tr><td>μ</td><td>Friction coefficient</td></tr><tr><td>ν</td><td>Poisson ratio</td></tr><tr><td>α<sub>n</sub></td><td>Normal damping coefficient</td></tr><tr><td>ε<sub>n</sub></td><td>Restitution coefficient</td></tr><tr><td><i>h</i><sub>i</sub></td><td>Thickness of the mixture deposit before calendering (μm)</td></tr><tr><td><i>h</i><sub>f</sub></td><td>Thickness of the mixture deposit after calendering (μm)</td></tr><tr><td><i>h</i><sub>g</sub></td><td>Gap between cylinder and current collector (μm)</td></tr></table>
<table><tr><td>v<sub>cal</sub></td><td>Calendering speed (m/min)</td></tr><tr><td>d<sub>50</sub></td><td>Mean size of AM particles (µm)</td></tr><tr><td>ρ</td><td>Density (g/cm<sup>3</sup>)</td></tr><tr><td>C<sub>r</sub></td><td>Applied thickness reduction ratio</td></tr><tr><td>C<sub>f</sub></td><td>Final thickness reduction ratio</td></tr><tr><td>ε</td><td>Porosity</td></tr><tr><td>φ</td><td>Solid fraction</td></tr><tr><td>σ</td><td>Stress measured on the calendering roll (MPa)</td></tr><tr><td>σ<sub>max</sub></td><td>Maximal stress measured on the calendering roll (MPa)</td></tr><tr><td>Z</td><td>Total coordination number</td></tr><tr><td>Z<sup>+</sup></td><td>Compressive coordination number</td></tr><tr><td>Z<sup>-</sup></td><td>Tensile coordination number</td></tr><tr><td>N<sub>c</sub></td><td>Total number of contacts</td></tr><tr><td>N<sub>p</sub></td><td>Total number of particles</td></tr><tr><td>N<sub>0</sub></td><td>Total number of floating particles</td></tr><tr><td>F</td><td>Fabric tensor</td></tr><tr><td>F'</td><td>Deviatoric fabric tensor</td></tr><tr><td>λ<sub>el</sub></td><td>Electronic conductivity (mS/cm)</td></tr><tr><td>λ<sub>ion</sub></td><td>Ionic conductivity (mS/cm)</td></tr></table>

### CRediT authorship contribution statement

Max Sonzogni: Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Methodology, Investigation, Data curation. Jean-Mathieu Vanson: Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization. Yvan Reynier: Writing – review & editing, Validation, Supervision, Investigation. Sébastien Martinet: Writing – review & editing, Validation, Supervision, Investigation. Katerina Ioannidou: Writing – review & editing, Supervision, Methodology, Investigation. Farhang Radjai: Writing – review & editing, Supervision, Methodology, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### References

[1] K. Mizushima, P.C. Jones, P.J. Wiseman, J.B. Goodenough, Li₂CoO₂ (0<x<1): A new cathode material for batteries of high energy density, Mater. Res. Bull. 15 (1980) 783-789, http://dx.doi.org/10.1016/0025-5408(80)90012-4.

[2] E. Karden, S. Ploumen, B. Fricke, T. Miller, K. Snyder, Energy storage devices for future hybrid electric vehicles, J. Power Sources 168 (2007) 2–11, http://dx.doi.org/10.1016/j.jpowsour.2006.10.090.

[3] A. Kwade, W. Haselrieder, R. Leithoff, A. Modlinger, F. Dietrich, K. Droeder, Current status and challenges for automotive battery production technologies, Nature Energy 3 (2018) 290–300, http://dx.doi.org/10.1038/s41560-018-0130-3.

[4] L. Mauler, F. Duffner, J. Leker, Economies of scale in battery cell manufacturing: The impact of material and process innovations, Appl. Energy 286 (2021) 116499, http://dx.doi.org/10.1016/j.apenergy.2021.116499.

[5] J. Smekens, R. Gopalakrishnan, N.V. den Steen, N. Omar, O. Hegazy, A. Hubin, J. Van Mierlo, Influence of electrode density on the performance of Li-Ion batteries: Experimental and simulation results, Energies 9 (2016).

[6] H. Kang, C. Lim, T. Li, Y. Fu, B. Yan, N. Houston, V.D. Andrade, F.D. Carlo, L. Zhu, Geometric and electrochemical characteristics of Li(Ni1/3Mn1/3Co1/3)O2 electrode with different calendering conditions, Electrochim. Acta 232 (2017) 431–438, http://dx.doi.org/10.1016/j.electacta.2017.02.151.

[7] X. Lu, S.R. Daemi, A. Bertei, M.D.R. Kok, K.B. O'Regan, L. Rasha, J. Park, G. Hinds, E. Kendrick, D.J.L. Brett, P.R. Shearing, Microstructural evolution of battery electrodes during calendering, Joule 4 (2020) 2746–2768, http://dx.doi.org/10.1016/j.joule.2020.10.010.

[8] M. Duquesnoy, C. Liu, D.Z. Dominguez, V. Kumar, E. Ayerbe, A.A. Franco, Machine learning-assisted multi-objective optimization of battery manufacturing from synthetic data generated by physics-based simulations, Energy Storage Mater. 56 (2023) 50–61, http://dx.doi.org/10.1016/j.ensm.2022.12.040.

[9] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique 29 (1979) 47–65, http://dx.doi.org/10.1680/geot.1979.29.1.47.

[10] C.S. Giménez, B. Finke, C. Schilde, L. Froböse, A. Kwade, Numerical simulation of the behavior of lithium-ion battery electrodes during the calendaring process via the discrete element method, Powder Technol. 349 (2019) 1–11, http://dx.doi.org/10.1016/j.powtec.2019.03.020.

[11] C. Sangrós Giménez, C. Schilde, L. Froböse, S. Ivanov, A. Kwade, Mechanical, electrical, and ionic behavior of lithium-ion battery electrodes via discrete element method simulations, Energy Technol. 8 (2020) 1900180, http://dx.doi.org/10.1002/ente.201900180.

[12] C. Sangrós Giménez, L. Helmers, C. Schilde, A. Diener, A. Kwade, Modeling the electrical conductive paths within all-solid-state battery electrodes, Chem. Eng. Technol. 43 (2020) 819–829, http://dx.doi.org/10.1002/ceat.201900501.

[13] A.C. Ngandjong, T. Lombardo, E.N. Primo, M. Chouchane, A. Shodiev, O. Arcelus, A.A. Franco, Investigating electrode calendering and its impact on electrochemical performance by means of a new discrete element method model: Towards a digital twin of Li-Ion battery manufacturing, J. Power Sources 485 (2021) 229320, http://dx.doi.org/10.1016/j.jpowsour.2020.229320.

[14] T. Lombardo, A.C. Ngandjong, A. Belhcen, A.A. Franco, Carbon-binder migration: A three-dimensional drying model for lithium-ion battery electrodes, Energy Storage Mater. 43 (2021) 337–347, http://dx.doi.org/10.1016/j.ensm.2021.09.015.

[15] I. Srivastava, D.S. Bolintineanu, J.B. Lechman, S.A. Roberts, Controlling binder adhesion to impact electrode mesostructures and transport, ACS Appl. Mater. Interfaces 12 (2020) 34919–34930, http://dx.doi.org/10.1021/acsmi.0c08251.

[16] C.S. Giménez, B. Finke, C. Nowak, C. Schilde, A. Kwade, Structural and mechanical characterization of lithium-ion battery electrodes via DEM simulations, Adv. Powder Technol. 29 (2018) 2312–2321, http://dx.doi.org/10.1016/j.apt.2018.05.014.

[17] A. Lundkvist, P.-L. Larsson, E. Olsson, A discrete element analysis of the mechanical behaviour of a lithium-ion battery electrode active layer, Powder Technol. 425 (2023) 118574, http://dx.doi.org/10.1016/j.powtec.2023.118574.

[18] J. Zhang, J. Sun, H. Huang, Z. Yuan, Influence of calendering process on the structural mechanics and heat transfer characteristics of lithium-ion battery electrodes via DEM simulations, Particuology 85 (2024) 252–267, http://dx.doi.org/10.1016/j.partic.2023.06.015.

[19] D. Schreiner, A. Klinger, G. Reinhart, Modeling of the calendering process for lithium-ion batteries with DEM simulation, Procedia CIRP 93 (2020) 149–155, http://dx.doi.org/10.1016/j.procir.2020.05.158.

[20] R. Ge, D.J. Cumming, R.M. Smith, Discrete element method (DEM) analysis of lithium ion battery electrode structures from X-ray tomography-the effect of calendering conditions, Powder Technol. 403 (2022) 117366, http://dx.doi.org/10.1016/j.powtec.2022.117366.

[21] R. Ge, A.M. Boyce, Y.S. Zhang, P.R. Shearing, D.J. Cumming, R.M. Smith, Discrete element method and electrochemical modelling of lithium ion cathode structures characterised by X-ray computed tomography, Chem. Eng. J. 465 (2023) 142749, http://dx.doi.org/10.1016/j.cej.2023.142749.

[22] D. Mayer, A.-K. Wurba, B. Bold, J. Bernecker, A. Smith, J. Fleischer, Investigation of the mechanical behavior of electrodes after calendering and its influence on singulation and cell performance, Processes 9 (2021) http://dx.doi.org/10.3390/pr9112009.

[23] Y. Chen, M.B. Jaksa, B.T. Scott, Y.-L. Kuo, Numerical investigation of the performance of the 3-sided impact roller, Comput. Geotech. 157 (2023) 105331, http://dx.doi.org/10.1016/j.compgeo.2023.105331.

[24] H. Yasuhiko, M. Yoshiyuki, Nonaqueous electrolyte secondary battery and nonaqueous electrolyte secondary battery fabricating method, 2011.

[25] V. Richefeu, P. Villard, Modeling Gravity Hazards from Rockfalls to Landslides, Elsevier, 2016, http://dx.doi.org/10.1016/B978-1-78548-076-8.

[26] M. Sonzogni, J.-M. Vanson, K. Ioannidou, Y. Reynier, S. Martinet, F. Radjai, Dynamic compaction of cohesive granular materials: Scaling behavior and bonding structures, Soft Matter (2024) http://dx.doi.org/10.1039/D3SM01116J.

[27] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017, http://dx.doi.org/10.1093/oso/9780198803195.001.0001.

[28] J.C. Bonaldo, S. Mazerat, S. Romero-Baivier, C.L. Martin, Microstructure-based discrete simulations of the compaction of refractory powder composites, Powder Technol. 407 (2022) 117577, http://dx.doi.org/10.1016/j.powtec.2022.117577.

[29] B.L. Trembacki, D.R. Noble, V.E. Brunini, M.E. Ferraro, S.A. Roberts, Mesoscale effective property simulations incorporating conductive binder, J. Electrochem. Soc. 164 (2017) E3613–E3626, http://dx.doi.org/10.1149/2.0601711jes.

[30] B. Derjaguin, V. Muller, Y. Toporov, Effect of contact deformations on the adhesion of particles, J. Colloid Interface Sci. 53 (1975) 314–326, http://dx.doi.org/10.1016/0021-9797(75)90018-1.

[31] F. Radjai, F. Dubois, Discrete-Element Modeling of Granular Materials, ISTE, 2011.
[32] F. Radjaï, I. Preechawuttipong, R. Peyroux, Cohesive granular texture, in: P.A. Vermeer, H.J. Herrmann, S. Luding, W. Ehlers, S. Diebels, E. Ramm (Eds.), Continuous and Discontinuous Modelling of Cohesive-Frictional Materials, Springer Berlin Heidelberg, Berlin, Heidelberg, 2001, pp. 149–162.

[33] M. Pasha, S. Dogbe, C. Hare, A. Hassanpour, M. Ghadi, A linear model of elastoplastic and adhesive contact deformation, Granul. Matter 16 (2014) 152–162, http://dx.doi.org/10.1007/s10035-013-0476-y.

[34] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Engrg. 157 (1998) 69–94.

[35] T. Calvet, J.-M. Vanson, R. Masson, A DEM/FFT approach to simulate the effective thermal conductivity of granular media, Int. J. Therm. Sci. 172 (2022) 107339.

[36] E.N. Primo, M. Touzin, A.A. Franco, Calendering of Li(Ni0.33Mn0.33Co0.33)O2-based cathodes: Analyzing the link between process parameters and electrode properties by advanced statistics, Batter. Supercaps 4 (2021) 834–844, http://dx.doi.org/10.1002/batt.202000324.

[37] E.N. Primo, M. Chouchane, M. Touzin, P. Vazquez, A.A. Franco, Understanding the calendering processability of Li(Ni0.33Mn0.33Co0.33)O2-based cathodes, J. Power Sources 488 (2021) 229361, http://dx.doi.org/10.1016/j.jpowsour.2020.229361.

[38] R. Amin, Y.-M. Chiang, Characterization of electronic and ionic transport in Li1-xNi0.33Mn0.33Co0.33O2 (NMC333) and Li1-xNi0.50Mn0.20Co0.30O2 (NMC523) as a function of li content, J. Electrochem. Soc. 163 (2016) A1512, http://dx.doi.org/10.1149/2.0131608jes.

[39] M.E. Spahr, R. Gilardi, D. Bonacchi, Carbon black for electrically conductive polymer applications, in: R. Rothon (Ed.), Fillers for Polymer Applications, Springer International Publishing, Cham, 2017, pp. 375–400.

[40] J. Puértolas, J. García-García, F. Pascual, J. González-Domínguez, M. Martínez, A. Ansón-Casaos, Dielectric behavior and electrical conductivity of PVDF filled with functionalized single-walled carbon nanotubes, Compos. Sci. Technol. 152 (2017) 263–274, http://dx.doi.org/10.1016/j.compscitech.2017.09.016.

[41] M. Marcinek, J. Syzdek, M. Marczewski, M. Piszcz, L. Niedzicki, M. Kalita, A. Plewa-Marczewska, A. Bitner, P. Wieczorek, T. Trzeciak, M. Kasprzyk, P. Lezak, Z. Zukowska, A. Zalewska, W. Wieczorek, Electrolytes for Li-ion transport – review, Solid State Ion. 276 (2015) 107–126, http://dx.doi.org/10.1016/j.ssi.2015.02.006.

[42] A. Taboada, N. Estrada, F. Radjaï, Additive decomposition of shear strength in cohesive granular media from grain-scale interactions, Phys. Rev. Lett. 97 (2006) 098302.

[43] N. Estrada, E. Azéma, F. Radjai, A. Taboada, Identification of rolling resistance as a shape parameter in sheared granular media, Batteries 84 (2011) 011306, http://dx.doi.org/10.1103/PhysRevE.84.011306, URL: https://link.aps.org/doi/10.1103/PhysRevE.84.011306.