Modeling the effect of calendering pressure on the microstructure and properties of Li-ion battery electrodes

Abstract

The microstructure changes of lithium-ion batteries significantly influence the performance, with calendaring being a crucial manufacturing step linked to the final microstructure. However, the current adjustments of compression parameters are empirical. This work aims to explain the impacts of calendaring force on the microstructure and electrical characteristics by a comprehensive model combining Discrete Element Method (DEM) and Finite Element Analysis (FEA), giving a validated model. A worthy correspondence between simulation and experimental findings was achieved. This work offers a validated combined model, intertwining the microstructure and electrochemical performance, furnishing a numerical simulation foundation for optimizing calendaring parameters for the electrode.

Keywords:

• Calendering
• lithium-ion battery
• electrochemical analysis
• simulations
• discrete element method

1. Introduction

Calendaring, which is a process of compressing electrode material to enhance its density and electrical connectivity, plays a critical role in the manufacture and performance optimization of lithium-ion batteries [1]. Lithium-ion batteries have become a dominant energy storage technology due to the high energy density, high power density, long cycle life, and low self-discharge rate, especially in electric vehicles, portable devices, and renewable energy systems fields [2,3]. The global focus on renewable energy and environmental protection has significantly boosted the development of the electric vehicle market, subsequently increasing the demand for improved battery performance [4,5]. As universally recognized, the macro-properties of materials are inherently influenced by the composition and microstructure [6]. Figure 1 shows that the performance caliber of lithium-ion batteries is intricately intertwined with their constituent components, encompassing electrode materials, electrolytes, and separators.

Figure 1. Schematic diagram of lithium-ion battery.

Variations in the microstructure will lead to great effects on the transport characteristics and electrochemical behavior of electrodes [7,8]. Among the various steps of manufacturing, the calendering process is crucial in controlling the electrode’s microstructure [9]. Early works have proved that the calendering process exerts diverse effects on battery performance [10]. However, too much calendaring pressure can potentially cause diffusion-induced stress (DIS) within batteries [11]. To achieve peak performance, an optimal porosity and corresponding calendering conditions must be found, although this optimization has traditionally been empirical [12].

Using the Discrete Element Method (DEM), we simulate and analyze the motions and forces of electrode particles, focusing on particle contacts, stress, and porosity [13]. This method is applied to study structural changes in electrodes during calendaring [14]. Concurrently, Finite Element Analysis (FEA) examines electrochemical performance changes during calendaring, using a one-dimensional model based on a pseudo-2-dimensional (P2D) battery model [15]. This combined approach will offer an in-depth understanding of how microstructural alterations induced by the calendaring impact electrochemical performance, encompassing rate capability and cycle stability.

In this study, DEM and FEA analyses indicate that calendaring significantly decreases electrode porosity and modifies its internal structure, subsequently affecting electrochemical properties such as capacity, rate capability, and cycle stability. Specifically, an increase in calendaring pressure was observed to lead to a more tortuous electrode structure and alterations in lithium-ion migration pathways. While these changes had minimal effects on electronic conductivity, they exerted notable influences on charge-discharge characteristics and rate performance. These findings, which align with experimental observations, underscore the crucial role of optimizing calendaring conditions to strike a balance between the structural integrity and electrochemical properties of lithium-ion battery electrodes.

2. Materials and methods

Figure 2 illustrates the research workflow and corresponding methodologies employed in this work.

Figure 2. Flow chart of analyzing the impact of calendaring force on the NMC111 electrode.

This study initiates with the establishment of an electrode sample model, justifying the selection of the Edinburgh Elasto-Plastic Adhesion (EEPA) and Bonding models. Subsequently, detailed parameters and steps pertaining to the DEM model setup are elaborated. Furthermore, methodologies employed for analyzing porosity and tortuosity of the microstructure using Python are described. Within the FEA segment, specific parameters, and simulation steps crucial for simulating electrochemical characteristics are documented.

2.1. Sample and DEM

In DEM, the translational and rotational movements of a singular particle within the structure can be accurately characterized as follows: (1) $$\begin{array}{c}{m}^{\mathrm{k}}\frac{\mathrm{d}{\mathit{v}}^{\mathbf{k}}}{\mathrm{d}t}={\mathit{F}}_{\mathrm{b}}^{\mathrm{k}}+{\mathit{F}}_{\mathrm{p}}^{\mathrm{k}}\end{array}$$(1) (2) $$\begin{array}{c}{I}^{\mathrm{k}}\frac{\mathrm{d}{\mathit{w}}^{\mathbf{k}}}{\mathrm{d}t}={\mathit{M}}_{\mathrm{b}}^{\mathrm{k}}+{\mathit{M}}_{\mathrm{P}}^{\mathrm{k}}\end{array}$$(2)

$I^{\mathrm{k}}$ and $m^{\mathrm{k}}$ refer to the moment of inertia and mass of the $\mathit{\text{kth}}$ particle, respectively. The variables ${\mathit{v}}^{\mathbf{k}}$ and ωk represent the particle’s translational and angular velocities. The bond torque and bond force, symbolized as ${\mathit{M}}_{\mathrm{b}}^{\mathrm{k}}$ and ${\mathit{F}}_{\mathrm{b}}^{\mathrm{k}},$ are used to describe the interactions in the CBD phase. The contact force with neighboring particles is denoted as ${\mathit{F}}_{\mathrm{p}}^{\mathrm{k}},$ and ${\mathit{M}}_{\mathrm{P}}^{\mathrm{k}}$ signifies the contact torque resulting from the tangential rolling friction and contact force among adjacent particles. To simplify the following sections, the superscript $\mathrm{k}$ specifying an individual particle is omitted.

2.1.1. Base contact model

This study employs a simulation methodology to investigate deeply into the calendering process. The EEPA model serves to explain the interactions among AM particles [16]. The contact normal force is expressed as the cumulative effect of the hysteretic spring force and the damping force: (3) $$\begin{array}{c}{\mathit{F}}_{\mathrm{n},\mathrm{p}}=\left(f_{\mathrm{\text{hys}}}+f_{\mathrm{\text{nd}}}\right)n\end{array}$$(3) where $\mathit{n}$ represents the unit normal vector extending from the contact point directly to the center of the particle. The relationship between the force and displacement characteristics of $f_{\mathrm{\text{hys}}}$ can be formulated by: (4) $$\begin{array}{c}{f}_{\mathrm{\text{hys}}}=\{\begin{array}{l}{f}_0+k_1{\delta }^n\mathit{\text{ if }}{k}_2\left({\delta }^n-{\delta }_p^n\right)\ge k_1{\delta }^n\\ f_0+k_2\left({\delta }^n-{\delta }_p^n\right)\mathit{\text{if }}{\mathrm{ }k}_1{\delta }^n>k_2\left({\delta }^n-{\delta }_p^n\right)>-k_{\mathrm{\text{adh}}}{\delta }^n\\ f_0-k_{\mathrm{\text{adh}}}{\delta }^n\mathit{\text{ if }}{k}_{\mathrm{\text{adh}}}{\delta }^n\ge k_2\left({\delta }^n-{\delta }_p^n\right)\end{array}\end{array}$$(4) where $f_0$ is constant adhesion force, an attractive force that is always present between particles during their contact and separation ${\delta }^n$ represents the total amount of normal overlap, and ${\delta }_p^n$ is the total overlap resulting from plastic deformation. ${\mathrm{ }k}_1$ is the loading stiffness, $k_2$ is the unloading stiffness, and $k_{\mathrm{\text{adh}}}$ is the adhesive stiffness, $n$ is the non-linear index parameter. The schematic representations of particle contact along with the normal force-overlap curve pertinent to this model are illustrated in Figure 3.

Figure 3. Normal contact force-displacement function for the EEPA contact model [17].

The loading stiffness $k_1$ is calculated by: (5) $$\begin{array}{c}{k}_1=\frac{4}{3}\sqrt{R^{\mathrm{*}}}{E}^{\mathrm{*}}\end{array}$$(5) where $E^{\mathrm{*}}$ and $R^{\mathrm{*}}$ refer to the equivalent Young’s modulus and radius. $E^{\mathrm{*}}$ is a representative 'averaged’ or 'composite’ value of the Young’s moduli of two different materials in contact. It essentially represents the overall stiffness of the contact interface between the two particles. In particle interactions, the contact region is not a point but an area. This area can be approximated as a circle with a radius referred to as the contact radius. When particles of varying shapes or sizes come into contact, the Equivalent Contact Radius provides a suitable descriptor of this contact radius [15].

The contact plasticity ratio, ${\lambda }_p$ is calculated as: (6) $$\begin{array}{c}{\lambda }_p=\left(1-\frac{k_1}{k_2}\right)\end{array}$$(6)

The damping force $f_{\mathrm{\text{nd}}}$ can be calculated by: (7) $$\begin{array}{c}{f}_{\mathrm{\text{nd}}}=-{\beta }_n{v}_n\end{array}$$(7)

$f_{\mathrm{\text{nd}}}$ is the resistive force arising from the relative motion between particles, serving to dissipate the associated energy and reduce relative movement. It is directly proportional to the relative velocity ($v_n$) and modulated by the damping coefficient (${\beta }_n$). In simulations, the damping force facilitates realistic particle interactions by preventing excessive motion or vibration, where ${\beta }_n$ is the normal damping coefficient: (8) $$\begin{array}{c}{\beta }_n=\sqrt{\frac{4m^{\mathrm{*}}{k}_1}{1+{\left(\frac{\pi }{\mathrm{\text{ln}}\hspace{0.17em}e}\right)}^2}}\end{array}$$(8) where the coefficient of restitution is denoted by e.

The contact tangential force ${\mathit{F}}_{\mathrm{t},\mathrm{p}}$ is determined as the aggregate of tangential spring force, ${\mathit{f}}_{\mathrm{\text{ts}}}$ and tangential damping force, ${\mathit{f}}_{\mathrm{\text{td}}}$: (9) $$\begin{array}{c}{\mathit{F}}_{\mathrm{t},\mathrm{p}}=\left({\mathit{f}}_{\mathrm{\text{ts}}}+{\mathit{f}}_{\mathrm{\text{td}}}\right)\end{array}$$(9)

The tangential spring force is formulated as follows: (10) $$\begin{array}{c}{\mathit{f}}_{\mathrm{\text{ts}}}=\left({\mathit{f}}_{\mathrm{\text{ts}}\left(\mathrm{n}-1\right)}+\Delta {\mathit{f}}_{\mathrm{\text{ts}}}\right)\end{array}$$(10) where the force of the previous time step is represented by ${\mathit{f}}_{\mathrm{\text{ts}}(n-1)},$ while the incremental tangential force $\Delta {\mathit{f}}_{\mathrm{\text{ts}}}$ is determined using the following calculation: (11) $$\begin{array}{c}{\Delta \mathit{f}}_{\mathrm{\text{ts}}}=-k_t{\mathit{v}}_{\mathit{t}}\Delta t\end{array}$$(11) where the tangential stiffness, designated as $k_t,$ represents the resistance to deformation in the tangential direction. The velocity in this tangential direction is designated as ${\mathit{v}}_{\mathit{t}},$ while $\Delta t$ refers to the duration of each time step. The tangential stiffness can be determined through the application of the subsequent equation: (12) $$\begin{array}{c}{k}_t={\zeta }_{\mathrm{\text{tm}}}\{\begin{array}{c}{k}_1\mathit{\text{ if n}}=1\\ 8G^{\mathrm{*}}\sqrt{R^{\mathrm{*}}{\delta }_n}\mathit{\text{ if n}}>1\end{array}\end{array}$$(12) where the tangential damping force, denoted as ${\mathit{f}}_{\mathrm{\text{td}}},$ is calculated using the tangential stiffness multiplier ${\zeta }_{\mathrm{\text{tm}}}$ and the equivalent shear modulus $G^{\mathrm{*}}$. (13) $$\begin{array}{c}{\mathit{f}}_{\mathrm{\text{td}}}=-{\beta }_t{\mathit{v}}_{\mathit{t}}\end{array}$$(13)

The tangential damping coefficient ${\beta }_t$ is determined through the following calculation: (14) $$\begin{array}{c}{\beta }_t=\sqrt{\frac{4m^{\mathrm{*}}{k}_t}{1+{\left(\frac{\pi }{\mathrm{\text{ln}}\hspace{0.17em}e}\right)}^2}}\end{array}$$(14)

The critical tangential force value, designated as $f_{\mathrm{\text{ct}}},$ is established to be equivalent to: (15) $$\begin{array}{c}{f}_{\mathrm{\text{ct}}}\le \mu \left|f_{\mathrm{\text{hys}}}+k_{\mathrm{\text{adh}}}{\delta }^n-f_0\right|\end{array}$$(15) where $\mu$ is the friction parameter.

2.1.2. The binder phase interactions

The Bonding model serves to explain the mechanical responses exhibited by the inter-particle binder phase [18]. This model establishes bonds among AM particles, aiming to capture the mechanical behavior of the CBD within electrode structures. These bonds are assumed to arise when the center-to-center distance between two neighboring particles, designated as PA and PB, satisfies a specific condition [19]: (16) $$\begin{array}{c}\mid \mid \mathit{\text{PA}}-\mathit{\text{PB}}\mid \mid \le \zeta b\cdot \left(\mathit{\text{RA}}+\mathit{\text{RB}}\right)\end{array}$$(16) where the bond radius is determined by considering the particle radius, $\mathit{\text{RA}}$ and $\mathit{\text{RB}}.$ (17) $$\begin{array}{c}\mathit{\text{Rb}}=\lambda b\cdot \mathrm{\text{min}}\hspace{0.17em}\left(\mathit{\text{RA}},\mathit{\text{RB}}\right)\end{array}$$(17)

Two multipliers were established as $\zeta b$=1.67, $\lambda b$=0.5 [19]. To ensure the aggregation of all particles without any unbound ones, these parameters were specifically selected. Before compression, inter-particle bonds were formed among the particles, indicating that these multipliers were exclusively used for creating the initial bonded assembly and were inactive during the compression process. During compression, the inter-particle bonds exhibited deformation or breakage, closely aligning with the actual mechanical properties of the CBD phase [3,4].

The normal force ${\mathit{F}}_{\mathrm{n},\mathrm{b}}$ can be represented as the cumulative bond force from the preceding time step ${\mathit{F}}_{\mathrm{n}\left(n-1\right),\mathrm{b}}$ (18) $$\begin{array}{c}{\mathit{F}}_{\mathrm{n},\mathrm{b}}=\left({\mathit{F}}_{\mathrm{n}\left(\mathrm{n}-1\right),\mathrm{b}}+-k_{\mathrm{n},\mathrm{b}}{\mathit{v}}_{\mathbf{n},\mathbf{b}}\Delta t\right)\end{array}$$(18) where $\mathrm{ }{k}_{\mathrm{n},\mathrm{b}}$ denotes the normal stiffness of the bond, and ${\mathit{v}}_{\mathbf{n},\mathbf{b}}$ signifies the velocity in the normal direction.

The tangential spring force, designated as ${\mathit{F}}_{\mathrm{t},\mathrm{b}},$ is calculated as follows: (19) $$\begin{array}{c}{\mathit{F}}_{\mathrm{t},\mathrm{b}}=\left({\mathit{F}}_{\mathrm{t}\left(\mathrm{n}-1\right),\mathrm{b}}+-k_{\mathrm{t},\mathrm{b}}{\mathit{v}}_{\mathbf{t},\mathbf{b}}\Delta t\right)\end{array}$$(19) where the bond tangential force at the preceding time step is denoted as ${\mathit{F}}_{\mathrm{t}\left(n-1\right),\mathrm{b}},$ while $k_{\mathrm{t},\mathrm{b}}$ represents the tangential stiffness of the bond, and ${\mathit{v}}_{\mathbf{t},\mathbf{b}}$ corresponds to the tangential velocity.

The bond tangential moments ${\mathit{M}}_{\mathrm{t},\mathrm{b}}$ and bond normal moments ${\mathit{M}}_{\mathrm{n},\mathrm{b}}$ can be formulated by considering the incremental values obtained from the preceding time step, expressed as follows: (20) $$\begin{array}{c}{\mathit{M}}_{\mathrm{t},\mathrm{b}}=\left({\mathit{M}}_{\mathrm{t}\left(\mathrm{n}-1\right),\mathrm{b}}+\Delta {\mathit{M}}_{\mathrm{t},\mathrm{b}}\right)\end{array}$$(20) (21) $$\begin{array}{c}{\Delta \mathit{M}}_{\mathrm{t},\mathrm{b}}=-k_{\mathrm{t},\mathrm{b}}{A_b\mathit{\omega }}_{\mathbf{t},\mathbf{b}}\Delta t{I}_b \end{array}$$(21) (22) $$\begin{array}{c}{\mathit{M}}_{\mathrm{n},\mathrm{b}}=\left({\mathit{M}}_{\mathrm{n}\left(\mathrm{n}-1\right),\mathrm{b}}+\Delta {\mathit{M}}_{\mathrm{n},\mathrm{b}}\right)\end{array}$$(22) (23) $$\begin{array}{c}{\Delta \mathit{M}}_{\mathrm{n},\mathrm{b}}=-k_{\mathrm{n},\mathrm{b}}{A_b\mathit{\omega }}_{\mathbf{n},\mathbf{b}}\Delta t{J}_b\end{array}$$(23) where the cross-sectional area of the bond is denoted as $A_b,$ while $I_b$ represents the moment of inertia and $J_b$ stands for the polar moment of inertia of the bond. The relative angular velocities in the tangential and normal directions are indicated by ${\mathit{\omega }}_{\mathbf{t},\mathbf{b}}$ and ${\mathit{\omega }}_{\mathbf{n},\mathbf{b}},$ respectively.

The bond strength is determined utilizing the beam theory, formulated as follows: (24) $$\begin{array}{c}-\frac{{\mathit{F}}_{\mathrm{n},\mathrm{b}}}{A_{\mathrm{b}}}+\frac{\left|{\mathit{M}}_{\mathrm{t},\mathrm{b}}\right|R_{\mathrm{b}}}{I_{\mathrm{b}}}<{\sigma }_{\mathrm{n}}\end{array}$$(24) (25) $$\begin{array}{c}\frac{\left|{\mathit{F}}_{\mathrm{t},\mathrm{b}}\right|}{A_{\mathrm{b}}}+\frac{\left|{\mathit{M}}_{\mathrm{n},\mathrm{b}}\right|R_b}{J_{\mathrm{b}}}<{\sigma }_{\mathrm{t}}\end{array}$$(25) where the tensile strength of the bond is represented by ${\sigma }_{\mathrm{n}},$ while ${\sigma }_{\mathrm{t}}$ denotes its shear strength.

2.1.3. Structure generation and simulation setup

In this investigation, Altair EDEM 2021 software was employed to conduct the DEM simulation. Due to the significant sphericity exhibited by the NMC particles, they were effectively modeled as spherical entities. The subdomain dimensions were precisely set to 150 μm × 150 μm × 75 μm. The system setup comprises a steel top plate simulating the calendering roller, alongside an aluminum bottom plate replicating the function of the current collector. Periodic boundary conditions were imposed on the lateral sides of the simulation environment, ensuring a seamless and representative spatial framework. A comprehensive listing of the specific DEM model parameters can be found in Table 1.

Table 1. DEM Simulation parameters used.

Description (Unit)	Parameter	Value	Source
Particle Young’s modulus (GPa)	EPar	144	[5,6]
Particle density (kg/m^3)	ρPar	4750
Particle Poisson’s ratio (–)	υPar	0.2
Particle-particle static friction (–)	μp,s	0.25
Particle-particle rolling friction (–)	μp,r	0.01
Coefficient of restitution (–)	e	0.2
Plate Young’s modulus (Steel) (GPa)	ESt	208	[20]
Plate density (Steel) (kg/m³)	ρSt	7900
Plate Poission’s ratio (Steel) (–)	υSt	0.3
Plate Young’s modulus (Aluminum) (GPa)	EAl	68.9
Plate density (Aluminum) (kg/m³)	ρAl	2700
Plate Poission’s ratio (Aluminum) (–)	υAl	0.3
Speed of calendering (m/s)		0.01
Bond stiffness (N/m^3)	kn,b	1.50e + 14
Bond strength (Pa)	σn	3.00e + 08
Cross-section (μm)		150 × 150
Plasticity ratio	λp	0.5	[21]
Surface energy (J/m²)	Δγ	5
Pull-off force (N)	f0	0
Tensile exponent		2.5
Tangential stiffness multiplier	ζtm	0.4

2.2. Microstructure analysis

Before and after the calendering process, the microstructure of the electrode was analyzed through the calculation of porosity and tortuosity. These changes in properties have the potential to significantly influence the battery’s performance during the calendering stage.

2.2.1. Porosity analysis

In DEM simulations, the use of a soft particle model unavoidably leads to deformations and overlaps among AM particles during the calendering process [22]. This complexity prevents the direct application of analytical methods for accurate porosity determination. To address this, a slicing technique was adopted in this study. The accuracy of this analysis method is constrained by the resolution of the images.

2.2.2. Tortuosity analysis

Tortuosity in the electrolyte phase signifies the lithium-ion transport potential, which is crucial for the charge-discharge rate of the battery electrode. An image processing approach was employed to quantify tortuosity, showcasing the impact of the calendering process on the battery microstructure. This metric was implemented from an open-source MATLAB code [23], simplified in Python.

2.3. Electrochemical characterization and FEA

2.3.1. Finite element implementation

The microscopic information of the electrodes at four different porosities was directly extracted from the ‘Result Data’ feature in EDEM 2021, and subsequently integrated into COMSOL Multiphysics 6.1 for one-dimensional model development based on a P2D physio chemical battery model framework. To compare with experimental results, battery components not generated in DEM, including CBD, anode, and electrolyte were modeled based on experimental conditions studying the effects of calendaring [24]. In COMSOL, grid division and generation were carried out directly. Discretized transport and electrode kinetic equations were addressed using an implicit solver, leveraging the backward differentiation formula for temporal integration. The solving process also employed the same formula for time stepping, with an initial size set at 0.001s.

2.3.2. Electrochemical characterization

In this work, we simulated the construction of lithium-ion batteries, utilizing the mean particle diameter and porosity of the NMC111 cathode sourced from EDEM. The remaining simulation parameters are detailed in Table 2. In this model construct, lithium foil was designated as the counter electrode, Celgard 2400 served as the separator, and a 1 mol/dm³ solution of LiPF6 in EC + DEC (1:1) was employed as the electrolyte. The true cross-sectional area of the half-cell was referenced from previous experimental works, set as a circle with a diameter of 12.7 mm [27]. All experimental protocols were conducted at a consistent temperature of 303.15K. Initially, under the stipulated voltage constraints of 4.5 V for charging and 3.0 V for discharging, the first charge-discharge cycle across varying cathode porosities was meticulously examined at a charge and discharge rate of 0.1 C. This was followed by an extensive assessment involving 25 cycles at C/10 charge and discharge rates for cathodes of differentiated porosities. After the formation cycles, the rate performances for cathodes of distinct porosities were evaluated, with the discharge capacity, as derived from the final formation cycle, proving instrumental in determining the C-rate for each respective cell. The rate performance assessments encompassed exhaustive discharges at rates of C/10, C/5, C/2, 1, 2, 5, and 10 C down to a terminal voltage of 3.0 V, each invariably succeeded by a C/10 charge up to 4.5 V [28].

Table 2. Electrochemical simulation parameters used.

Description (Unit)	Parameter	Value	Source
Positive electrode thickness (μm)	δPos	160	[25]
Separator thickness (μm)	δSep	24
Volume fraction of active material	εam	0.57
Volume fraction of inert material	εim	0.12
Porosity of separator ive electrode	εsep	0.5
Maximum solid phase concentration (mol/m^3)	cp,max	48897
Initial solid phase concentration (mol/m^3)	cs0	26480	[26]
Initial electrolyte concentration (mol/m^3)	ce0	1000
Radius of the particle (μm)	R0	5.36
Transfer coefficient for anodic current	aa	0.5
Transfer coefficient for cathodic current	ac	0.5
Gas constant (J/(mol.K))		8.314
Li^+ effective diffusion coefficient of NMC (m^2/s)	Dnmc	1.5e-13
Geometrical factor in solid phase	γ	5
Li^+ diffusion coefficient in electrolyte (m^2/s)	De	1.5e-13
Li^+ transference number	t+	0.38

In summary, this chapter thoroughly elaborates on all the simulation and modeling steps involved in this study. Initially, it introduces the establishment steps of the DEM model, including the rationale behind the selection of the EEPA and Bonding contact models, along with the detailed parameters and steps for setting up the DEM model. Concurrently, it also discusses the characterization methods for microstructural aspects such as porosity and tortuosity. Following this, the FEA section covers the specific parameters and characterization steps for simulating electrochemical properties.

3. Results and discussion

In this section, a comprehensive validation of the integrated model combining DEM and FEA electrochemical evaluations is conducted. The validation focuses on three primary dimensions: macro-mechanical properties during calendaring, microstructural details, and electrochemical characteristics. Initially, a comparison is made between the simulated and experimental displacement and pressure curves, followed by an assessment of porosity distribution and pathway tortuosity. To compare with experimental data from other studies, the dynamic calendaring process is simplified to represent electrodes at four distinct calendaring levels. These levels are characterized by four porosities: 0.6, 0.58, 0.55, and 0.37. The rationale for selecting these values stems from the chosen comparative experimental data results [19]. After the inclusion of the CBD phase, an in-depth analysis is facilitated through electrochemical modeling to decipher the changes in battery behaviors. The aim is to thoroughly validate this integrated model, establishing a connection between calendaring dynamics, nuances of microstructure, and its electrochemical properties.

3.1. DEM model

The results of the establishment of the DEM model are presented first. Initially, the particles in the uncalendered structure of NMC111 were studied, with their particle size distribution illustrated in Figure 4.

Figure 4. Particle size distribution in DEM simulation.

To simulate the mechanical response of the CBD phase within the electrode structures, inter-particle bonds were formed among the AM particles, as shown in Figure 5b. The conditions for bond formation and the thickness of the bond torque are specified in Equations (16)(16) $$\begin{array}{c}\mid \mid \mathit{\text{PA}}-\mathit{\text{PB}}\mid \mid \le \zeta b\cdot \left(\mathit{\text{RA}}+\mathit{\text{RB}}\right)\end{array}$$(16) and (17)(17) $$\begin{array}{c}\mathit{\text{Rb}}=\lambda b\cdot \mathrm{\text{min}}\hspace{0.17em}\left(\mathit{\text{RA}},\mathit{\text{RB}}\right)\end{array}$$(17) . To produce a representative simulation of AM contacts, AM particles were first generated within an extended subdomain using the 'Particle Factory’ function of EDEM under 0 gravity conditions. This ensured no contacts or relative movement in the initial state. Subsequently, bonding between particles was established based on the given strength and contact radius. This approach allows for a more accurate representation of AM contacts, preventing localized concentrations of AM particles. The calendaring process was simulated by lowering the top plate at a constant velocity of 0.01 m/s by ≈75μm, leading to the gradual compression of the electrode structure, as depicted in Figure 5c. Once calendared to the target height (≈ 35 μm), the entire calendaring process simulation concluded, as shown in Figure 5d. The time step was set to 10⁻¹¹ s.

Figure 5. DEM simulation process: (a) constructed structure in DEM; (b) structure with inter particle bonds; (c) DEM simulation setup; and (d) calendering simulation completed.

3.2. Microstructure analysis

3.2.1. Displacement and pressure analysis

The displacement and pressure during the compaction process is further analyzed through DEM. As shown in Figure 6, a discrepancy emerged from the 26um to 38um phase, which can be attributed to the mechanical response of the whole electrode pressure change gradually shifting from being dominated by Bonding fractures to being dominated by interparticle friction. In the experiment, the particle size distribution is more uneven, and the difference in diameter size is larger, which will cause variations in the interparticle forces [29].

Figure 6. Comparison of normal displacement under different compression pressures, the experimental data is from Ge et al. [19].

3.2.2. Porosity analysis

Figure 7 illustrates the images obtained from the slicing procedure, which are subsequently converted into binary images.

Figure 7. Processing procedure of the cross-section image of the AM phase.

In this investigation, the resolution of the cross-section is maintained at a pixel length of 0.5 μm to ensure the reliability of the results [19]. The overlapping regions are intuitively displayed in the cross-sectional images. In the following analysis, the porosity of pore and binder phases, denoted as $\mathit{\epsilon },$ is represented as: (26) $$\begin{array}{c}\epsilon =1-{\mathit{\phi }}_{\mathit{\text{Par}}}\end{array}$$(26) where ${\mathit{\phi }}_{\mathbf{\text{Par}}}$ is defined as the volume fraction of AM particles.

The porosity variations observed in the DEM simulations align closely with the corresponding experimental results, as exposed in Figure 8, exhibiting a nearly identical trend.

Figure 8. Comparison of porosity ε under different compression pressures, the tomography data is from Ebner et al. [30].

A noticeable discrepancy can be discerned during the initial phase, which can likely be attributed to two primary factors. Firstly, the particles used in the experiments are anisotropic imperfect spherical particles with a more uneven particle size distribution; this can lead to a loss in particle volume. Secondly, the mechanical response between CBD in the experiments is more complex than the Bonding model used in the simulations, which can result in different pressures being required to achieve the same porosity.

3.2.3. Tortuosity analysis

Figure 9 exhibits the process of quantifying tortuosity from a slice, which is converted to grayscale, binarized using Otsu’s method, and skeletonized.

Figure 9. Processing procedure of the cross-section image of the AM phase.

A breadth-first search (BFS) is utilized to find the shortest path on the skeleton, which is distinctly highlighted in red. This metric dynamically demonstrates the effect of the calendering process on the battery microstructure. In this work, we computed the average tortuosity factor $\mathit{\tau }$_Pore+binder for the porous phase that includes both pores and binder [19]. The tortuosity factor $\mathit{\tau }$ can be defined by the following formula [23]: (27) $$\begin{array}{c}\mathit{\tau }=\frac{\mathit{\text{Le}}}{\mathit{L0}\mathrm{ }} \end{array}$$(27) where $\mathit{\text{Le}}$ represents the actual path length, and $\mathit{L}0$ represents the straight-line length.

As illustrated in Figure 10, the changes in tortuosity from DEM simulations align closely with the corresponding experimental results.

Figure 10. Comparison of tortuosity $\mathrm{\tau }$ under different compression pressures, the tomography data is from Ge et al. [19].

Minor differences can be observed at the four distinct pressure stages, which could likely be due to two main factors. Firstly, the particle size distribution in the experiments is more uneven, which could lead to a larger tortuosity in the experimental results. Secondly, the tortuosity calculation method is constrained by the pixel size of the two-dimensional cross-section, which could result in minor differences in the tortuosity results.

3.3. Electrochemical analysis

3.3.1. Electronic conductivity

Figure 11 presents a comparison of simulated and experimental results for the electronic conductivity of the electrode.

Figure 11. Comparison of electronic conductivity, (a) simulation data; (b) experimental data is from Zheng et al. [31].

In the simulation, the electronic conductivity of the NMC111 electrode layer at diverse porosities is about 0.29 S/cm, while the experimental result is around 0.30 S/cm. Moreover, this value does not change significantly either in the simulation or in the experimental process. Hence, it’s concluded that the calendaring process does not significantly affect the electronic conductivity. The electronic conductivity is related to the type of material, experimental conditions, battery dimensions, and microstructure. The electrode materials and conditions in the simulation are consistent well with the compared experimental conditions.

3.3.2. First charge–discharge profiles

Figure 12 compares the simulation and experimental results of the first charge-discharge cycle.

Figure 12. Comparison of first charge–discharge profiles (0.1 C 303.15k), (a) simulation data; (b) experimental data is from Zheng et al. [31].

To better comprehend these disparities, the First Coulombic Efficiency (FCE) was introduced, serving as a vital criterion for assessing electrode quality. During electrochemical cycling, a high electrochemical efficiency translates to reduced lithium consumption.

Within the porosity range of 0.6-0.37, both simulation and experiment yield a reversible capacity of about 180mAh/g, indicating that compaction does not bring about noticeable changes to the reversible capacity. In the charging simulation process, the half-cell with a porosity of 0.37 exhibits a smaller specific capacity difference, which might be due to the poorer wettability between the electrode and the electrolyte at a lower porosity. The diminished porosity might curtail the rapid migration paths of lithium ions, leading to a reduction in specific capacity. The enhanced compactness of the electrode could detrimentally affect the wettability of the electrolyte. When the electrolyte doesn’t sufficiently wet the electrode material, it further hinders the transport of lithium ions, especially inside the active material. Such wettability issues are likely more pronounced in electrodes with the least porosity. Another point of consideration is the utilization of the active material [32]. Electrodes subjected to rigorous calendaring may present regions of the active material that aren’t adequately in contact with the electrolyte. Such regions might remain electrochemically inactive during cycling, thereby diminishing the overall specific capacity. (28) $$\begin{array}{c}\mathit{\text{FCE}}=\frac{\mathbf{\text{First}}\mathrm{ }\mathbf{\text{Charge}}\mathrm{ }\mathbf{\text{Capacity}}\left(\frac{\mathbf{\text{mAh}}}{\mathbf{g}}\right)}{\mathbf{\text{First}}\mathrm{ }\mathbf{\text{Discharge}}\mathrm{ }\mathbf{\text{Capacity}}\left(\frac{\mathbf{\text{mAh}}}{\mathbf{g}}\right)} \end{array}$$(28) where the FCE is defined as the ratio of the first discharge capacity to the first charge capacity.

As depicted in Figure 13, an inverse correlation is observed between porosity and the FCE. This decrement in FCE with increasing porosity is presumably attributed to the decrease in the electrode’s specific surface area. For numerous lithium storage materials, it has been established that, under certain conditions, the magnitude of SEI formation reactions at the electrode/electrolyte interface scales directly with the electrode’s specific surface area [33]. Consequently, the reduced SEI formation reactions resulting from the decreased specific surface area caused by calendering explain the enhancement in FCE observed at lower porosities.

Figure 13. The first coulombic efficiency of the Li [Ni1/3Mn1/3Co1/3] O₂ cathode at different porosities.

Here we see a compromise for a high-quality electrode laminate, both a high specific capacity and high FCE are preferable. Thus, calendering process must be accurately regulated to achieve the ideal electrode porosity.

3.3.3. Cycling performance

Figure 14 compares the cycling performance results from simulation and experiment.

Figure 14. Comparison of the cycling performance, (a) simulation data; (b) experimental data is from Zheng et al. [31].

Both the experiment and simulation exhibit good cycling performance, with a capacity loss of less than 5% within 0.37-0.6 porosities. The cycling performance is believed to be related to the electrode volume changes associated with lithium-ion insertion and extraction [34]. The volume change of the NMC111 electrode during the charge-discharge process is about 3%. The same electrode material leads to similar cycling performance, resulting in similar comparative results.

3.3.4. Rate capability

Figure 15 displays the trend of rate performance variations between simulation and experiment. The rate performance obtained from simulation decays faster than the experimental results at high rates beyond 1 C, which might be associated with the differences between the utilized P2D model and the real experiment.

Figure 15. Comparison of the rate capability, (a) simulation data; (b) experimental data is from Zheng et al. [31].

In the P2D model, the particle size of active materials is simplified to be uniform, and the tortuosity effect brought about by a higher degree of calendering is smaller; under high-rate operation conditions, this discrepancy is amplified as lithium ions have higher transmission efficiency in media with lesser tortuosity [35]. Therefore, there’s more capacity decay under high rates. Additionally, the P2D model only considers lithium-ion kinetics motion in the X direction, this difference is also amplified at high rates, thereby contributing to the variation in decay speed [36].

3.4. Summary

In this section, a comprehensive validation of the integrated model, which combines DEM and FEA electrochemical evaluations, is conducted. The validation primarily focuses on three main dimensions: macro-mechanical properties during calendaring, microstructural details, and electrochemical characteristics. Initially, a comparison is made between simulated and experimental displacement and pressure curves, followed by an assessment of porosity distribution and pathway tortuosity. Subsequently, an in-depth analysis is facilitated through electrochemical modeling to decipher the changes in battery behaviors. The aim is to thoroughly validate this integrated model, establishing a connection between calendaring dynamics, nuances of microstructure, and its electrochemical properties. Through comparing experimental and simulated data, the accuracy and reliability of the model in mechanical properties, microstructure, and electrochemical characteristics are demonstrated, providing a reliable integrated model for the mutual influence between microstructure and electrochemical performance during the electrode calendaring process. At the same time, this comparison also reveals some discrepancies between simulation and experimental results, which might be caused by experimental conditions, material properties, and model simplifications.

4. Conclusions

In this work, the DEM combined with a bonded particle model was employed to investigate the microstructural evolution of the NMC111 lithium-ion cathode during calendaring. The results demonstrated a pronounced decline in electrode porosity throughout the calendaring process. Specifically, the porosity reduced from an initial value of 0.6 and reached 0.37 upon the application of a calendaring pressure of 200 MPa, highlighting the increasing compactness of the electrode structure. Moreover, as the calendaring pressure was elevated from 60 MPa to 200 MPa, the tortuosity of the internal pathways in the electrode distinctly rose from 1.45 to 1.65. Such a trend is likely associated with an augmentation in inter-particle friction within the electrode, especially at elevated calendaring pressures. Furthermore, when comparing the idealized DEM structures devoid of the CBD phase with experimental characterizations, the DEM microstructure changes during calendaring matched accurately, with initial discrepancies in porosity attributed to a more intricate particle size distribution in real electrodes.

Electrochemical analyses on the half-cell featuring NMC111 as the cathode were conducted utilizing the FEA approach. The initial charge-discharge characteristics of the electrode manifested some alterations during calendaring. Notably, when the electrode’s porosity stood at 0.37, its reversible capacity slightly diminished compared to other electrodes, but the first coulombic efficiency saw an uptick, mirroring the trends observed in experimental findings. Delving deeper, the diminished reversible capacity can be attributed to decreased electrolyte wettability and enhanced tortuosity in lithium-ion migration pathways within the electrode. In terms of cycling performance, electrodes across various porosities registered a capacity loss of less than 5% post 25 electrochemical cycles, indicating a minor impact of the calendaring process on cycling performance. Conversely, in the realm of rate capability, electrodes subjected to higher calendaring exhibited diminished performance at elevated charge-discharge rates, emphasizing the pivotal role of electrolyte wettability and internal lithium-ion transport at higher rates. Such observations might also be linked to the inherent simplifications of the P2D model utilized in the simulations. The P2D model predominantly focuses on the kinetic behavior of lithium ions in the X-direction, whereas in real electrode architectures, the migration paths of lithium ions might be more convoluted. Additionally, the P2D model’s assumption of a uniform particle size distribution might have also impacted the migration efficiency of lithium ions, particularly at higher rates.

The analysis established a direct connection between calendaring, microstructure, and electrochemical properties of the electrode, providing a numerical basis for the design and optimization of calendaring parameters for the NMC111 electrode. Foreseeable extensions of this study will focus on parameters such as calendaring velocity, thermal dynamics, and exerted pressures. Within the purview of this model, spherical NMC entities were harnessed, given their pervasive presence in the cathodes of lithium-ion batteries. For electrode materials like graphite, showcasing non-spherical or anomalously contoured entities, the versatility of the DEM framework permits depictions through a multi-sphere configuration, underscoring the model’s adaptability.

Data availability statement

The data used to support the findings of this study are included within the article.

Disclosure statement

The authors report there are no competing interests to declare.

Additional information

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Notes on contributors

Junye Pan

Junye Pan, UNSW 2023/11 Master graduation. His research interests include DEM Simulation,

Yudong Zou

Yudong Zou is now a PhD at School of Materials Science and Engineering. His research interests include DEM simulation.