# A Level Set Discrete Element Model (LS-DEM) for sintering with an optimization-based contact detection

Brayan Paredes-Goyes, David Jauffres, Christophe L Martin

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

## To cite this version:

Brayan Paredes-Goyes, David Jauffres, Christophe L Martin. A Level Set Discrete Element Model (LS-DEM) for sintering with an optimization-based contact detection. 2024. hal-04594920

HAL Id: hal-04594920

https://hal.science/hal-04594920

Preprint submitted on 30 May 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Level Set Discrete Element Model for sintering with an optimization-based contact detection

Brayan Paredes-Goyes, David Jauffres, Christophe L Martin

► To cite this version:

Brayan Paredes-Goyes, David Jauffres, Christophe L Martin. A Level Set Discrete Element Model for sintering with an optimization-based contact detection. Journal of Theoretical, Computational and Applied Mechanics, 2025, 10.46298/jtcam.13721. hal-04594920v2

HAL Id: hal-04594920

https://hal.science/hal-04594920v2

Submitted on 24 Apr 2025

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
# A Level Set Discrete Element Model (LS-DEM) for sintering with an optimization-based contact detection

Brayan PAREDES-GOYES<fnref n="1" />,<fnref n="2" />, David JAUFFRES<fnref n="1" />, and Christophe L. MARTIN<fnref n="1" />

Sintering is a high temperature process for the consolidation of ceramic, metal and polymer powders. The Discrete Element Method (DEM) has been effectively used to model the sintering process at the particle scale considering spherical particles. However, standard manufacturing processes rarely deal with spherical particles. As sintering is a curvature-controlled process, it is important to take into account the deviation from sphericity. This study presents a DEM sintering model for non-spherical particles. The description and dynamic evolution of arbitrary shape particles is achieved by using the Level Set discrete element method (LS-DEM). The original LS-DEM approach uses boundary nodes on the particles to detect contacts. We employ an optimization-based contact detection approach. This improves the capture of small contacts, which is important for a correct description of sintering evolution with reasonable CPU-time consumption. A Newton-Raphson scheme is employed for the optimization algorithm. The normal force and neck size evolution expressions of spherical particles are adapted for arbitrary shape particles by using the local curvature at the contact. The developed model is validated for elastic contacts on superquadric ellipsoids. It is compared with standard DEM on spheres for sintering. The model is applied to investigate the consolidation kinetics of a packing of ellipsoidal particles. It is shown, that a deviation from sphericity is beneficial for both prolate and oblate ellipsoids. An optimum aspect ratio is evaluated, demonstrating that particles that are too elongated slow down densification kinetics.

Keywords discrete element method, level set, sintering, non-spherical particles, contact detection

## 1 Introduction

Sintering is a prominent high-temperature process to manufacture ceramic, metallic and polymeric materials by consolidating powders. The driving force to transform an initial particulate material into a bulk material is the reduction of the interfacial energy of the system. The seminal experimental work of Petzow and Exner showed that particle rearrangement is an important feature of solid-state sintering for crystalline and amorphous powders Petzow et al. 1976. More recent studies using X-ray tomography confirmed that translational, rolling and intrinsic rotation movements of particles play an important role at an early stage of sintering Grupp et al. 2011. These works substantiated the argument that during sintering, powders cannot be considered as a continuum and that the discrete nature of the initial material is maintained. This motivated the use of discrete simulations to take explicitly into account the particulate nature of the processed material. The Discrete Element Method (DEM), is well adapted to this task as it can handle the interactions of a large number of discrete particles. It has been used extensively to model sintering over the past 20 years Henrich et al. 2007; C. L. Martin, Schneider, et al. 2006; Nosewicz, Rojek, Pietrzak, et al. 2013; Wang et al. 2013; S. Martin, Guessasma, et al. 2014; Besler et al. 2015; S. Martin, Navarro, et al. 2016; Rasp et al. 2017; Nosewicz, Rojek, and Chmielewski 2020; Teixeira et al. 2021; Paredes-goyes et al. 2021; Paredes-Goyes et al. 2022 For example, the evolution of important microstructural parameters such as pore networks, coordination number or grain sizes can be investigated at the length scale of particles Wang et al. 2013; S. Martin, Guessasma, et al. 2014; Besler et al. 2015; S. Martin, Navarro, et al. 2016; Paredes-goyes et al. 2021. More recently, the beneficial effect of high heating-rates has been better understood thanks to DEM simulations Teixeira et al. 2021; Paredes-Goyes et al. 2022 For simplicity and to limit CPU time, all these studies represent particles as spheres. The main reason is that contact detection between two spheres is fast and can be fairly easily optimized, even for broad particle size distribution Ogarko et al. 2012.

<footnote n="1">¹ Univ. Grenoble-Alpes, CNRS, Grenoble INP, SIMaP, Grenoble F-38000, France</footnote>

<footnote n="2">² IMT Mines Albi, Université de Toulouse, CNRS, Centre RAPSODEE, Albi F-81013, France. bparedes@mines-albi.fr</footnote>
Identifiers

DOI 10.46298/jtcam.13721

HAL hal-04594920v2

History

Received Jun 4, 2024

Accepted Dec 16, 2024

Published Apr 25, 2025

**Associate Editor**

Phu NGUYEN

Reviewers

Tushar Kanti MANDAL

Anonymous

Open Review
HAL hal-05021530

Supplementary Material

Software

SWH dp3D

Data

DOI 10.5281/zenodo.15168970

Licence

CC BY 4.0

©The Authors

# A Level Set Discrete Element Model for sintering with an optimization-based contact detection

Brayan PAREDES-GOYES¹,², David JAUFFRES¹, and Christophe MARTIN¹

¹ Université Grenoble-Alpes, CNRS, Grenoble INP, SIMaP, Grenoble, France

² IMT Mines Albi, Université de Toulouse, CNRS, Centre RAPSODEE, Albi, France

Sintering is a high temperature process for the consolidation of ceramic, metal and polymer powders. The Discrete Element Method (DEM) has been effectively used to model the sintering process at the particle scale considering spherical particles. However, standard manufacturing processes rarely deal with spherical particles. As sintering is a curvature-controlled process, it is important to take into account the deviation from sphericity. This study presents a DEM sintering model for non-spherical particles. The description and dynamic evolution of arbitrary shape particles is achieved by using the Level Set Discrete Element Method (LS-DEM). The original LS-DEM approach uses boundary nodes on the particles to detect contacts. We employ an optimization-based contact detection approach. This improves the capture of small contacts, which is important for a correct description of sintering evolution with reasonable CPU-time consumption. A Newton-Raphson scheme is employed for the optimization algorithm. The normal force and neck size evolution expressions of spherical particles are adapted for arbitrary shape particles by using the local curvature at the contact. The developed model is validated for elastic contacts on superquadric ellipsoids. It is compared with standard DEM on spheres for sintering. The model is applied to investigate the consolidation kinetics of a packing of ellipsoidal particles. It is shown, that a deviation from sphericity is beneficial for both prolate and oblate ellipsoids. An optimum aspect ratio is evaluated, demonstrating that particles that are too elongated slow down densification kinetics.

Keywords: discrete element method, level set, sintering, non-spherical particles, contact detection

## 1 Introduction

Sintering is a prominent high-temperature process to manufacture ceramic, metallic and polymeric materials by consolidating powders. The driving force to transform an initial particulate material into a bulk material is the reduction of the interfacial energy of the system. The seminal experimental work of Petzow and Exner showed that particle rearrangement is an important feature of solid-state sintering for crystalline and amorphous powders (Petzow and Exner 1976). More recent studies using X-ray tomography confirmed that translational, rolling and intrinsic rotation movements of particles play an important role at an early stage of sintering (Grupp et al. 2011). These works substantiated the argument that during sintering, powders cannot be considered as a continuum and that the discrete nature of the initial material is maintained. This motivated the use of discrete simulations to take explicitly into account the particulate nature of the processed material. The Discrete Element Method (DEM), is well adapted to this task as it can handle the interactions of a large number of discrete particles. It has been used extensively to model sintering over the past 20 years (Henrich et al. 2007; Martin et al. 2006; Martin and Bordia 2009; Nosewicz et al. 2013; Wang and Chen 2013; Martin et al. 2014; Besler et al. 2015; Martin et al. 2016; Rasp et al. 2017; Nosewicz et al. 2020; Teixeira et al. 2021; Paredes-Goyes et al. 2021; Paredes-Goyes et al. 2022) For example, the evolution of important microstructural parameters such as pore networks, coordination number or grain sizes can be investigated at the length scale of particles (Wang and Chen 2013; Martin et al. 2014; Besler et al. 2015; Martin et al. 2016; Paredes-Goyes et al. 2021). More recently, the beneficial effect of high heating-rates has been
Nevertheless, actual particles are not spherical. For example, SEM images of typical alumina powders show that the particle shape is far from being spherical (Fig. 1a). Even particles in advanced alumina manufacturing process, where the size, size distribution and shape are controlled, are not perfect spheres (Fig. 1b). The shape can influence the microstructural evolution as sintering is a curvature-driven process. For instance, kinetic Monte-Carlo simulations have found that packings with higher aspect ratio particles reach higher densities Yan et al. 2017; Bjørk 2022.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** SEM images of powders of a) alumina from general manufacturing process, b) alumina from advanced manufacturing process (hydrolysis of aluminum alkoxide) Venkatesh 2023.</figcaption></figure>

In the DEM framework, the most common representations of non-spherical particles, as summarized by Lu et al. Lu et al. 2015, are multi-spheres, ellipsoids, polyhedral, superquadrics, the combination of geometric elements, and potential particles. Multi-spheres is a widely used method for representing arbitrary shape particles. Its use has been attempted for sintering Hugonnet et al. 2020 but it is not convenient since the local curvature may not be correctly captured. More recently, outside the field of sintering, new developments have been proposed in terms of level-set method Kawamoto et al. 2016, surface meshes Zhan et al. 2020 and Fourier series Lai et al. 2020.

The level-set method was pioneered by Dervieux and Thomasset Dervieux et al. 1980 and by Osher and Sethian Osher et al. 1988. It uses a scalar function to represent a close surface in 3D. The level-set function is zero for any point on the particle surface, negative inside and positive outside. Coupling DEM and level-set method (LS-DEM) Kawamoto et al. 2016 is an interesting approach that captures arbitrary shapes using level-sets, while keeping the discontinuous framework of DEM. This method allows real particle shapes to be obtained directly from 3D tomography images. Although LS-DEM computational cost is substantially higher than DEM, the method is tractable for tens or hundreds of thousands of particles Kawamoto et al. 2016. A detailed comparison between LS-DEM and DEM on accuracy and computational cost has been recently proposed Duriez et al. 2021. LS-DEM has already been applied for triaxial compaction tests Kawamoto et al. 2016, breakage mechanics Harmon, Arthur, et al. 2020; Pazmiño et al. 2022, electrostatic cohesion Bustamante et al. 2020, prediction of shear banding Kawamoto et al. 2018, the investigation of incremental behavior of granular materials Karapiperis et al. 2020 and particle bonding Harmon, Karapiperis, et al. 2021. The above studies have in common that they focus on geomaterials, which feature elastic interactions between particles. For engineering materials, specific contact laws adapted for DEM have been proposed for sintering spherical particles C. L. Martin, Schneider, et al. 2006; Henrich et al. 2007; Paredes-goyes et al. 2021. They introduce material parameters such as diffusion coefficients and surface energies. So far, and to the best of our knowledge, LS-DEM has not been applied to the sintering process. The aim of this work is to present a LS-DEM framework that is compatible both with elasto-plastic interactions and sintering interactions.

Contact detection is the most challenging stage of simulating non-spherical particles. LS-DEM, as proposed originally by Kawamoto et al. Kawamoto et al. 2016 performs this task by creating surface nodes on the particle surface and evaluating if they are inside another particle. We
better understood thanks to DEM simulations (Teixeira et al. 2021; Paredes-Goyes et al. 2022) For simplicity and to limit CPU time, all these studies represent particles as spheres. The main reason is that contact detection between two spheres is fast and can be fairly easily optimized, even for broad particle size distribution (Ogarko and Luding 2012).

Nevertheless, actual particles are not spherical. For example, SEM images of typical alumina powders show that the particle shape is far from being spherical, see Figure 1(a). Even particles in advanced alumina manufacturing process, where the size, size distribution and shape are controlled, are not perfect spheres, see Figure 1(b). The shape can influence the microstructural evolution as sintering is a curvature-driven process. For instance, kinetic Monte-Carlo simulations have found that packings with higher aspect ratio particles reach higher densities (Yan et al. 2017; Bjørk 2022).

<figure><img src="image_2.png" /><figcaption>**Figure 1** SEM images of powders of (a) alumina from general manufacturing process, (b) alumina from advanced manufacturing process (hydrolysis of aluminum alkoxide) (Venkatesh 2023).</figcaption></figure>

In the DEM framework, the most common representations of non-spherical particles, as summarized by Lu et al. (2015), are multi-spheres, ellipsoids, polyhedral, superquadrics, the combination of geometric elements, and potential particles. Multi-spheres is a widely used method for representing arbitrary shape particles. Its use has been attempted for sintering (Hugonnet et al. 2020) but it is not convenient since the local curvature may not be correctly captured. More recently, outside the field of sintering, new developments have been proposed in terms of level-set method (Kawamoto et al. 2016), surface meshes (Zhan et al. 2020) and Fourier series (Lai et al. 2020).

The level-set method was pioneered by Dervieux and Thomasset (1980) and Osher and Sethian (1988). It uses a scalar function to represent a close surface in 3D. The level-set function is zero for any point on the particle surface, negative inside and positive outside. Coupling DEM and level-set method (LS-DEM) (Kawamoto et al. 2016) is an interesting approach that captures arbitrary shapes using level-sets, while keeping the discontinuous framework of DEM. This method allows real particle shapes to be obtained directly from 3D tomography images. Although LS-DEM computational cost is substantially higher than DEM, the method is tractable for tens or hundreds of thousands of particles (Kawamoto et al. 2016). A detailed comparison between LS-DEM and DEM on accuracy and computational cost has been recently proposed (Duriez and Bonelli 2021). LS-DEM has already been applied for triaxial compaction tests (Kawamoto et al. 2016), breakage mechanics (Harmon et al. 2020; Pazmiño et al. 2022), electrostatic cohesion (Bustamante et al. 2020), prediction of shear banding (Kawamoto et al. 2018), the investigation of incremental behavior of granular materials (Karapiperis et al. 2020) and particle bonding (Harmon et al. 2021). The above studies have in common that they focus on geomaterials, which feature elastic interactions between particles. For engineering materials, specific contact laws adapted for DEM have been proposed for sintering spherical particles (Martin et al. 2006; Henrich et al. 2007; Paredes-Goyes et al. 2021). They introduce material parameters such as diffusion coefficients and surface energies. So far, and to the best of our knowledge, LS-DEM has not been applied to the sintering process. The aim of this work is to present a LS-DEM framework that is compatible both with elasto-plastic interactions and sintering interactions.

Contact detection is the most challenging stage of simulating non-spherical particles. LS-DEM, as proposed originally by Kawamoto et al. (2016) performs this task by creating surface nodes on the particle surface and evaluating if they are inside another particle. We have tested this
have tested this technique for a packing of relative density 0.59 with slightly non-uniform size
particles, Fig. 2 indicates that even 40 000 nodes only detect around 70% of the existing contacts.
For denser packings (0.64 and 0.69 relative density) that lead to larger indentations, we observed
that 4 000 nodes lead to 99% of the contacts detected. This indicates that the contact detection
algorithm mainly misses small contacts. These results are coherent with those of Duriez et
al. Duriez et al. 2021 who have shown that 10 000 nodes lead to a 15% underestimation of the
macroscopic pressure (the number of missed contacts was not indicated). In a packing where
interactions are elastic, missing small contacts (30% for 40 000 nodes) has only a limited impact
on the macroscopic pressure since the pressure is linearly related to contact forces (Love equation
Christoffersen et al. 1981). However, for sintering, small contacts are associated with relatively
large tensile forces driven by surface energy minimization (as detailed in section 2.4). Thus, small
contact detection is critical for a correct description of the packing macroscopic behavior, in
particular, at the initial stage of sintering. At this stage, where defects may arise (in the form of
non-homogeneous sintering), it is important to correctly describe the densification rate and the
rearrangement of particles that are mainly controlled by the tensile forces between particles.

Increasing the number of surface nodes is not a viable option as it becomes CPU prohibitive
above 10 000 nodes. A recent LS-DEM study Davis et al. 2021 applied an optimization-based
contact detection in order to decrease the computational cost and eliminate the dependence
of the force on the number of surface nodes. While they used a derivative-free optimization
algorithm, we opt for a Newton-Raphson method because of its fast local convergence and its
ease of calculating derivatives from the discrete Level Set function, which will also be used for
curvature calculation. This approach is based on the two-contact points search proposed by
Houlsby Houlsby 2009 in the context of what he is defining as 2D potential particles (and what
we define as shape function). The author developed it for convex particles, however, we add here
a multi-start strategy for contact detection among convex particles, that could also be used to find
multiple contact points in non-convex particles.

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Contact detection accuracy as a function of the number of surface nodes using the original LS-DEM method for a packing of spherical particles with a density of 0.59. The accuracy is checked against a standard DEM scheme with double precision.</figcaption></figure>

Our work presents a LS-DEM model of sintering with an optimization-based contact detection.
Section 2 describes the model, detailing the contact detection scheme. Section 3 validates the
developed model for elastic and sintering interactions. Section 4 models the sintering of a
two-particle system and of a packing of particles. As the aim is to provide a proof of concept of
our scheme, we have limited applications to ellipsoid shapes.
technique for a packing of relative density 0.59 with slightly non-uniform size particles, Figure 2 indicates that even 40 000 nodes only detect around 70% of the existing contacts. For denser packings (0.64 and 0.69 relative density) that lead to larger indentations, we observed that 4 000 nodes lead to 99% of the contacts detected. This indicates that the contact detection algorithm mainly misses small contacts, i.e., contacts that present a small inter-particle indentation. These results are coherent with those of Duriez and Bonelli (2021) who have shown that 10 000 nodes lead to a 15% underestimation of the macroscopic pressure (the number of missed contacts was not indicated). In a packing where interactions are elastic, missing small contacts (30% for 40 000 nodes) has only a limited impact on the macroscopic pressure since the pressure is linearly related to contact forces (Love equation (Christoffersen et al. 1981)). However, for sintering, small contacts are associated with relatively large tensile forces driven by surface energy minimization (as detailed in Section 2.4). Thus, small contact detection is critical for a correct description of the packing macroscopic behavior, in particular, at the initial stage of sintering. At this stage, where defects may arise (in the form of non-homogeneous sintering), it is important to correctly describe the densification rate and the rearrangement of particles that are mainly controlled by the tensile forces between particles.

Increasing the number of surface nodes is not a viable option as it becomes CPU prohibitive above 10 000 nodes. A recent LS-DEM study (Davis et al. 2021) applied an optimization-based contact detection in order to decrease the computational cost and eliminate the dependence of the force on the number of surface nodes. While they used a derivative-free optimization algorithm, we opt for a Newton-Raphson method because of its fast local convergence and its ease of calculating derivatives from the discrete Level Set function, which will also be used for curvature calculation. This approach is based on the two-contact points search proposed by Houlsby (2009) in the context of what he is defining as 2D potential particles (and what we define as shape function). The author developed it for convex particles, however, we add here a multi-start strategy for contact detection among convex particles, that could also be used to find multiple contact points in non-convex particles.

<figure><img src="image_3.png" /><figcaption>**Figure 2** Contact detection accuracy as a function of the number of surface nodes using the original LS-DEM method for a packing of spherical particles with a density of 0.59. The accuracy is checked against a standard DEM scheme with double precision.</figcaption></figure>

Our work presents a LS-DEM model of sintering with an optimization-based contact detection. Section 2 describes the model, detailing the contact detection scheme. Section 3 validates the developed model for elastic and sintering interactions. Section 4 models the sintering of a two-particle system and of a packing of particles. As the aim is to provide a proof of concept of our scheme, we have limited applications to ellipsoid shapes.

# 2 Model description

## 2.1 Equations of motion

A dynamic DEM approach is adopted. The evolution of the position of the center of mass $\vec{x}_i$, for each particle $i$ follows Newton's second law

$$m_i \ddot{\vec{x}}_i = \sum \vec{F}_{ij} \qquad (1)$$
# 2 Model description

## 2.1 Equations of motion

A dynamic DEM approach is adopted. The evolution of the position of the center of mass $\vec{x}_i$, for each particle $i$ follows Newton's second law:

$$m_i \ddot{\vec{x}}_i = \sum \vec{F}_{ij} \qquad (1)$$

where $m_i$ is the particle mass and $\vec{F}_{ij}$ the total force (with normal and tangential components in the local contact framework) exerted by a particle $j$ in contact. The particle position is updated by imposing in the first half time step an affine displacement which follows the macroscopic imposed strain-rate Thornton et al. 1998 and then integrating Eq. (1) in the second half time step using a Velocity Verlet algorithm as detailed in C. L. Martin and Bordia 2009.

For non-spherical particles, even in the absence of frictional forces, rotations need to be considered. Although we treat the motion of particles dynamically using Eq. (1), we approximate sintering as a quasi-static process, in this case the nonlinear term of the general equation of rotational motion can be neglected obtaining the following expression for the angular position $\vec{\theta}_i$ Hart et al. 1988:

$$I_i \ddot{\vec{\theta}}_i = \sum \vec{M}_{ij} \qquad (2)$$

with the particle moment inertia $I_i$ and the torque $\vec{M}_{ij}$. In addition, the quasi-static assumption allows the non-diagonal terms of the moment of inertia tensor to be neglected Hart et al. 1988; Salque 2017. The particle mass and moment of inertia are estimated as proposed in the original LS-DEM methodology Kawamoto et al. 2016.

Quaternions are used to track rotations of particles. The advantage (as compared to Euler angles) is to avoid singularities at small angles Džiugys et al. 2001. A quaternion $\mathbf{q}$, defined by an angle $\theta$ and an orientation vector $\vec{u}$, is attached to each particle:

$$\mathbf{q} = (q_0, q_1, q_2, q_3)^T = \cos(\theta(t)/2) + \sin(\theta(t)/2) [u_x \vec{x} + u_y \vec{y} + u_z \vec{z}] \qquad (3)$$

At each time step $\Delta t$, its variation is calculated depending on the rotational velocity at full time step $\vec{\theta}_{t+\Delta t}$:

$$\mathbf{q}_{\Delta t} = \cos\left(\frac{|\vec{\theta}_{t+\Delta t}|\Delta t}{2}\right) + \sin\left(\frac{|\vec{\theta}_{t+\Delta t}|\Delta t}{2}\right) \frac{\vec{\theta}_{t+\Delta t}}{|\vec{\theta}_{t+\Delta t}|} \qquad (4)$$

The updated orientation is computed from the quaternion product: $\mathbf{q}_{t+\Delta t} = \mathbf{q}_t \mathbf{q}_{\Delta t}$ Wachs et al. 2012.

## 2.2 Discrete level-set function and derivatives

The particle shape is described by the signed distance to the particle surface (level-set function) stored at grid points Kawamoto et al. 2016. The level-set function can be computed analytically from the distance equation to a 3D shape or obtained from tomography images of the packing of particles Vlahinić et al. 2014. The bounding sphere of the particle is used to build a local uniform Cartesian grid with $ncell_{int}$ number of cells in each direction (Fig. 3). $ncell_{ext}$ extra cells of the same size are added in all directions for contact detection purposes (see section 2.3). The level-set values are stored at all these grid points (green dots in Fig. 3).

During the contact detection and calculation stages, the level-set values and their derivatives are required between the grid points. Linear interpolation is used to calculate the level-set value $\phi$ at any point $p$ from the surrounding $abc$ grid points with level-set values $\phi_{abc}$ Kawamoto et al. 2016:

$$\phi(p) = \sum_{a=0}^{1} \sum_{b=0}^{1} \sum_{c=0}^{1} \phi_{abc} [(1-a)(1-x) + ax] [(1-b)(1-y) + by] [(1-c)(1-z) + cz] \qquad (5)$$
where $m_i$ is the particle mass and $\vec{F}_{ij}$ the total force (with normal and tangential components in the local contact framework) exerted by a particle j in contact. The particle position is updated by imposing in the first half time step an affine displacement which follows the macroscopic imposed strain-rate (Thornton and Ning 1998) and then integrating Equation (1) in the second half time step using a Velocity Verlet algorithm as detailed in (Martin and Bordia 2009).

For non-spherical particles, even in the absence of frictional forces, rotations need to be considered. Although we treat the motion of particles dynamically using Equation (1), we approximate sintering as a quasi-static process, in this case the nonlinear term of the general equation of rotational motion can be neglected thus obtaining the expression

$$I_i \ddot{\theta}_i = \sum \vec{M}_{ij} \qquad (2)$$

for the angular position $\vec{\theta}_i$ (Hart et al. 1988), with the particle moment inertia $I_i$ and the torque $\vec{M}_{ij}$. In addition, the quasi-static assumption allows the non-diagonal terms of the moment of inertia tensor to be neglected (Hart et al. 1988; Salque 2017). The particle mass and moment of inertia are estimated as proposed in the original LS-DEM methodology (Kawamoto et al. 2016).

Quaternions are used to track rotations of particles. The advantage (as compared to Euler angles) is to avoid singularities at small angles (Džiugys and Peters 2001). A quaternion $q$, defined by an angle $\theta$ and an orientation vector $\vec{u}$, is attached to each particle:

$$q = (q_0, q_1, q_2, q_3)^T = \cos(\theta(t)/2) + \sin(\theta(t)/2) [u_x \vec{x} + u_y \vec{y} + u_z \vec{z}]. \qquad (3)$$

At each time step $\Delta t$, its variation is calculated depending on the rotational velocity at full time step $\dot{\vec{\theta}}_{t+\Delta t}$:

$$q_{\Delta t} = \cos\left(\frac{|\dot{\vec{\theta}}_{t+\Delta t} | \Delta t}{2}\right) + \sin\left(\frac{|\dot{\vec{\theta}}_{t+\Delta t} | \Delta t}{2}\right) \frac{\dot{\vec{\theta}}_{t+\Delta t}}{|\dot{\vec{\theta}}_{t+\Delta t}|}. \qquad (4)$$

The updated orientation is computed from the quaternion product $q_{t+\Delta t} = q_t q_{\Delta t}$ (Wachs et al. 2012).

## 2.2 Discrete level-set function and derivatives

The particle shape is described by the signed distance to the particle surface (level-set function) stored at grid points (Kawamoto et al. 2016). The level-set function can be computed analytically from the distance equation to a 3D shape or obtained from tomography images of the packing of particles (Vlahinić et al. 2014). The bounding sphere of the particle is used to build a local uniform Cartesian grid with $ncell_{int}$ number of cells in each direction (Figure 3). $ncell_{ext}$ extra cells of the same size are added in all directions for contact detection purposes, see Section 2.3. The level-set values are stored at all these grid points indicated as green dots in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3 2D schematic of a local grid of a particle with grid points in green. The grid is composed by $ncell_{int}$ + $ncell_{ext}$ cells in each direction. The bounding sphere of the particle is represented by the dashed circle.</figcaption></figure>

During the contact detection and calculation stages, the level-set values and their derivatives are required between the grid points. Linear interpolation is used to calculate the level-set value
<figure><img src="image_5.png" /><figcaption>Figure 3: 2D schematic of a local grid of a particle with grid points in green. The grid is composed by ncell$_{int}$+ncell$_{ext}$ cells in each direction. The bounding sphere of the particle is represented by the dashed circle.</figcaption></figure>

Here for simplicity, the first and second derivatives of the level-set functions are computed by central finite differences (FD). To prevent the second derivative from vanishing due to the first-degree polynomial $\phi(p)$, a FD step size greater than the grid cell size is adopted. We verified that these approximations do not have any noticeable impact on the final results.

In contrast with other LS-DEM works, our contact law needs the input of the local curvature (see section 2.4). We use the mean curvature instead of the Gaussian curvature to avoid singularities of the curvature radius as pointed out by Podlozhnyuk et al. Podlozhnyuk et al. 2017 for superquadric particles. The mean curvature $\kappa$ of level-set functions is used by adopting the formula proposed in Goldman 2005 for implicit surfaces. $\kappa$ is written in terms of first and second derivatives of a particle i:

$$ \kappa = \frac{\phi_{ix}^2 \phi_{iyy} - 2\phi_{ix}\phi_{iyy}\phi_{ixy} + \phi_{iy}^2 \phi_{ixx} + \phi_{ix}^2 \phi_{izz} - 2\phi_{ix}\phi_{iz}\phi_{ixz} + \phi_{iz}^2 \phi_{ixx} + \phi_{iy}^2 \phi_{izz} - 2\phi_{iy}\phi_{iz}\phi_{iyz} + \phi_{iz}^2 \phi_{iyy}}{2|\nabla\phi_i|^3} \quad (6) $$

where $\phi_{ik} = \frac{\partial\phi_i}{\partial k}$ and $\phi_{ikl} = \frac{\partial^2\phi_i}{\partial k \partial l}$.

## 2.3 Contact detection

The search algorithm for contacts is based on an optimization approach. We extend to 3D the method developed by Houlsby Houlsby 2009 for 2D particles. In his case, the shape of particles was described by an analytical expression. In our case, we adapt the method to the level-set function discretized on the grid.

For a given particle, its level-set function is negative inside the particle, zero on the surface, and positive outside as shown in Fig. 4. Let i and j be two particles in potential contact. The potential contact list is given by the list of contacts between bounding spheres. The first step is to find (if it exists) the innermost point in j while imposing that it belongs to the surface of i ($\vec{x}_{op,\phi i}$
$\phi$ at any point $p$ from the surrounding $abc$ grid points with level-set values $\phi_{abc}$ (Kawamoto et al. 2016):

$$ \phi(p) = \sum_{a=0}^{1} \sum_{b=0}^{1} \sum_{c=0}^{1} \phi_{abc} [(1-a)(1-x)+ax] [(1-b)(1-y)+by] [(1-c)(1-z)+cz]. \quad (5) $$

Here for simplicity, the first and second derivatives of the level-set functions are computed by central finite differences (FD). To prevent the second derivative from vanishing due to the first-degree polynomial $\phi(p)$, a FD step size greater than the grid cell size is adopted. We verified that these approximations do not have any noticeable impact on the final results.

In contrast with other LS-DEM works, our contact law needs the input of the local curvature (see Section 2.4). We use the mean curvature instead of the Gaussian curvature to avoid singularities of the curvature radius as pointed out by Podlozhnyuk et al. (2017) for superquadric particles. The mean curvature $\kappa$ of level-set functions is used by adopting the formula proposed for implicit surfaces in (Goldman 2005) and employing the sign convention of Osher et al. (2004) in the context of level-sets. $\kappa$ is written in terms of first and second derivatives of a particle $i$:

$$ \kappa = \frac{\phi_{ix}^2 \phi_{iyy} - 2\phi_{ix}\phi_{iyy}\phi_{iyy} + \phi_{iyy}^2 \phi_{ixx} + \phi_{ixx}^2 \phi_{iyy} - 2\phi_{ix}\phi_{iyy}\phi_{iyy} + \phi_{iyy}^2 \phi_{ixx} + \phi_{ixx}^2 \phi_{iyy} - 2\phi_{iyy}\phi_{iyy}\phi_{iyy} + \phi_{iyy}^2 \phi_{ixx} + \phi_{ixx}^2 \phi_{iyy} - 2\phi_{iyy}\phi_{iyy}\phi_{iyy} + \phi_{iyy}^2 \phi_{ixx} + \phi_{ixx}^2 \phi_{iyy}}{2|\nabla \phi_i|^3} \quad (6) $$

where $\phi_{ik} = \frac{\partial \phi_i}{\partial k}$ and $\phi_{ikl} = \frac{\partial^2 \phi_i}{\partial k \partial l}$.

## 2.3 Contact detection

The search algorithm for contacts is based on an optimization approach. We extend to 3D the method developed by Houlsby (2009) for 2D particles. In his case, the shape of particles was described by an analytical expression. In our case, we adapt the method to the level-set function discretized on the grid.

For a given particle, its level-set function is negative inside the particle, zero on the surface, and positive outside as shown in Figure 4. Let $i$ and $j$ be two particles in potential contact. The

<figure><img src="image_5.png" /><figcaption>**Figure 4** Optimization-based contact detection for two particles $i$ and $j$ described by level-set function $\phi$. $\vec{x}_{op,\phi_i}$ is the innermost point in $j$ on the surface of $i$ (first step). $\vec{x}_{op,\phi_j}$ is the innermost point in $i$ on the surface of $j$ (second step).</figcaption></figure>

potential contact list is given by the list of contacts between bounding spheres. The first step is to find (if it exists) the innermost point in $j$ while imposing that it belongs to the surface of $i$ ($\vec{x}_{op,\phi_i}$ in Figure 4). This point can be found by solving the constrained optimization problem

$$ \begin{array}{ll} \min_{x,y,z} & \phi_j(x, y, z) \\ \text{s.t.} & \phi_i(x, y, z) = 0 \end{array} \quad (7) $$

where s.t. stands for subject to. If the point found is inside particle $j$, i.e., $\phi_j(\vec{x}_{op,\phi_i}) < 0$, the two particles are in contact. In that case, a similar optimization problem is solved to find the point $\vec{x}_{op,\phi_j}$ (Figure 4):

$$ \begin{array}{ll} \min_{x,y,z} & \phi_i(x, y, z) \\ \text{s.t.} & \phi_j(x, y, z) = 0. \end{array} \quad (8) $$
<figure><img src="image_6.png" /><figcaption>Figure 4: Optimization-based contact detection for two particles i and j described by level-set function $\phi$. $\vec{x}_{op,\phi i}$ is the innermost point in j on the surface of i (first step). $\vec{x}_{op,\phi j}$ is the innermost point in i on the surface of j (second step).</figcaption></figure>

in Fig. 4). This point can be found by solving the following constrained optimization problem:

$$ 
\begin{array}{ll}
\min_{x,y,z} & \phi_j(x, y, z) \\
\text{s.t.} & \phi_i(x, y, z) = 0
\end{array}
\qquad (7)
$$

where s.t. stands for subject to. If the point found is inside particle j, i.e., $\phi_j(\vec{x}_{op,\phi i}) < 0$, the two particles are in contact. In that case, a similar optimization problem is solved to find the point $\vec{x}_{op,\phi j}$ (Fig. 4):

$$ 
\begin{array}{ll}
\min_{x,y,z} & \phi_i(x, y, z) \\
\text{s.t.} & \phi_j(x, y, z) = 0
\end{array}
\qquad (8)
$$

This second optimization step is needed to compute the required contact parameters Houlsby 2009 (section 2.4). It is also used as a back-up to ensure that the contact actually exists. This is because in a few cases, inaccuracies can arise due to the discretized character of the level-set function on the grid.

Solving the optimization problems in Eqs. (7) and (8) separately, Houlsby 2009 applied in 2D the method of Lagrangian multipliers to obtain a set of nonlinear equations, which he solved with the Newton-Raphson (N-R) method. Here, the same approach is adopted. The linear system of one N-R iteration that we have obtained for the optimization problem (7) is in 3D:

$$ 
\left( 
\begin{array}{c}
\phi_{ix} \\
\phi_{ixx}\phi_{jy} + \phi_{ix}\phi_{jxy} - \\
\phi_{jxx}\phi_{iy} - \phi_{jx}\phi_{ixy} \\
\phi_{ixx}\phi_{jz} + \phi_{ix}\phi_{jxz} - \\
\phi_{jxx}\phi_{iz} - \phi_{jx}\phi_{ixz}
\end{array}
\right)
=
\left( 
\begin{array}{c}
\phi_{iy} \\
\phi_{ixy}\phi_{jy} + \phi_{ix}\phi_{jyy} - \\
\phi_{jxy}\phi_{iy} - \phi_{jx}\phi_{iyy} \\
\phi_{ixy}\phi_{jz} + \phi_{ix}\phi_{jyz} - \\
\phi_{jxy}\phi_{iz} - \phi_{jx}\phi_{iyz}
\end{array}
\right)
-
\left( 
\begin{array}{c}
\phi_{iz} \\
\phi_{ixz}\phi_{jy} + \phi_{ix}\phi_{jyz} - \\
\phi_{jxz}\phi_{iy} - \phi_{jx}\phi_{iyz} \\
\phi_{ixz}\phi_{jz} + \phi_{ix}\phi_{jzz} - \\
\phi_{jxz}\phi_{iz} - \phi_{jx}\phi_{izz}
\end{array}
\right)
\qquad (9)
$$
This second optimization step is needed to compute the required contact parameters (Houlsby 2009), Section 2.4. It is also used as a back-up to ensure that the contact actually exists. This is because in a few cases, inaccuracies can arise due to the discretized character of the level-set function on the grid.

Solving the optimization problems in Equations (7) and (8) separately, Houlsby (2009) applied in 2D the method of Lagrangian multipliers to obtain a set of nonlinear equations, which he solved with the Newton-Raphson (N-R) method. Here, the same approach is adopted. In 3D, The linear system of one N-R iteration obtained for the optimization problem (7) is

$$ \begin{pmatrix} \phi_{ix} & \phi_{iy} & \phi_{iz} \\ [\phi_{ixx}\phi_{jy} + \phi_{ix}\phi_{jxy} & [\phi_{ixy}\phi_{jy} + \phi_{ix}\phi_{jyy} & [\phi_{ixz}\phi_{jy} + \phi_{ix}\phi_{jyz} \\ -\phi_{jxx}\phi_{iy} - \phi_{jx}\phi_{ixy}] & -\phi_{jxy}\phi_{iy} - \phi_{jx}\phi_{iyy}] & -\phi_{jxz}\phi_{iy} - \phi_{jx}\phi_{iyz}] \\ [\phi_{ixx}\phi_{jz} + \phi_{ix}\phi_{jxz} & [\phi_{ixy}\phi_{jz} + \phi_{ix}\phi_{jyz} & [\phi_{ixz}\phi_{jz} + \phi_{ix}\phi_{jzz} \\ -\phi_{jxx}\phi_{iz} - \phi_{jx}\phi_{ixz}] & -\phi_{jxy}\phi_{iz} - \phi_{jx}\phi_{iyz}] & -\phi_{jxz}\phi_{iz} - \phi_{jx}\phi_{izz}] \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta y \\ \Delta z \end{pmatrix} = \begin{pmatrix} -\phi_i \\ -\phi_{ix}\phi_{jy} + \phi_{jx}\phi_{iy} \\ -\phi_{ix}\phi_{jz} + \phi_{jx}\phi_{iz} \end{pmatrix} \quad (9) $$

or, in compact form, $\vec{J}\vec{\Delta}_x = -\Phi$ where $\vec{J}$ is the Jacobian matrix. At iterations $nr$, the new point is $\vec{x}_{nr+1} = \vec{x}_{nr} + \alpha_{nr}\vec{\Delta}_x$ with $\alpha_{nr}$ a scalar ($\alpha_{nr} \le 1$). The optimal point $\vec{x}_{op,\phi_i}$ is obtained after a few N-R iterations. The scalar parameter $\alpha_{nr}$ ensures stability and is calculated as in (Podlozhnyuk et al. 2017). The initial point $\vec{x}_{0sph}$ is chosen as the barycenter of the two spheres that circumscribe the particles. The barycenter initial point for N-R is kept as long as no contact is found. Once a contact is detected, the last value of $\vec{x}_{op,\phi_i}$ is kept for the next time step as the initial point for N-R.

In very few instances, we have observed that the above procedure does not detect existing contacts. This problem is due to the possibility of local Newton-Raphson convergence. To make contact detection more robust, we implement a simple multi-start strategy, see Figure 5. If the $\vec{x}_{0sph}$ initial point does not lead to a contact, $N_{ms}$ new initial points are generated and tested for the N-R search until a contact is detected. If no contact is found after $N_{ms}$ attempts, we consider that the contact does not exist. The $N_{ms}$ new initial points are generated with the constraints that they are inside a sphere with origin in $\vec{x}_{0sph}$ and radius $r_{ms} = 0.5 \min(r_i, r_j)$ and inside particle $i$ (and similarly inside particle $j$ to find the $\vec{x}_{op,\phi_j}$ point). There is a compromise between the effectiveness of the method and the computational cost. Our simulations have shown that $N_{ms} = 10$ gives a good balance for ellipsoids. The total number of grid cells, $ncell_{int} + ncell_{ext}$ in Figure 3, also affects this choice since the grid size impacts the precision of the N-R search.

The choice of $N_{ms}$ is more challenging for particles beyond ellipsoids. One possible approach is to carry out tests for representative pairs of particles to assess the ability to detect a contact and the computation time for different predefined $N_{ms}$ values. The use of circumscribed boxes instead of spheres could also reduce the number of $N_{ms}$ required by using the additional information provided by the axes of the boxes without reducing the robustness of contact detection.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Multi-start strategy for ensuring contact detection. If the initial point $\vec{x}_{0sph}$ (barycenter of the two encompassing spheres) does not lead to a contact, 10 new initial points are generated inside the red sphere with the constraint that they are located inside particle $i$.</figcaption></figure>

Similarly, the contact detection of a particle-plane pair is based on optimization. From the general equation of the plane, a level-set function can be deduced. This allows to formulate the optimization problem in the form of Equation (7), where $i$ is the particle and $j$ the plane.
or in compact form:

$$J\vec{\Delta}_x = -\Phi \qquad (10)$$

where J is the Jacobian matrix. At iterations nr, the new point is $\vec{x}_{nr+1} = \vec{x}_{nr} + \alpha_{nr}\vec{\Delta}_x$ with $\alpha_{nr}$ a scalar ($\alpha_{nr} \le 1$). The optimal point $\vec{x}_{op,\phi i}$ is obtained after a few N-R iterations. The scalar parameter $\alpha_{nr}$ ensures stability and is calculated as in Podlozhnyuk et al. 2017. The initial point $\vec{x}_{0sph}$ is chosen as the barycenter of the two spheres that circumscribe the particles. The barycenter initial point for N-R is kept as long as no contact is found. Once a contact is detected, the last value of $\vec{x}_{op,\phi i}$ is kept for the next time step as the initial point for N-R.

In very few instances, we have observed that the above procedure does not detect existing contacts. This problem is due to the possibility of local Newton-Raphson convergence. To make contact detection more robust, we implement a simple multi-start strategy (see Fig. 5). If the $\vec{x}_{0sph}$ initial point does not lead to a contact, $N_{ms}$ new initial points are generated and tested for the N-R search until a contact is detected. If no contact is found after $N_{ms}$ attempts, we consider that the contact does not exist. The $N_{ms}$ new initial points are generated with the constraints that they are inside a sphere with origin in $\vec{x}_{0sph}$ and radius $r_{ms} = 0.5 \min(r_i, r_j)$ and inside particle i (and similarly inside particle j to find the $\vec{x}_{op,\phi j}$ point). There is a compromise between the effectiveness of the method and the computational cost. Our simulations have shown that $N_{ms} = 10$ gives a good balance for ellipsoids. The total number of grid cells, $ncell_{int} + ncell_{ext}$ in Fig. 3, also affects this choice since the grid size impacts the precision of the N-R search.

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Multi-start strategy for ensuring contact detection. If the initial point $\vec{x}_{0sph}$ (barycenter of the two encompassing spheres) does not lead to a contact, 10 new initial points are generated inside the red sphere with the constraint that they are located inside particle i.</figcaption></figure>

Similarly, the contact detection of a particle-plane pair is based on optimization. From the general equation of the plane, a level-set function can be deduced. This allows to formulate the optimization problem in the form of Eq. (7), where i is the particle and j the plane. We consider only planes perpendicular to the main axes. For illustration, considering a plane perpendicular to the z axis, the following system of equations is obtained:

$$\begin{pmatrix} \phi_{ix} & \phi_{iy} & \phi_{iz} \\ \phi_{ixx} & \phi_{ixy} & \phi_{ixz} \\ \phi_{ixy} & \phi_{iyy} & \phi_{iyz} \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta y \\ \Delta z \end{pmatrix} = \begin{pmatrix} -\phi_i \\ -\phi_{ix} \\ -\phi_{iy} \end{pmatrix} \qquad (11)$$

The simplicity of the level-set function of a plane makes Eq. (11) simpler than particle-particle Eq. (9). Eq. (11) is solved in the same way with a multi-start procedure. The only difference is that for particle-plane it is not necessary to solve the second optimization problem (Eq. (8)).

## 2.4 Sintering contact law

Once a contact between two particles i and j is detected, the contact force is calculated in the next stage. Like other forces, the sintering force needs the values of the overlap or indentation $h_{ij}$
<figure><img src="image_7.png" /><figcaption>**Figure 6** Algorithm flowchart for detecting contacts. $\vec{x}_{0sph}$ is the initial point for the constrained optimization problem. It is the barycenter of the two encompassing spheres for $N_{ms} = 1$, see Figures 4 and 5.</figcaption></figure>

We consider only planes perpendicular to the main axes. For illustration, considering a plane perpendicular to the z axis, the following system of equations is obtained:

$$ \begin{pmatrix} \phi_{ix} & \phi_{iy} & \phi_{iz} \\ \phi_{ixx} & \phi_{ixy} & \phi_{ixz} \\ \phi_{ixy} & \phi_{iyy} & \phi_{iyz} \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta y \\ \Delta z \end{pmatrix} = \begin{pmatrix} -\phi_i \\ -\phi_{ix} \\ -\phi_{iy} \end{pmatrix}. \quad (10) $$

The simplicity of the level-set function of a plane makes Equation (10) simpler than particle-particle Equation (9). Equation (10) is solved in the same way with a multi-start procedure. The only difference is that for particle-plane it is not necessary to solve the second optimization problem, see Equation (8).

## 2.4 Sintering contact law

Once a contact between two particles i and j is detected, the contact force is calculated in the next stage. Like other forces, the sintering force needs the values of the overlap or indentation $h_{ij}$ and of the unit normal vector at the contact. Unlike the original LS-DEM method (Kawamoto et al. 2016), we have two points at each contact ($\vec{x}_{op,\phi_i}$ and $\vec{x}_{op,\phi_j}$, Figure 4) instead of one. Therefore, the contact variables are calculated based on these two points. We follow the procedure of potential particles proposed in (Houlsby 2009) that suggests to average the values. The indentation between two particles is given by

$$ h_{ij} = - \frac{\phi_i(\vec{x}_{op,\phi_j}) + \phi_j(\vec{x}_{op,\phi_i})}{2} \quad (11) $$
and of the unit normal vector at the contact. Unlike the original LS-DEM method Kawamoto et al. 2016, we have two points at each contact ($\vec{x}_{op,\phi i}$ and $\vec{x}_{op,\phi j}$, Fig. 4) instead of one. Therefore, the contact variables are calculated based on these two points. We follow the procedure of potential particles proposed in Houlsby 2009 that suggests to average the values. The indentation between two particles is given by:

$$h_{ij} = -\frac{\phi_i(\vec{x}_{op,\phi j}) + \phi_j(\vec{x}_{op,\phi i})}{2} \qquad (12)$$

while the unit normal vector is calculated from:

$$\vec{n}_{ij} = \frac{\nabla\phi_i(\vec{x}_{op,\phi i}) - \nabla\phi_j(\vec{x}_{op,\phi j})}{|\nabla\phi_i(\vec{x}_{op,\phi i}) - \nabla\phi_j(\vec{x}_{op,\phi j})|} \qquad (13)$$

Following the suggestion of Li et al. Li et al. 2019, our force calculation depends on the maximum indentation, instead of the average indentation, which would make it more mesh dependent.

For two spherical particles, the normal force $N_{ij}$ Bouvard et al. 1996; C. L. Martin and Bordia 2009 is:

$$N_{ij} = \frac{\pi a_{ij}^4}{8\Delta_{GB}} \frac{dh_{ij}}{dt} - \frac{\alpha}{2} \pi R_{eq,ij} \gamma_S \qquad (14)$$

where $\alpha$ is a parameter that depends on the ratio between surface and grain boundary diffusion coefficients. It is set to 2.46 for the case of a typical ceramic material (alumina) Bouvard et al. 1996; Paredes-goyes et al. 2021. $\gamma_S$ is the surface energy and $\Delta_{GB} = \frac{\Omega}{k_b T} D_{GB} \delta_{GB}$, with the atomic volume $\Omega$, the Boltzmann constant $k_b$, the temperature $T$, the grain boundary diffusivity $D_{GB}$, the grain boundary thickness $\delta_{GB}$. $R_{eq,ij} = r_i r_j / (r_i + r_j)$ is the equivalent radius of particles i and j with radii $r_i$ and $r_j$. The first term is a viscous component that counteracts the relative motion of the two particles and the second term is the sintering attractive force. Eq. (14) introduces the rate of approach of the two particles ($\frac{dh_{ij}}{dt}$). This value is computed using the relative velocity of the two particles (at the last time step) and the unit normal vector $\vec{n}_{ij}$.

Eq. (14) has been developed for two spherical particles (radii $r_i$ and $r_j$), with $R_{eq,ij} = r_i r_j / (r_i + r_j)$, denominated as the equivalent radius. Sintering is a process driven by the local curvature gradient. Thus, we propose to generalize Eq. (14) for non-spherical particles by ascribing the mean local curvature radius to $r_i$ and $r_j$. A similar approach was employed for the elastic Hertzian law in Podlozhnyuk et al. 2017 with good results when compared with FEM simulations. The mean local curvatures $\kappa_i$ and $\kappa_j$ are calculated from Eq. (6) on the contact points $\vec{x}_{op,\phi i}$ and $\vec{x}_{op,\phi j}$. The curvature radii are given by $r_i = 1/\kappa_i$ and $r_j = 1/\kappa_j$. The sintering neck size $a_{ij}$ in Eq. (14) is related to the indentation by Coble's geometric model Coble 1958; C. L. Martin and Bordia 2009:

$$a_{ij}^2 = 4R_{eq,ij}h_{ij} \qquad (15)$$

Tangential forces are neglected. Thus, the force between two sintering particles is given by: $\vec{F}_{ij} = N_{ij} \vec{n}_{ij}$. Unlike spheres, a normal force applied on a non-spherical particle creates a torque. When large contacts develop, it is reasonable to consider that this torque will be counteracted by a resisting torque at the contact. We have tested several conditions to stop the rotation of the two particles that pertain to a contact for which $a_{ij} / (2R_{eq,ij}) > \tilde{a}_c$. We have verified that the value of $\tilde{a}_c$ has a negligible effect in the interval [0.1, 0.5]. All simulations presented hereafter are run with $\tilde{a}_c = 0.1$. For simplicity, grain growth is not considered. Typically, grain growth becomes dominant at high relative density (> 0.85). We will restrict our simulations to the first sintering stage (< 0.80).

## 2.5 Analytical model

A simple analytical equation of the mean normalized indentation as a function of time can be derived from Eq. (14) by assuming that the packing equilibrium is attained if contacts fulfill
while the unit normal vector is calculated from

$$ \vec{n}_{ij} = \frac{\nabla \phi_i(\vec{x}_{op,\phi_i}) - \nabla \phi_j(\vec{x}_{op,\phi_j})}{|\nabla \phi_i(\vec{x}_{op,\phi_i}) - \nabla \phi_j(\vec{x}_{op,\phi_j})|} \qquad (12) $$

Following the suggestion of Li et al. (2019), our force calculation depends on the maximum indentation, instead of the average indentation, which would make it more mesh dependent.

For two spherical particles, the normal force $N_{ij}$ (Bouvard and McMeeking 1996; Martin and Bordia 2009) is

$$ N_{ij} = \frac{\pi a_{ij}^4}{8\Delta_{GB}} \frac{dh_{ij}}{dt} - \frac{\alpha}{2} \pi R_{eq,ij} \gamma_S \qquad (13) $$

where $\alpha$ is a parameter that depends on the ratio between surface and grain boundary diffusion coefficients. It is set to 2.46 for the case of a typical ceramic material (alumina) (Bouvard and McMeeking 1996; Paredes-Goyes et al. 2021). $\gamma_S$ is the surface energy and $\Delta_{GB} = \frac{\Omega}{k_b T} D_{GB} \delta_{GB}$, with the atomic volume $\Omega$, the Boltzmann constant $k_b$, the temperature $T$, the grain boundary diffusivity $D_{GB}$, the grain boundary thickness $\delta_{GB}$. $R_{eq,ij} = r_i r_j / (r_i + r_j)$ is the equivalent radius of particles i and j with radii $r_i$ and $r_j$. The first term is a viscous component that counteracts the normal component of the relative motion of the two particles and the second term is the sintering attractive force. Equation (13) introduces the rate of approach of the two particles ($dh_{ij}/dt$). This value is the dot product of the relative velocity of the two particles (at the last time step) and the unit normal vector $\vec{n}_{ij}$.

Equation (13) has been developed for two spherical particles (radii $r_i$ and $r_j$), with $R_{eq,ij} = r_i r_j / (r_i + r_j)$, denominated as the equivalent radius. Sintering is a process driven by the local curvature gradient. Thus, we propose to generalize Equation (13) for non-spherical particles by ascribing the mean local curvature radius to $r_i$ and $r_j$. A similar approach was employed for the elastic Hertzian law in (Podlozhnyuk et al. 2017) with good results when compared with FEM simulations. The mean local curvatures $\kappa_i$ and $\kappa_j$ are calculated from Equation (6) on the contact points $\vec{x}_{op,\phi_i}$ and $\vec{x}_{op,\phi_j}$. The curvature radii are given by $r_i = 1/\kappa_i$ and $r_j = 1/\kappa_j$, which are positive (as required in sintering equations) since the mean curvature values are positive for convex particles following the sign convention of Equation (6). The sintering neck radius $a_{ij}$ in Equation (13) is related to the indentation by Coble's geometric model (Coble 1958; Martin and Bordia 2009) (see (Paredes-Goyes et al. 2021, Fig. 1) for a schematic of the contact):

$$ a_{ij}^2 = 4R_{eq,ij}h_{ij}. \qquad (14) $$

Tangential forces are neglected. Our experience is that it leads to a slight overestimation of particle rearrangement and consequently of the sintering kinetics. Tangential forces are non-negligible when tackling constrained sintering (Martin and Bordia 2009). Thus, the force between two sintering particles is simply given by $\vec{F}_{ij} = N_{ij} \vec{n}_{ij}$. Unlike spheres, a normal force applied on a non-spherical particle creates a torque. When large contacts develop, it is reasonable to consider that this torque will be counteracted by a resisting torque at the contact. We have tested several conditions to stop the rotation of the two particles that pertain to a contact for which $a_{ij}/(2R_{eq,ij}) > \tilde{a}_c$. We have verified that the value of $\tilde{a}_c$ has a negligible effect in the interval [0.1, 0.5]. All simulations presented hereafter are run with $\tilde{a}_c = 0.1$. For simplicity, grain growth is not considered. Typically, grain growth becomes dominant at high relative density (> 0.85). We will restrict our simulations to the first sintering stage (< 0.80).

## 2.5 Analytical model

A simple analytical equation of the mean normalized indentation as a function of time can be derived from Equation (13) by assuming that the packing equilibrium is attained if contacts fulfill $N_{ij} = 0$. This simplistic assumption, together with assuming that all particles are spherical with the same radius $R = 2R_{eq}$, and that the initial indentation is zero leads to

$$ \left(\frac{h}{R}\right)^3 = \frac{3}{2} \alpha \frac{t}{t} \qquad (15) $$
$N_{ij} = 0$. This simplistic assumption, together with assuming that all particles are spherical with the same radius $R = 2R_{eq}$, and that the initial indentation is zero leads to:

$$ \left(\frac{h}{R}\right)^3 = \frac{3}{2} \alpha \bar{t} \qquad (16) $$

with:

$$ \tilde{t} = \frac{R^4}{\Delta_{GB}\gamma_S} \qquad (17) $$

Eq. (16) exhibits the classic dependence on time ($\frac{h}{R} \propto t^{1/3}$) proposed by Coblenz 1980 for sintering and provides a useful time normalization parameter $\tilde{t}$. The ratio $\frac{h}{R}$ is a good proxy for densification. Assuming that the indentation $h$ is zero at the initial relative density $\rho_i$ (typically smaller than the Random Close Packing) and there is no particle rearrangement, the ratio $\frac{h}{R}$ is simply related to the isostatic densification of a monomodal spherical powder Arzt 1982; Storåkers et al. 1999:

$$ \frac{h}{2R} = 1 - \left(\frac{\rho_i}{\rho}\right)^{\frac{1}{3}} \qquad (18) $$

From Eqs. (16) and (18), density evolution from the initial relative density $\rho_i$ can be related to time:

$$ \frac{\rho}{\rho_i} = \left[ 1 - \frac{1}{2} \left( \frac{3}{2} \alpha \bar{t} \right)^{\frac{1}{3}} \right]^{-3} \qquad (19) $$

Eq. (19) is valid for the first stage of sintering (up to 0.85 relative density) Arzt 1982; Storåkers et al. 1999. It has been developed for monomodal packings of spheres and thus provides a simple and useful reference for the sintering of ellipsoids.

## 3 Level-set model validation

In this and the following section, ellipsoidal particles defined by semi-axes $a, b, c$ are used for validating the Level-Set methodology. In particular, we work with different types of spheroids (ellipsoid of revolution) where $c$ is the long-axis. Both prolate ($a = b \le c$) and oblate ($a \le b = c$) are used. We have tested that a mesh with $ncell_{int} = 100$ and $ncell_{ext} = 20$, was a good balance between contact detection accuracy, CPU-time consumption and memory use. This mesh is used for all subsequent simulations.

A level-set function for these spheroids must be provided to fill the grid points as described by Eq. (5). There is no simple analytical equation for determining the level-set function of a spheroid. Numerical solutions have been proposed by Eberly David Eberly n.d. (numerical resolution to find the roots of polynomial equations) and Pope Pope 2008 (quadratic minimization problem with a constraint).

A simpler and much faster alternative consists of using an approximate of the level-set value, $\tilde{\phi}(p)$ by considering simply the distance between the center of the ellipse and the point $p$ at which $\tilde{\phi}(p)$ needs to be calculated. $p_0$ is the intersection between $\overrightarrow{Op}$ and the ellipse surface:

$$ \tilde{\phi}(p) = \|\overrightarrow{Op}\| - \|\overrightarrow{Op_0}\| \qquad (20) $$

The level-set value is approximated by the distance between point $p$ and $p_0$ where Fig. 6 shows schematically that the value is exact at the apex of the ellipsoid (and more importantly at any point on the surface of the ellipsoid) and overestimates the absolute value of $\phi(p)$ otherwise. For an ellipsoid with semi-axes $a = 0.5, b = 0.5, c = 1$, we checked that the maximum error is of the order of 10% of the larger semi-axis for $\phi(p) > 0$ with the standard grid used around the particle (see Fig. 3). The exact solution was computed using the algorithms for ellipsoids provided in Pope 2008. Thus, the use of the approximate level-set value $\tilde{\phi}(p)$ does not modify the surface of the ellipsoid for contact detection (because $\tilde{\phi}(p) = \phi(p) = 0$ on the ellipsoid surface).
with

$$ \tilde{t} = \frac{R^4}{\Delta_{GBYS}}. \qquad (16) $$

Equation (15) exhibits the classical dependence on time $h/R \propto t^{1/3}$ proposed by Coblenz et al. (1980) for sintering and provides a useful time normalization parameter $\tilde{t}$. The ratio $h/R$ is a good proxy for densification. Assuming that the indentation $h$ is zero at the initial relative density $\rho_i$ (typically smaller than the Random Close Packing) and there is no particle rearrangement, the ratio $h/R$ is simply related to the isostatic densification of a monomodal spherical powder (Arzt 1982; Storåkers et al. 1999):

$$ \frac{h}{2R} = 1 - \left(\frac{\rho_i}{\rho}\right)^{1/3}. \qquad (17) $$

From Equations (15) and (17), density evolution from the initial relative density $\rho_i$ can be related to time:

$$ \frac{\rho}{\rho_i} = \left[ 1 - \frac{1}{2} \left( \frac{3}{2} \alpha \bar{z} \frac{t}{\tilde{t}} \right)^{1/3} \right]^{-3}. \qquad (18) $$

Equation (18) is valid for the first stage of sintering (up to 0.85 relative density) (Arzt 1982; Storåkers et al. 1999). It has been developed for monomodal packings of spheres and thus provides a simple and useful reference for the sintering of ellipsoids.

## 3 Level-set model validation

In this and the following section, ellipsoidal particles defined by semi-axes $a, b, c$ are used for validating the Level-Set methodology. In particular, we work with different types of spheroids (ellipsoid of revolution) where $c$ is the long-axis. Both prolate ($a = b \le c$) and oblate ($a \le b = c$) are used. We have tested that a mesh with $ncell_{int} = 100$ and $ncell_{ext} = 20$, was a good balance between contact detection accuracy, CPU-time consumption and memory use. This mesh is used for all subsequent simulations.

A level-set function for these spheroids must be provided to fill the grid points as described by Equation (5). There is no simple analytical equation for determining the level-set function of a spheroid. Numerical solutions have been proposed by Eberly (2020) (numerical resolution to find the roots of polynomial equations) and Pope (2008) (quadratic minimization problem with a constraint).

A simpler and much faster alternative consists of using an approximate of the level-set value, $\tilde{\phi}(p)$ by considering simply the distance between the center of the ellipse and the point $p$ at which $\tilde{\phi}(p)$ needs to be calculated. $p_0$ is the intersection between $\overrightarrow{Op}$ and the ellipse surface:

$$ \tilde{\phi}(p) = \|\overrightarrow{Op}\| - \|\overrightarrow{Op_0}\|. \qquad (19) $$

The level-set value is approximated by the distance between point $p$ and $p_0$ where Figure 7 shows that the value is exact at the apex of the ellipsoid (and more importantly at any point on the surface of the ellipsoid) and overestimates the absolute value of $\phi(p)$ otherwise. For an ellipsoid with semi-axes $a = 0.5, b = 0.5, c = 1$, we checked that the maximum error is of the order of 10% of the larger semi-axis for $\phi(p) > 0$ with the standard grid used around the particle, see Figure 3. The exact solution was computed using the algorithms for ellipsoids provided in (Pope 2008). Thus, the use of the approximate level-set value $\tilde{\phi}(p)$ does not modify the surface of the ellipsoid for contact detection (because $\tilde{\phi}(p) = \phi(p) = 0$ on the ellipsoid surface). However, once detected the approximate level-set contact slightly overestimates the indentation (positive error), making the contact stiffer as compared to the exact level-set. For example, Figure 7(c) shows that for an elongated prolate ellipsoid ($a = 0.5, b = 0.5, c = 1$), the error is smaller than 0.1% for small contacts (typically elastic contacts) and smaller than 5% for very large indentation ($\phi(p) \approx -0.3$) that arise at the end of sintering.
However, once detected the approximate level-set contact slightly overestimates the indentation (positive error), making the contact stiffer as compared to the exact level-set. For example, Fig. 6c shows that for an elongated prolate ellipsoid ($a = 0.5, b = 0.5, c = 1$), the error is smaller than 0.1% for small contacts (typically elastic contacts) and smaller than 5% for very large indentation ($\phi(p) \approx -0.3$) that arise at the end of sintering.

<figure><img src="image_8.png" /><figcaption>Figure 6: Schematic of exact, $\phi(p)$, and approximate $\tilde{\phi}(p)$ level-sets for ellipsoids. (a) for a point $p$ outside and (b) inside the ellipsoid. Relative error of the approximate solution $\phi(p)$ for an elongated prolate ellipsoid ($a = 0.5, b = 0.5, c = 1$).</figcaption></figure>

## 3.1 Elastic contact between two ellipsoidal particles

The first test aims at validating our model against other DEM formulation and FEM solution Podlozhnyuk et al. 2017; Zheng et al. 2013. For this, a simple system of two equal-size ellipsoids (Fig. 7) aligned on the z axis with elastic contact is chosen. The contact force is calculated as a function of the mutual indentation using the Hertz equation:

$$N_{Hertz} = \frac{2}{3} \frac{E}{1 - v^2} R_{eq}^{1/2} h^{3/2} \qquad (21)$$

where E is the Young's modulus and v the Poisson's coefficient.

<figure><img src="image_9.png" /><figcaption>Figure 7: Resulting normalized elastic force versus the indentation between two aligned ellipsoids for LS-DEM, DEM superquadrics Podlozhnyuk et al. 2017 and FEM Zheng et al. 2013. The force is normalized by $Eb^2$, where E is the Young's modulus and b is the short axis.</figcaption></figure>

A complete description of the problem is found in Podlozhnyuk et al. 2017. The two ellipsoids are identical with $a = 5, b = 2.5, c = 2.5$ in mm. Fig. 7 shows the results of our LS-DEM calculation
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7** Exact, $\phi(p)$, and approximate $\tilde{\phi}(p)$ level-sets for ellipsoids: (a) for a point $p$ outside and (b) inside the ellipsoid. (c) Relative error of the approximate solution $\tilde{\phi}(p)$ for an elongated prolate ellipsoid ($a = 0.5, b = 0.5, c = 1$).</figcaption></figure>

### 3.1 Elastic contact between two ellipsoidal particles

The first test aims at validating our model against other DEM formulation and FEM solution (Podlozhnyuk et al. 2017; Zheng et al. 2013). For this, a simple system of two equal-size ellipsoids, Figure 8, aligned on the z axis with elastic contact is chosen. The contact force is calculated as a function of the mutual indentation using the Hertz equation

$$N_{\text{Hertz}} = \frac{2}{3} \frac{E}{1 - \nu^2} R_{\text{eq}}^{1/2} h^{3/2} \qquad (20)$$

where E is Young's modulus and $\nu$, Poisson's coefficient.

A complete description of the problem is found in (Podlozhnyuk et al. 2017). The two ellipsoids are identical with $a = 5$, $b = 2.5$, $c = 2.5$ in mm. Figure 8 shows the results of our LS-DEM calculation compared with DEM superquadrics (Podlozhnyuk et al. 2017) and FEM calculation (Zheng et al. 2013). The overlay of the DEM superquadrics and our curve

<figure><img src="image_10.png" /><figcaption>**Figure 8** Resulting normalized elastic force versus the indentation between two aligned ellipsoids for LS-DEM, DEM superquadrics (Podlozhnyuk et al. 2017) and FEM (Zheng et al. 2013). The force is normalized by $Eb^2$, where E is the Young modulus and b, the short axis.</figcaption></figure>

demonstrates that the calculation of the local curvature and the first and second derivatives of our Level Set function is correct. The close results to FEM values support the previously discussed hypothesis of using the local curvature radii in the Hertzian law. Additionally, it is possible to confirm the ability to detect very small contacts, which will be discussed further in the following sections.

### 3.2 Jamming of a packing of ellipsoidal particles

The relative density of packings prior to sintering is usually between 0.5-0.6. This is close but smaller than the density of a random close packing or a maximally random jammed state (Torquato et al. 2000). This packing can be obtained numerically by jamming an initial random gas of
compared with DEM superquadrics Podlozhnyuk et al. 2017 and FEM calculation Zheng et al. 2013. The overlay of the DEM superquadrics and our curve demonstrates that the calculation of the local curvature and the first and second derivatives of our Level Set function is correct. The close results to FEM values support the previously discussed hypothesis of using the local curvature radii in the Hertzian law. Additionally, it is possible to confirm the ability to detect very small contacts, which will be discussed further in the following sections.

## 3.2 Jamming of a packing of ellipsoidal particles

The relative density of packings prior to sintering is usually between 0.50-0.6. This is close but smaller than the density of a random close packing or a maximally random jammed state Torquato et al. 2000. This packing can be obtained numerically by jamming an initial random gas of particles. The procedure consists of randomly locating particles into the simulation box with the constraint that there is no contact between them. The jamming itself is a stress-controlled simulation with a very small control pressure *P* as compared to the Young's modulus of the particles (*P*/*E* = 10⁻⁶). This ensures that mutual indentation between particles is kept very small and does not contribute to densification (only particle rearrangement contributes to densification). With such a scheme, the maximum packing fraction is asymptotically approached, while the pressure tends toward the control pressure. Details on this process for spherical particles can be found in C. L. Martin and Bordia 2008.

The jamming process is a relevant test of the proposed algorithm for contact detection, as it tests its ability to detect small contacts. We first tested the ability of LS-DEM of detecting all contacts in a jamming procedure as compared to standard DEM for spheres. A periodic box with 1 000 spheres was used. For such a packing, the contacts exhibit a small normalized indentation (*h*/*R* ≤ 4 10⁻⁴) due to the very small control pressure. Fig. 8a shows that LS-DEM accounts for the same coordination number evolution as the standard DEM, demonstrating its ability to account for all contacts.

For ellipsoids, the packing procedure is performed for two packings of 1 000 non-uniform size particles composed by prolate ellipsoids with aspect ratio *c*/*a* = *c*/*b* = 1/0.6 and oblate ellipsoids with the same aspect ratio *c*/*a* = *b*/*a* = 1/0.6. The simulation box is periodic. A bounded box with stiff planes that interact with particles (Eq.(21)) was also tested with the same qualitative conclusions. Fig. 8c shows the jamming of prolate and oblate ellipsoids starting from a relative density of 0.22 and 0.36, respectively, to a final density of 0.60. The normalized mean indentation shows a base line of the order of *h*/*R* = 10⁻⁴ with spikes around 5.10⁻³ (Fig. 8b).

<figure><img src="image_10.png" /><figcaption>Figure 8: (a) Coordination number during jamming of spheres. (b) Mean indentation of ellipsoids as a function of relative density, *ρ*, during jamming. (c) Initial and final microstructures obtained.</figcaption></figure>

## 3.3 Sintering of a packing of spherical particles

The tests with elastic contacts in the preceding sections have validated the contact detection scheme. The objective of this section is to validate the simulations that implement sintering contacts. A packing of LS spherical particles is chosen to compare with the analytical solution and DEM results. The initial packing contains 1 000 particles and exhibits a relative density *ρ*ᵢ = 0.6, that is, a packing with very small indentations. The simulation is carried out until *h*/*R*=0.1,
particles. The procedure consists of randomly locating particles into the simulation box with the constraint that there is no contact between them. The jamming itself is a stress-controlled simulation with a very small control pressure *P* as compared to the Young's modulus of the particles (*P*/*E* = 10⁻⁶). This ensures that mutual indentation between particles is kept very small and does not contribute to densification (only particle rearrangement contributes to densification). With such a scheme, the maximum packing fraction is asymptotically approached, while the pressure tends toward the control pressure. Details on this process for spherical particles can be found in (Martin and Bordia 2008).

The jamming process is a relevant test of the proposed algorithm for contact detection, as it tests its ability to detect small contacts. We first tested the ability of LS-DEM of detecting all contacts in a jamming procedure as compared to standard DEM for spheres. A periodic box with 1000 spheres was used. For such a packing, the contacts exhibit a small normalized indentation (*h*/*R* ≤ 4·10⁻⁴) due to the very small control pressure. Figure 9(a) shows that LS-DEM accounts for the same coordination number evolution as the standard DEM, demonstrating its ability to account for all contacts.

For ellipsoids, the packing procedure is performed for two packings of 1000 non-uniform size particles composed by prolate ellipsoids with aspect ratio *c*/*a* = *c*/*b* = 1/0.6 and oblate ellipsoids with the same aspect ratio *c*/*a* = *b*/*a* = 1/0.6. The simulation box is periodic. A bounded box with stiff planes that interact with particles, see Equation (20), was also tested with the same qualitative conclusions. Figure 9(c) shows the jamming of prolate and oblate ellipsoids starting from a relative density of 0.22 and 0.36, respectively, to a final density of 0.6. In Figure 9(b), the normalized mean indentation shows a base line of the order of *h*/*R* = 10⁻⁴ with spikes around 5 · 10⁻³.

<figure><img src="image_11.png" /><figcaption>**Figure 9** Sintering kinetics of ellipsoids: (a) Coordination number during jamming of spheres. (b) Mean indentation of ellipsoids as a function of relative density during jamming. (c) Initial and final microstructures obtained.</figcaption></figure>

### 3.3 Sintering of a packing of spherical particles

The tests with elastic contacts in the preceding sections have validated the contact detection scheme. The objective of this section is to validate the simulations that implement sintering contacts. A packing of LS spherical particles is chosen to compare with the analytical solution and DEM results. The initial packing contains 1000 particles and exhibits a relative density *ρ*ᵢ = 0.6, that is, a packing with very small indentations. The simulation is carried out until *h*/*R* = 0.1, which is the most critical stage for contact detection and calculation.

The left axis of Figure 10(a) shows the evolution of the mean indentation with normalized time, see Equation (16), for standard DEM and LS-DEM simulations. Both curves show the same evolution, with LS-DEM leading to slightly larger values of *h*/*R*. The right axis of Figure 10(a) plots (*h*/*R*)³ vs normalized time. Both the LS-DEM and DEM simulations lead to a linear relation, as predicted by the analytical solution in Equation (15) which leads to larger indentation as compared to the DEM and LS-DEM. This is because it imposes isostatic densification of particles simply by imposing force equilibration at the contact level instead of force equilibration at the level of particles for DEM.

Figure 10(b) shows the evolution of the mean coordination number with time for the standard
which is the most critical stage for contact detection and calculation.

The left axis of Fig. 9a shows the evolution of the mean indentation with normalized time (Eq. (17)) for standard DEM and LS-DEM simulations. Both curves show the same evolution, with LS-DEM leading to slightly larger values of $h/R$. The right axis of Fig. 9a plots $(h/R)^3$ vs normalized time. Both the LS-DEM and DEM simulations lead to a linear relation, as predicted by the analytical solution (Eq.(16)). Eq. (16) leads to larger indentation as compared to the DEM and LS-DEM. This is because it imposes isostatic densification of particles simply by imposing force equilibration at the contact level instead of force equilibration at the level of particles for DEM.

<figure><img src="image_11.png" /><figcaption>Figure 9: Sintering of spheres with standard DEM and LS-DEM. (a) Evolution of the mean normalized indentation ($h/R$) with normalized time ($t/\tilde{t}$, Eq.(17)). Comparison with the analytical solution (Eq. (16)). (b) Evolution of the average coordination number with normalized time.</figcaption></figure>

Fig. 9b shows the evolution of the mean coordination number with time for the standard DEM and the LS-DEM. Good agreement is also found, although the LS-DEM leads to a slightly larger number of contacts at the end of the sintering. The same agreement was observed for the neck size $a$. Therefore, a good correlation was found between LS-DEM and DEM simulations of spheres, both for the densification kinetics ($h$) and the consolidation kinetics ($a$).

# 4 Sintering applications

## 4.1 Sintering between two ellipsoidal particles

To better understand the sintering behavior of ellipsoidal particles, a simple system of two equal-sized particles is examined. Two prolate ellipsoids in contact along their long axis or along their short axis are studied to evaluate the influence of curvature on sintering rates. The aspect ratio is chosen to be small ($a = b = 0.5 \times c$) to depart significantly from the sphere and accentuate the effects of curvature. Simulations are also performed with LS spheres ($a = b = c$) to assess the influence of the aspect ratio of the particles. Comparing ellipsoids of different aspect ratios (here between elongated ellipsoids and spheres) in the context of sintering requires a certain amount of caution. Because the sintering characteristic time scales with $(1/R^4)$ (Eq. (16)), the particle size plays a dominant role. We have therefore chosen to compare ellipsoids of different aspect ratios while maintaining the same particle volume, in line with other works studying the effect of morphology in powder sintering. Yan et al. 2017; Bjørk 2022. In other words, in all subsequent simulations, we define a normalizing radius, $\tilde{R}$, that fulfills:

$$ \frac{4}{3}\pi\tilde{R}^3 = \frac{4}{3}\pi abc \qquad (22) $$

Fig. 10a shows the evolution of the normalized indentation, $\frac{h}{R}$. As in Fig. 9, Eq. (16) provides a correct approximation of the sintering kinetics for two spheres.

Although simplistic, the sintering of two elongated ellipsoids offers important insight into the effect of departure from sphericity. The vertical ($V$) configuration (see Fig. 10b) densifies faster
<figure><img src="image_12.png" /><figcaption>**Figure 10** Sintering of spheres with standard DEM and LS-DEM. (a) Mean normalized indentation $h/R$ versus normalized time $t/\tilde{t}$, see Equation (16). Comparison with the analytical solution, Equation (15). (b) Average coordination number versus normalized time. The finer grid ($n_{cell, int} = 150$, $n_{cell, ext} = 40$) shows essentially the same result.</figcaption></figure>

DEM and the LS-DEM. Good agreement is also found, although the LS-DEM leads to a slightly larger number of contacts at the end of the sintering. In addition, we performed five simulations with five different random seeds for generating the initial packings. Figure 10(b) shows the worst comparison of the five tests. No clear trend was found in the five tests with regard to over- or underestimation of the number of contacts by the LS DEM compared with the standard DEM. Also, we tested a finer grid in Figure 3, but Figure 10(b) shows that the number of contacts is not affected. The same type of agreement was observed for the neck size $a$. Therefore, a good correlation was found between LS-DEM and DEM simulations of spheres, both for the densification kinetics $h$ and the consolidation kinetics $a$.

# 4 Sintering applications

## 4.1 Sintering between two ellipsoidal particles

To better understand the sintering behavior of ellipsoidal particles, a simple system of two equal-sized particles is examined. Two prolate ellipsoids in contact along their long axis or along their short axis are studied to evaluate the influence of curvature on sintering rates. The aspect ratio is chosen to be small ($a = b = 0.5 \times c$) to depart significantly from the sphere and accentuate the effects of curvature. Simulations are also performed with LS spheres ($a = b = c$) to assess the influence of the aspect ratio of the particles. Comparing ellipsoids of different aspect ratios (here between elongated ellipsoids and spheres) in the context of sintering requires a certain amount of caution. Because the sintering characteristic time scales with $(1/R^4)$, see Equation (15), the particle size plays a dominant role. We have therefore chosen to compare ellipsoids of different aspect ratios while maintaining the same particle volume, in line with other works studying the effect of morphology in powder sintering (Yan et al. 2017; Bjørk 2022). In other words, in all subsequent simulations, we define a normalizing radius, $\tilde{R}$, that fulfills

$$ \frac{4}{3}\pi\tilde{R}^3 = \frac{4}{3}\pi abc. \qquad (21) $$

Figure 11(a) shows the evolution of the normalized indentation, $h/\tilde{R}$. As in Figure 10, Equation (15) provides a correct approximation of the sintering kinetics for two spheres.

Although simplistic, the sintering of two elongated ellipsoids offers important insight into the effect of departure from sphericity. The vertical ($V$) configuration, see Figure 11(b), densifies faster than the spherical configuration, which itself densifies slightly faster than the horizontal ($H$) configuration. Figure 11(a) shows that the cubic relation between densification and time
<figure><img src="image_12.png" /><figcaption>Figure 10: (a) Densification and (b) consolidation kinetics of two prolate ellipsoids ($a = b = 0.5 \times c$) aligned along their vertical long axis (V) or horizontal long axis (H). Comparison with two spheres and with Eq. (16). All particles have the same volume.</figcaption></figure>

than the spherical configuration, which itself densifies slightly faster than the horizontal (H) configuration. Fig. 10a shows that the cubic relation between densification and time remains valid for ellipsoids. This is because Eq. (14) uses the equivalent radius $R_{eq}$ that stays approximately constant all along the simulation (no particle rotation). The equivalent radius dictates the contact size (Eq. (15)):

$$ \frac{a}{\tilde{R}} = \sqrt{2} \sqrt{\frac{2R_{eq}}{\tilde{R}}} \sqrt{\frac{h}{\tilde{R}}} \qquad (23) $$

and differs markedly for the three configurations:

$$ V, 2R_{eq}/\tilde{R} \approx 0.4 - S, 2R_{eq}/\tilde{R} \approx 1.0 - H, 2R_{eq}/\tilde{R} \approx 1.3. $$

A smaller contact radius $a$ leads to a lower viscous term in Eq. (14), which scales with $a^4$. Fig. 10b shows the evolution of contact size, with an inverse hierarchy to Fig. 10a, as confirmed by Eq. (23).

## 4.2 Sintering of a packing of ellipsoidal particles

The previous section, limited to two ellipsoidal particles, showed that densification (indentation) and consolidation (contact size) deviate from the case of spherical particles. In particular, it demonstrated that the local curvature at the contact plays a central role in the densification kinetics. Simulations of a particle packing can provide more realistic information. Several packings with aspect ratio in the range of 0.50 to 1.0 are sintered with periodic conditions, starting from a 0.60 relative density. Each packing contains 1 000 prolate or oblate ellipsoidal particles. These packings are slightly non-monomodal (±5% deviation on the bounding sphere radius). Oblate, prolate, and spherical particles have all the same volume (Eq. (22)).

Fig. 11 shows the normalized time to sinter each packing from density 0.60 to 0.70. Some simulations were conducted up to larger densities but the qualitative conclusions are the same as those drawn in Fig. 11. For each simulation, at least three different random seeds are chosen to generate different initial randomly packed configurations, thus providing a rough estimate of dispersion. Fig. 11 shows that the oblate and prolate ellipsoids behave similarly, with the oblate particles sintering slightly slower for the elongated ellipsoids. In the range of aspect ratio tested here, [0.5-1], our simulations indicate that it is beneficial to depart from spheres. More precisely, Fig. 11 shows the existence of an optimum aspect ratio (around 0.8-0.85), which minimizes the sintering time.

This result can be understood by recalling the role of the curvature radius demonstrated in the preceding section and by observing the evolution of the distribution of equivalent radii $R_{eq}$ with the aspect ratio. Fig. 12 shows this evolution for the normalized value $2R_{eq}/\tilde{R}$ for three
<figure><img src="image_13.png" /><figcaption>**Figure 11** Jamming: (a) Densification and (b) consolidation kinetics of two prolate ellipsoids ($a = b = 0.5 \times c$) aligned along their vertical long axis (V) or horizontal long axis (H). Comparison with two spheres and with Equation (15). All particles have the same volume.</figcaption></figure>

remains valid for ellipsoids. This is because Equation (13) uses the equivalent radius $R_{eq}$ that stays approximately constant all along the simulation (no particle rotation). The equivalent radius dictates the contact radius, see Equation (14):

$$ \frac{a}{\tilde{R}} = \sqrt{2} \sqrt{\frac{2R_{eq}}{\tilde{R}}} \sqrt{\frac{h}{\tilde{R}}} \qquad (22) $$

and differs markedly for the three configurations: $V, 2R_{eq}/\tilde{R} \approx 0.4 - S, 2R_{eq}/\tilde{R} \approx 1.0 - H, 2R_{eq}/\tilde{R} \approx 1.3$. A smaller contact radius $a$ leads to a lower viscous term in Equation (13), which scales with $a^4$. Figure 11(b) shows the evolution of contact radius, with an inverse hierarchy to Figure 11(a), as confirmed by Equation (22).

## 4.2 Sintering of a packing of ellipsoidal particles

The previous section, limited to two ellipsoidal particles, showed that densification (indentation) and consolidation (contact radius) deviate from the case of spherical particles. In particular, it demonstrated that the local curvature at the contact plays a central role in the densification kinetics. Simulations of a particle packing can provide more realistic information. Several packings with aspect ratio in the range of 0.5 to 1.0 are sintered with periodic conditions, starting from a 0.6 relative density. Each packing contains 1000 prolate or oblate ellipsoidal particles. These packings are slightly non-monomodal (±5% deviation on the bounding sphere radius). Oblate, prolate, and spherical particles have all the same volume, see Equation (21).

Figure 12 shows the normalized time to sinter each packing from density 0.6 to 0.7. Some simulations were conducted up to larger densities but the qualitative conclusions are the same as those drawn in Figure 12. For each simulation, at least three different random seeds are chosen to generate different initial randomly packed configurations, thus providing a rough estimate of dispersion. Figure 12 shows that the oblate and prolate ellipsoids behave similarly, with the oblate particles sintering slightly slower for the elongated ellipsoids. In the range of aspect ratio tested here, [0.5 - 1], our simulations indicate that it is beneficial to depart from spheres. More precisely, Figure 12 shows the existence of an optimum aspect ratio (around 0.8-0.85), which minimizes the sintering time.

This result can be understood by recalling the role of the curvature radius demonstrated in the preceding section and by observing the evolution of the distribution of equivalent radii $R_{eq}$ with the aspect ratio. Figure 13 shows this evolution for the normalized value $2R_{eq}/\tilde{R}$ for three packings: spheres and two oblate configurations with aspect ratios of 0.8 and 0.5. All distributions lead to an average value around $2R_{eq}/\tilde{R} = 1$. However, the average value does not provide a
<figure><img src="image_13.png" /><figcaption>**Figure 11:** Sintering of packings of ellipsoids with different aspect ratios: normalized time (Eq. (17)) to reach the same packing density (0.7). Both prolate and oblate ellipsoids are tested.</figcaption></figure>

packings: spheres and two oblate configurations with aspect ratios of 0.8 and 0.5. All distributions lead to an average value around $2R_{eq}/\tilde{R} = 1$. However, the average value does not provide a complete picture. The distribution for spheres is very narrow, since only the ±5% size distribution leads to some deviation from unity. For non-spherical particles, the distribution widens as the aspect ratio increases. For the 0.8 aspect ratio, the distribution is nearly symmetric with small values of $2R_{eq}/\tilde{R}$, which are beneficial for sintering kinetics, and large values of $2R_{eq}/\tilde{R}$, which are detrimental for sintering kinetics. On the balance, this distribution is beneficial, since small $2R_{eq}/\tilde{R}$ values lead to faster densification kinetics, while moderately large $2R_{eq}/\tilde{R}$ values hardly hinder densification. For an aspect ratio of 0.5, the distribution is no longer symmetric, with very high $2R_{eq}/\tilde{R}$ values severely hampering densification kinetics. The detrimental effect of large $R_{eq}$ is confirmed by Eq. (16), which shows that $h/R$ is not inversely linear to $R$ but scales with $1/R^{4/3}$.

The optimum aspect ratio at 0.85 can also be attributable to the additional degrees of freedom that departure from the spherical shape allows. The initial density before sintering is 0.6 in our simulations to allow for a fair comparison between all packings. This density is smaller than the Random Close Packing (RCP) of spheres (0.63). Furthermore, Donev et al. 2004 have demonstrated that ellipsoids can randomly pack more densely than spheres with an RCP of around 0.70 for an optimal aspect ratio of 0.65 for both oblate and prolate ellipsoids. Although their optimal aspect ratio is smaller than ours, this would indicate that packings with aspect ratios around 0.85 have additional possibilities to rearrange as compared to packings of spheres, thus promoting densification.

Yan et al. Yan et al. 2017 and Bjørk 2022 have already suggested that, for the same volume, elongated particles sinter faster than spherical ones. Both studies used Monte Carlo simulations and included grain growth. They did not observe an optimal aspect ratio, as highly elongated particles (up to aspect ratios of 0.5 for oblate and prolate ellipsoids) continue to exhibit shorter sintering times.

# 5 Concluding remarks

The modeling of sintering of non-spherical particles is challenging. The sintering LS-DEM model developed in this work looks promising to study the densification of realistic packings composed of non-spherical particles. The framework developed can be applied to sintering of arbitrary shaped particles. Here, the sintering of ellipsoidal particles has been simulated as a proof of concept. The natural next step would be to use particle shapes derived from 3D imaging of real particles (X-ray tomography).
<figure><img src="image_14.png" /><figcaption>**Figure 12** Sintering of packings of ellipsoids with different aspect ratios: normalized time, Equation (16), to reach the same packing density (0.7). Both prolate and oblate ellipsoids are tested.</figcaption></figure>

complete picture. The distribution for spheres is very narrow, since only the ±5 % size distribution leads to some deviation from unity. For non-spherical particles, the distribution widens as the aspect ratio increases. For the 0.8 aspect ratio, the distribution is nearly symmetric with small values of $2R_{eq}/\tilde{R}$, which are beneficial for sintering kinetics, and large values of $2R_{eq}/\tilde{R}$, which are detrimental for sintering kinetics. On the balance, this distribution is beneficial, since small $2R_{eq}/\tilde{R}$ values lead to faster densification kinetics, while moderately large $2R_{eq}/\tilde{R}$ values hardly hinder densification. For an aspect ratio of 0.5, the distribution is no longer symmetric, with very high $2R_{eq}/\tilde{R}$ values severely hampering densification kinetics. The detrimental effect of large $R_{eq}$ is confirmed by Equation (15), which shows that $h/R$ is not inversely linear to $R$ but scales with $1/R^{4/3}$.

The optimum aspect ratio at 0.85 can also be attributable to the additional degrees of freedom that departure from the spherical shape allows. The initial density before sintering is 0.6 in our simulations to allow for a fair comparison between all packings. This density is smaller than the Random Close Packing (RCP) of spheres (0.63). Furthermore, Donev et al. (2004) have demonstrated that ellipsoids can randomly pack more densely than spheres with an RCP of around 0.70 for an optimal aspect ratio of 0.65 for both oblate and prolate ellipsoids. Although their optimal aspect ratio is smaller than ours, this would indicate that packings with aspect ratios around 0.85 have additional possibilities to rearrange as compared to packings of spheres, thus promoting densification.

Yan et al. (2017); Bjørk (2022) have already suggested that, for the same volume, elongated particles sinter faster than spherical ones. Both studies used Monte Carlo simulations and included grain growth. They did not observe an optimal aspect ratio, as highly elongated particles (up to aspect ratios of 0.5 for oblate and prolate ellipsoids) continue to exhibit shorter sintering times.

# 5 Concluding remarks

The modeling of sintering of non-spherical particles is challenging. The sintering LS-DEM model developed in this work looks promising to study the densification of realistic packings composed of non-spherical particles. The framework developed can be applied to sintering of arbitrary shaped particles. Here, the sintering of ellipsoidal particles has been simulated as a proof of concept. The natural next step would be to use particle shapes derived from 3D imaging of real particles (X-ray tomography).

Contact detection is one of the most critical stages when modeling non-spherical particles. The presented optimization-based contact detection is effective for detecting small contacts, that are important for sintering due to the high attractive forces that come with small contacts. The
<figure><img src="image_14.png" /><figcaption>**Figure 12:** Distribution of the equivalent radius (normalized by $\frac{1}{2}\tilde{R}$) for spheres and oblate ellipsoids (aspect ratios 0.8 and 0.5). All plots contain 100 bins. Typical contact configurations for small and large values of $2R_{eq}/\tilde{R}$ are shown for ellipsoids.</figcaption></figure>

Contact detection is one of the most critical stages when modeling non-spherical particles. The presented optimization-based contact detection is effective for detecting small contacts, that are important for sintering due to the high attractive forces that come with small contacts. The use of a multi-start optimization algorithm was proposed in order to decrease the number of missed contacts. A similar technique can be used in future work to detect multiple contact points in non-convex particles.

Our optimization-based contact detection scheme consumes less CPU-time than the surface node technique, as the number of evaluations of the Level-Set function is lower. However, our simulations show that the number of grid cells and therefore the RAM consumption is higher in the first case. A detailed comparison of the overall computational cost between both methods is needed. If RAM consumption is a limiting factor, one solution may be to switch to more efficient grid structures, such as the octree proposed in Duriez et al. 2021. Another possible approach can be to use a density-based grid cell size, as a fine grid is mainly needed at low densities when contacts are very small.

The LS-DEM model described here was fully integrated into the dp3D code C. L. Martin 2024. Thus, all the contact laws that exist in dp3D should a priori be available for use with LS-DEM. The dp3D code is parallelized with a fine-grained method (at the loop level) with openMP directives. The LS-DEM contact detection scheme is within the main parallelized loop over contacts. However, it is clear that some optimization is needed. In particular, the access to the memory is an important bottleneck when parallelizing a code and it is an added issue for the current version of the LS-DEM, which uses many very large arrays to access grid points.

Moreover, here we extend the normal force and neck size expressions of spherical to non-spherical particles by simply using the local curvature radii instead of the sphere radii. A validation of this assumption is needed. Alternatively, a new formulation of an expression of the sintering force for non-spherical particles should be proposed. A formulation in terms of the overlap volume can be envisaged by using the recently developed VLS-DEM methodology van der Haven et al. 2023. This was already proposed for elasticity with ellipsoidal particles Zheng et al. 2013 or arbitrary shaped particles Feng et al. 2012; Feng 2021.
<figure><img src="image_15.png" /><figcaption>**Figure 13** Distribution of the equivalent radius (normalized by $\tilde{R}/2$) for spheres and oblate ellipsoids (aspect ratios 0.8 and 0.5). All plots contain 100 bins. Typical contact configurations for small and large values of $2R_{eq}/\tilde{R}$ are shown for ellipsoids.</figcaption></figure>

use of a multi-start optimization algorithm was proposed in order to decrease the number of missed contacts. A similar technique can be used in future work to detect multiple contact points in non-convex particles.

Our optimization-based contact detection scheme consumes less CPU-time than the surface node technique, as the number of evaluations of the Level-Set function is lower. However, our simulations show that the number of grid cells and therefore the RAM consumption is higher in the first case. A detailed comparison of the overall computational cost between both methods is needed. If RAM consumption is a limiting factor, one solution may be to switch to more efficient grid structures, such as the octree proposed in (Duriez and Bonelli 2021). Another possible approach can be to use a density-based grid cell size, as a fine grid is mainly needed at low densities when contacts are very small.

The LS-DEM model described here was fully integrated into the dp3D code (Martin 2024). Thus, all the contact laws that exist in dp3D should a priori be available for use with LS-DEM. The dp3D code is parallelized with a fine-grained method (at the loop level) with openMP directives. The code is available upon request. The LS-DEM contact detection scheme is within the main parallelized loop over contacts. However, it is clear that some optimization is needed. In particular, the access to the memory is an important bottleneck when parallelizing a code and it is an added issue for the current version of the LS-DEM, which uses many very large arrays to access grid points.

Moreover, here we extend the normal force and neck size expressions of spherical to non-spherical particles by simply using the local curvature radii instead of the sphere radii. A validation of this assumption is needed. Alternatively, a new formulation of an expression of the sintering force for non-spherical particles should be proposed. A formulation in terms of the overlap volume can be envisaged by using the recently developed VLS-DEM methodology (Haven et al. 2023). This was already proposed for elasticity with ellipsoidal particles (Zheng et al. 2013) or arbitrary shaped particles (Feng et al. 2012; Feng 2021).

# 6 References

Arzt, E. (1982). The influence of an increasing particle coordination on the densification of spherical powders. *Acta Metallurgica* 30(10):1883–1890. [DOI], [OA].

Besler, R., M. Rossetti da Silva, J. J. do Rosario, M. Dosta, S. Heinrich, and R. Janssen (2015). Sintering simulation of periodic macro porous alumina. *Journal of the American Ceramic Society* 98(11):3496–3502. [DOI], [OA].

Bjørk, R. (2022). The sintering behavior of ellipsoidal particles. *Journal of the American Ceramic Society* 105(10):6427–6436. [DOI], [OA].
# 6 Backmatter information

## Funding

* This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement MATHEGRAM No 813202.

## Acknowledgements

* To Aatreya Manjulagiri Venkatesh for providing the SEM images in Fig. 1 from his PhD thesis.

## Competing interests

* The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

* The authors declare that they have no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a full copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## References

Arzt, E. (1982). "The influence of an increasing particle coordination on the densification of spherical powders". *Acta Met.* 30, pp. 1883-1890

Besler, R., M. Rossetti Da Silva, J. J. Do Rosario, M. Dosta, S. Heinrich, and R. Janssen (2015). "Sintering Simulation of Periodic Macro Porous Alumina". *Journal of the American Ceramic Society* 98.11, pp. 3496-3502. DOI: 10.1111/jace.13684

Bjørk, R. (2022). "The sintering behavior of ellipsoidal particles". *Journal of the American Ceramic Society* May, pp. 1-10. DOI: 10.1111/jace.18580

Bouvard, D. and R. M. McMeeking (1996). "Deformation of interparticle necks by diffusion controlled creep". *J. Am. Ceram. Soc.* 79, pp. 666-672

Bustamante, D., A. X. Jerves, and S. A. Pazmiño (2020). "A generalized three-dimensional discrete element method with electrostatic induced cohesion". *Granular Matter* 22.4, pp. 1-13. DOI: 10.1007/s10035-020-01048-4. URL: https://doi.org/10.1007/s10035-020-01048-4

Christoffersen, J., M. M. Mehrabadi, and S. Nemat-Nasser (1981). "A Micromechanical Description of Granular Material Behavior." 48.June, p. 67

Coble, R. L. (1958). "Initial Sintering of Alumina and Hematite". *Journal of the American Ceramic Society* 41.2, pp. 55-62

Coblenz William; Coble, R. (1980). "Initial stage solid state sintering models. A critical analysis and assessment". *Sintering Processes* 13.January 1980, pp. 141-157. DOI: 10.1007/978-1-4899-5301-8

David Eberly (n.d.). *Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid*. URL: https://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Davis, A., B. West, N. Frisch, D. O'Connor, and M. Parno (2021). "ParticLS: Object-oriented software for discrete element methods and peridynamics". *Computational Particle Mechanics*. DOI: 10.1007/s40571-021-00392-3. URL: https://doi.org/10.1007/s40571-021-00392-3

Dervieux, A. and F. Thomasset (1980). "A finite element method for the simulation of a Rayleigh-Taylor instability". *Approximation methods for Navier-Stokes problems*. Springer, pp. 145-158

Donev, A., I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin (2004). "Improving the Density of Jammed Disordered Packings Using Ellipsoids". *Science* 303.5660, pp. 990-993. DOI: 10.1126/science.1093010. eprint: https://www.science.org/doi/pdf/10.1126/science.1093010. URL: https://www.science.org/doi/abs/10.1126/science.1093010
Bouvard, D. and R. M. McMeeking (1996). Deformation of interparticle necks by diffusion controlled creep. *Journal of the American Ceramic Society* 79(3):666-672. [DOI].

Bustamante, D., A. X. Jerves, and S. A. Pazmiño (2020). A generalized three-dimensional discrete element method with electrostatic induced cohesion. *Granular Matter* 22(4):1-13. [DOI].

Christoffersen, J., M. M. Mehrabadi, and S. Nemat-Nasser (1981). A micromechanical description of granular material behavior. *Journal of Applied Mechanics* 48(2):339-344. [DOI].

Coble, R. L. (1958). Initial sintering of alumina and hematite. *Journal of the American Ceramic Society* 41(2):55-62. [DOI].

Coblenz, W. S., J. M. Dynys, R. M. Cannon, and R. L. Coble (1980). Initial stage solid state sintering models. A critical analysis and assessment. *Sintering Processes* 13(1980):141-157. [HAL].

Davis, A. D., B. A. West, N. J. Frisch, D. T. O'Connor, and M. D. Parno (2021). ParticLS: Object-oriented software for discrete element methods and peridynamics. *Computational Particle Mechanics* 9(1):1-13. [DOI], [ARXIV].

Dervieux, A. and F. Thomasset (1980). A finite element method for the simulation of a Rayleigh-Taylor instability. *Approximation Methods for Navier-Stokes Problems*. Springer, pp 145-158. [DOI].

Donev, A., I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin (2004). Improving the density of jammed disordered packings using ellipsoids. *Science* 303(5660):990-993. [DOI].

Duriez, J. and S. Bonelli (2021). Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM. *Computers and Geotechnics* 134:104033. [DOI], [OA].

Džiugys, A. and B. Peters (2001). An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. *Granular Matter* 3(4):231-266. [DOI].

Eberly, D. (2020). *Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid*. [HAL].

Feng, Y. (2021). An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles: Basic framework and general contact model. *Computer Methods in Applied Mechanics and Engineering* 373:113454. [DOI], [OA].

Feng, Y., K. Han, and D. Owen (2012). Energy-conserving contact interaction models for arbitrarily shaped discrete elements. *Computer Methods in Applied Mechanics and Engineering* 205-208(1):169-177. [DOI].

Goldman, R. (2005). Curvature formulas for implicit curves and surfaces. *Computer Aided Geometric Design* 22(7):632-658. [DOI].

Grupp, R., M. Nöthe, B. Kieback, and J. Banhart (2011). Cooperative material transport during the early stage of sintering. *Nature Communications* 2(1):298. [DOI], [OA].

Harmon, J. M., D. Arthur, and J. E. Andrade (2020). Level set splitting in DEM for modeling breakage mechanics. *Computer Methods in Applied Mechanics and Engineering* 365:112961. [DOI].

Harmon, J. M., K. Karapiperis, L. Li, S. Moreland, and J. E. Andrade (2021). Particle bonding within the level set discrete element method. *Computer Methods in Applied Mechanics and Engineering* 373:113486. [DOI].

Hart, R., P. Cundall, and J. Lemos (1988). Formulation of a three-dimensional distinct element model-Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts* 25(3):117-125. [DOI].

Haven, D. L. van der, I. S. Fragkopoulos, and J. A. Elliott (2023). A physically consistent Discrete Element Method for arbitrary shapes using Volume-interacting Level Sets. *Computer Methods in Applied Mechanics and Engineering* 414:116165. [DOI], [OA].

Henrich, B., A. Wonisch, T. Kraft, M. Moseler, and H. Riedel (2007). Simulations of the influence of rearrangement during sintering. *Acta Materialia* 55(2):753-762. [DOI].

Houlsby, G. (2009). Potential particles: a method for modelling non-circular particles in DEM. *Computers and Geotechnics* 36(6):953-959. [DOI].

Hugonnet, B., J.-M. Missiaen, C. Martin, and C. Rado (2020). Effect of contact alignment on shrinkage anisotropy during sintering: Stereological model, discrete element model and experiments on NdFeB compacts. *Materials & Design* 191:108575. [DOI], [OA].

Karapiperis, K., J. Harmon, E. Andò, G. Viggiani, and J. E. Andrade (2020). Investigating the
Duriez, J. and S. Bonelli (2021). "Precision and computational costs of Level Set-Discrete Element Method (LS-DEM) with respect to DEM". *Computers and Geotechnics* 134.January, p. 104033.
DOI: 10.1016/j.compgeo.2021.104033. URL: https://doi.org/10.1016/j.compgeo.2021.104033

Džiugys, A. and B. Peters (2001). "An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers". *Granular Matter* 3.4, pp. 231–265. DOI: 10.1007/plo0010918

Feng, Y. T. (2021). "An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles: Basic framework and general contact model". *Computer Methods in Applied Mechanics and Engineering* 373, p. 113454. DOI: 10.1016/j.cma.2020.113454. URL: https://doi.org/10.1016/j.cma.2020.113454

Feng, Y. T., K. Han, and D. R. Owen (2012). "Energy-conserving contact interaction models for arbitrarily shaped discrete elements". *Computer Methods in Applied Mechanics and Engineering* 205-208.1, pp. 169–177. DOI: 10.1016/j.cma.2011.02.010. URL: http://dx.doi.org/10.1016/j.cma.2011.02.010

Goldman, R. (2005). "Curvature formulas for implicit curves and surfaces". *Computer Aided Geometric Design* 22.7 SPEC. ISS. Pp. 632–658. DOI: 10.1016/j.cagd.2005.06.005

Grupp, R., M. Nothe, B. Kieback, and J. Banhart (2011). "Cooperative material transport during the early stage of sintering". *Nat. Commun.* 2, p. 298. DOI: 10.1038/ncomms1038

Harmon, J. M., K. Karapiperis, L. Li, and S. Moreland (2021). "Particle bonding within the level set discrete element method". *Computer Methods in Applied Mechanics and Engineering* 373, p. 113486. DOI: 10.1016/j.cma.2020.113486. URL: https://doi.org/10.1016/j.cma.2020.113486

Harmon, J. M., D. Arthur, and J. E. Andrade (2020). "Level set splitting in DEM for modeling breakage mechanics". *Computer Methods in Applied Mechanics and Engineering* 365, p. 112961. DOI: 10.1016/j.cma.2020.112961. URL: https://doi.org/10.1016/j.cma.2020.112961

Hart, R., P. A. Cundall, and J. Lemos (1988). "Formulation of a three-dimensional distinct element model-Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks". *International Journal of Rock Mechanics and Mining Sciences and* 25.3, pp. 117–125. DOI: 10.1016/0148-9062(88)92294-2

Henrich, B., A. Wonisch, T. Kraft, M. Moseler, and H. Riedel (2007). "Simulations of the influence of rearrangement during sintering". *Acta Materialia* 55.2, pp. 753–762. DOI: 10.1016/j.actamat.2006.09.005

Houlsby, G. T. (2009). "Potential particles: a method for modelling non-circular particles in DEM". *Computers and Geotechnics* 36.6, pp. 953–959. DOI: 10.1016/j.compgeo.2009.03.001. URL: http://dx.doi.org/10.1016/j.compgeo.2009.03.001

Hugonnet, B., J. M. Missiaen, C. L. Martin, and C. Rado (2020). "Effect of contact alignment on shrinkage anisotropy during sintering: Stereological model, discrete element model and experiments on NdFeB compacts." *Materials and Design* 191, p. 108575. DOI: 10.1016/j.matdes.2020.108575. URL: https://doi.org/10.1016/j.matdes.2020.108575

Karapiperis, K., J. Harmon, E. Andò, G. Viggiani, and J. E. Andrade (2020). "Investigating the incremental behavior of granular materials with the level-set discrete element method". *Journal of the Mechanics and Physics of Solids* 144, p. 104103. DOI: 10.1016/j.jmps.2020.104103. URL: https://doi.org/10.1016/j.jmps.2020.104103

Kawamoto, R., E. Andò, G. Viggiani, and J. E. Andrade (2016). "Level set discrete element method for three-dimensional computations with triaxial case study". *Journal of the Mechanics and Physics of Solids* 91, pp. 1–13. DOI: 10.1016/j.jmps.2016.02.021. URL: http://dx.doi.org/10.1016/j.jmps.2016.02.021

Kawamoto, R., E. Andò, G. Viggiani, and J. E. Andrade (2018). "All you need is shape: Predicting shear banding in sand with LS-DEM". *Journal of the Mechanics and Physics of Solids* 111, pp. 375–392. DOI: 10.1016/j.jmps.2017.10.003

Lai, Z., Q. Chen, and L. Huang (2020). "Fourier series-based discrete element method for computational mechanics of irregular-shaped particles". *Computer Methods in Applied Mechanics and Engineering* 362, p. 112873. DOI: 10.1016/j.cma.2020.112873. URL: https://doi.org/10.1016/j.cma.2020.112873
incremental behavior of granular materials with the level-set discrete element method. *Journal of the Mechanics and Physics of Solids* 144:104103. [DOI].

Kawamoto, R., E. Andò, G. Viggiani, and J. E. Andrade (2016). Level set discrete element method for three-dimensional computations with triaxial case study. *Journal of the Mechanics and Physics of Solids* 91:1–13. [DOI], [HAL].

Kawamoto, R., E. Andò, G. Viggiani, and J. E. Andrade (2018). All you need is shape: Predicting shear banding in sand with LS-DEM. *Journal of the Mechanics and Physics of Solids* 111:375–392. [DOI], [OA].

Lai, Z., Q. Chen, and L. Huang (2020). Fourier series-based discrete element method for computational mechanics of irregular-shaped particles. *Computer Methods in Applied Mechanics and Engineering* 362:112873. [DOI].

Li, L., E. Marteau, and J. E. Andrade (2019). Capturing the inter-particle force distribution in granular material using LS-DEM. *Granular Matter* 21(3):1–16. [DOI].

Lu, G., J. Third, and C. Müller (2015). Discrete element models for non-spherical particle systems: From theoretical developments to applications. *Chemical Engineering Science* 127:425–465. [DOI].

Martin, C. L. and R. K. Bordia (2008). Influence of adhesion and friction on the geometry of packings of spherical particles. *Physical Review E* 77(3):1–8. [DOI].

Martin, C. L. (2024). *dp3D: DEM for materials science*.

Martin, C. and R. Bordia (2009). The effect of a substrate on the sintering of constrained films. *Acta Materialia* 57(2):549–558. [DOI].

Martin, C., L. Schneider, L. Olmos, and D. Bouvard (2006). Discrete element modeling of metallic powder sintering. *Scripta Materialia* 55(5):425–428. [DOI], [HAL].

Martin, S., M. Guessasma, J. Léchelle, J. Fortin, K. Saleh, and F. Adenot (2014). Simulation of sintering using a Non Smooth Discrete Element Method. Application to the study of rearrangement. *Computational Materials Science* 84:31–39. [DOI], [HAL].

Martin, S., S. Navarro, H. Palancher, A. Bonnin, J. Léchelle, M. Guessasma, J. Fortin, and K. Saleh (2016). Validation of DEM modeling of sintering using an in situ X-ray microtomography analysis of the sintering of NaCl powder. *Computational Particle Mechanics* 3(4):525–532. [DOI], [HAL].

Nosewicz, S., J. Rojek, K. Pietrzak, and M. Chmielewski (2013). Viscoelastic discrete element model of powder sintering. *Powder Technology* 246:157–168. [DOI].

Nosewicz, S., J. Rojek, and M. Chmielewski (2020). Discrete element framework for determination of sintering and postsintering residual stresses of particle reinforced composites. *Materials* 13(18):4015. [DOI], [OA].

Ogarko, V. and S. Luding (2012). A fast multilevel algorithm for contact detection of arbitrarily polydisperse objects. *Computer Physics Communications* 183(4):931–936. [DOI].

Osher, S., R. Fedkiw, and K. Piechor (2004). Level set methods and dynamic implicit surfaces. *Applied Mechanics Reviews* 57(3):B15–B15. [DOI], [OA].

Osher, S. and J. A. Sethian (1988). Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. *Journal of Computational Physics* 79(1):12–49. [DOI], [OA].

Paredes-Goyes, B., D. Jauffres, J.-M. Missiaen, and C. L. Martin (2021). Grain growth in sintering: A discrete element model on large packings. *Acta Materialia* 218:117182. [DOI], [HAL].

Paredes-Goyes, B., A. M. Venkatesh, D. Jauffres, and C. L. Martin (2022). Two-step sintering of alumina nano-powders: A discrete element study. *Journal of the European Ceramic Society* 43(2):501–509. [DOI], [OA].

Pazmiño, S. A., A. X. Jerves, J. Dijkstra, D. A. Medina, and H. P. Jostad (2022). A generalized 3DLS-DEM scheme for grain breakage. *Computer Methods in Applied Mechanics and Engineering* 399:115383. [DOI].

Petzow, G. and H. E. Exner (1976). *Particle rearrangement in solid state sintering*. [DOI].

Podlozhnyuk, A., S. Pirker, and C. Kloss (2017). Efficient implementation of superquadric particles in Discrete Element Method within an open-source framework. *Computational Particle Mechanics* 4(1):101–118. [DOI], [OA].

Pope, S. (2008). *Algorithms for Ellipsoids*. Tech. rep. FDA-08-01. Cornell University. [HAL].
Li, L., E. Marteau, and J. E. Andrade (2019). "Capturing the inter-particle force distribution in granular material using LS-DEM". *Granular Matter* 21.3, pp. 1-16. DOI: 10.1007/s10035-019-0893-7

Lu, G., J. R. Third, and C. R. Müller (2015). "Discrete element models for non-spherical particle systems: From theoretical developments to applications". *Chemical Engineering Science* 127, pp. 425-465. DOI: 10.1016/j.ces.2014.11.050

Martin, C. L. and R. K. Bordia (2008). "Influence of adhesion and friction on the geometry of packings of spherical particles". *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics* 77.3, pp. 1-8. DOI: 10.1103/PhysRevE.77.031307

Martin, C. L. and R. K. Bordia (2009). "The effect of a substrate on the sintering of constrained films". *Acta Materialia* 57.2, pp. 549-558. DOI: 10.1016/j.actamat.2008.09.041

Martin, C. L., L. C. Schneider, L. Olmos, and D. Bouvard (2006). "Discrete element modeling of metallic powder sintering". *Scripta Materialia* 55.5, pp. 425-428. DOI: 10.1016/j.scriptamat.2006.05.017

Martin, C. L. (2024). *dp3D: DEM for materials science*. URL: https://simap.grenoble-inp.fr/fr/equipe/s/animations-discrete-element-method

Martin, S., M. Guessasma, J. Léchelle, J. Fortin, K. Saleh, and F. Adenot (2014). "Simulation of sintering using a Non Smooth Discrete Element Method. Application to the study of rearrangement". *Computational Materials Science* 84, pp. 31-39. DOI: 10.1016/j.commatsci.2013.11.050. URL: http://dx.doi.org/10.1016/j.commatsci.2013.11.050

Martin, S., S. Navarro, H. Palancher, A. Bonnin, J. Léchelle, M. Guessasma, J. Fortin, and K. Saleh (2016). "Validation of DEM modeling of sintering using an in situ X-ray microtomography analysis of the sintering of NaCl powder". *Computational Particle Mechanics* 3.4, pp. 525-532. DOI: 10.1007/s40571-015-0062-7

Nosewicz, S., J. Rojek, and M. Chmielewski (2020). "Discrete Element Framework for Determination of Sintering and Postsintering Residual Stresses of". *Materials* 13.4015, pp. 1-21. DOI: 10.3390/ma13184015

Nosewicz, S., J. Rojek, K. Pietrzak, and M. Chmielewski (2013). "Viscoelastic discrete element model of powder sintering". *Powder Technology* 246, pp. 157-168. DOI: 10.1016/j.powtec.2013.05.020. URL: http://dx.doi.org/10.1016/j.powtec.2013.05.020

Ogarko, V. and S. Luding (2012). "A fast multilevel algorithm for contact detection of arbitrarily polydisperse objects". *Computer Physics Communications* 183.4, pp. 931-936. DOI: 10.1016/j.cp.c.2011.12.019. URL: http://dx.doi.org/10.1016/j.cpc.2011.12.019

Osher, S. and J. A. Sethian (1988). "Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations". *Journal of Computational Physics* 79.1, pp. 12-49. DOI: 10.1016/0021-9991(88)90002-2

Paredes-goyes, B., D. Jauffres, J.-m. Missiaen, and C. L. Martin (2021). "Grain growth in sintering: A discrete element model on large packings". *Acta Materialia* 218, p. 117182. DOI: 10.1016/j.actamat.2021.117182. URL: https://doi.org/10.1016/j.actamat.2021.117182

Paredes-Goyes, B., A. M. Venkatesh, D. Jauffres, and C. L. Martin (2022). "Two-step sintering of alumina nano-powders: A discrete element study". *J. Eur. Ceram. Soc.* 43, pp. 501-509. DOI: https://doi.org/10.1016/j.jeurceramsoc.2022.10.001

Pazmiño, S. A., A. X. Jerves, J. Dijkstra, D. A. Medina, and H. P. Jostad (2022). "A generalized 3DLS-DEM scheme for grain breakage". *Computer Methods in Applied Mechanics and Engineering* 399, p. 115383. DOI: 10.1016/j.cma.2022.115383. URL: https://doi.org/10.1016/j.cma.2022.115383

Petzow, G. and H. E. Exner (1976). *Particle Rearrangement in Solid State Sintering*. DOI: 10.1515/ijm r-1976-670906

Podlozhnyuk, A., S. Pirker, and C. Kloss (2017). "Efficient implementation of superquadric particles in Discrete Element Method within an open-source framework". *Computational Particle Mechanics* 4.1, pp. 101-118. DOI: 10.1007/s40571-016-0131-6

Pope, S. (2008). *Algorithms for Ellipsoids*. Tech. rep. February. Cornell University, p. 47. URL: http://tcg.mae.cornell.edu/pubs/Pope_FDA_08.pdf

Rasp, T., C. Jamin, O. Guillon, and T. Kraft (2017). "Cracking and shape deformation of cylindrical cavities during constrained sintering". *Journal of the European Ceramic Society* 37.8, pp. 2907-
Rasp, T., C. Jamin, O. Guillon, and T. Kraft (2017). Cracking and shape deformation of cylindrical cavities during constrained sintering. *Journal of the European Ceramic Society* 37(8):2907-2917. [DOI].

Salque, B. (2017). Caractérisation mécanique de la respiration des hydrures pour une conception optimisée des réservoirs de stockage de l'hydrogène par voie solide. French. PhD thesis. France: Université Grenoble Alpes. [OA].

Storåkers, B., N. Fleck, and R. McMeeking (1999). The viscoplastic compaction of composite powders. *Journal of the Mechanics and Physics of Solids* 47(4):785-815. [DOI].

Teixeira, M. H. P., V. Skorych, R. Janssen, S. Y. G. González, A. De Noni Jr, J. B. Rodrigues Neto, D. Hotza, and M. Dosta (2021). High heating rate sintering and microstructural evolution assessment using the discrete element method. *Open Ceramics* 8:100182. [DOI], [OA].

Thornton, C. and Z. Ning (1998). A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres. *Powder Technology* 99(2):154-162. [DOI].

Torquato, S., T. M. Truskett, and P. G. Debenedetti (2000). Is random close packing of spheres well defined? *Physical Review Letters* 84(10):2064-2067. [DOI], [ARXIV].

Venkatesh, A. M. (2023). In-situ X-ray nano-tomography analysis of ceramic powder sintering. PhD thesis. France: Université Grenoble Alpes. [OA].

Vlahinić, I., E. Andò, G. Viggiani, and J. E. Andrade (2014). Towards a more accurate characterization of granular media: Extracting quantitative descriptors from tomographic images. *Granular Matter* 16(1):9-21. [DOI].

Wachs, A., L. Girolami, G. Vinay, and G. Ferrer (2012). Grains3D, a flexible DEM approach for particles of arbitrary convex shape - Part I: Numerical model and validations. *Powder Technology* 224:374-389. [DOI].

Wang, C. and S. Chen (2013). Application of the complex network method in solid-state sintering. *Computational Materials Science* 69:14-21. [DOI], [OA].

Yan, Z., S. Hara, and N. Shikazono (2017). Effect of powder morphology on the microstructural characteristics of La<sub>0.6</sub>Sr<sub>0.4</sub>Co<sub>0.2</sub>Fe<sub>0.8</sub>O<sub>3</sub> cathode: A Kinetic Monte Carlo investigation. *International Journal of Hydrogen Energy* 42(17):12601-12614. [DOI], [OA].

Zhan, L., C. Peng, B. Zhang, and W. Wu (2020). A surface mesh represented discrete element method (SMR-DEM) for particles of arbitrary shape. *Powder Technology* 377:760-779. [DOI], [OA].

Zheng, Q., Z. Zhou, and A. Yu (2013). Contact forces between viscoelastic ellipsoidal particles. *Powder Technology* 248:25-33. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** Brayan Paredes-Goyes contributed to conceive the study, to develop the methodology, to implement the code, to perform numerical simulations and to draft the manuscript. David Jauffres contributed to conceive the study, to develop the methodology, to co-supervise the project and to critically review the manuscript. Christophe L. Martin contributed to conceive the study, to perform numerical simulations, helped with implementation and numerical issues, to coordinate the project and to critically review the manuscript.

**Supplementary Material** A dataset containing three folders can be found at the permalink 10.5281/zenodo.15168970. The input folder includes templates of dp3D input files used in the manuscript simulations. The data folder includes all the input data of the figures. The workflows contains scripts to run the simulations and reproduce the figures.

**Funding** This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement MATHEGRAM No 813202.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.
2917. DOI: 10.1016/j.jeurceramsoc.2017.03.013. URL: http://dx.doi.org/10.1016/j.jeurceramsoc.2017.03.013

Salque, B. (June 2017). "Caractérisation mécanique de la respiration des hydrures pour uneconception optimisée des réservoirs de stockage de l'hydrogène par voie solide". Theses. Université Grenoble Alpes. URL: https://tel.archives-ouvertes.fr/tel-01688042

Storåkers, B., N. A. Fleck, and R. M. McMeeking (1999). "The viscoplastic compaction of composite powders". *Journal of the Mechanics and Physics of Solids* 47.4, pp. 785-815. DOI: 10.1016/S0022-5096(98)00076-3

Teixeira, M. H. P., V. Skorych, R. Janssen, S. Y. G. González, A. De Noni, J. B. Rodrigues Neto, D. Hotza, and M. Dosta (2021). "High heating rate sintering and microstructural evolution assessment using the discrete element method". *Open Ceramics* 8.September. DOI: 10.1016/j.oceram.2021.100182

Thornton, C. and Z. Ning (1998). "A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres". *Powder Technology* 99.2, pp. 154-162. DOI: 10.1016/S0032-5910(98)00099-0

Torquato, S., T. M. Truskett, and P. G. Debenedetti (2000). "Is random close packing of spheres well defined?" *Physical Review Letters* 84.10, pp. 2064-2067. DOI: 10.1103/PhysRevLett.84.2064

van der Haven, D. L., I. S. Fragkopoulos, and J. A. Elliott (2023). "A physically consistent Discrete Element Method for arbitrary shapes using Volume-interacting Level Sets". *Computer Methods in Applied Mechanics and Engineering* 414, p. 116165. DOI: https://doi.org/10.1016/j.cma.2023.116165. URL: https://www.sciencedirect.com/science/article/pii/S004578252300289X

Venkatesh, A. M. (2023). "Analyse du frittage de poudres céramiques par nano-tomographie aux rayons X in-situ". 2023GRALI055. PhD thesis. URL: http://www.theses.fr/2023GRALI055/document

Vlahinić, I., E. Andò, G. Viggiani, and J. E. Andrade (2014). "Towards a more accurate characterization of granular media: Extracting quantitative descriptors from tomographic images". *Granular Matter* 16.1, pp. 9-21. DOI: 10.1007/s10035-013-0460-6

Wachs, A., L. Girolami, G. Vinay, and G. Ferrer (2012). "Grains3D, a flexible DEM approach for particles of arbitrary convex shape - Part I: Numerical model and validations". *Powder Technology* 224, pp. 374-389. DOI: 10.1016/j.powtec.2012.03.023. URL: http://dx.doi.org/10.1016/j.powtec.2012.03.023

Wang, C. and S. Chen (2013). "Application of the complex network method in solid-state sintering". *Computational Materials Science* 69, pp. 14-21. DOI: 10.1016/j.commatsci.2012.11.020. URL: http://dx.doi.org/10.1016/j.commatsci.2012.11.020

Yan, Z., S. Hara, and N. Shikazono (2017). "Effect of powder morphology on the microstructural characteristics of Lao.6Sro.4Coo.2Feo.8O3 cathode: A Kinetic Monte Carlo investigation". *International Journal of Hydrogen Energy* 42.17, pp. 12601-12614. DOI: 10.1016/j.ijhydene.2017.03.136. URL: http://dx.doi.org/10.1016/j.ijhydene.2017.03.136

Zhan, L., C. Peng, B. Zhang, and W. Wu (2020). "A surface mesh represented discrete element method (SMR-DEM) for particles of arbitrary shape". *Powder Technology* 377, pp. 760-779. DOI: 10.1016/j.powtec.2020.09.046. URL: https://doi.org/10.1016/j.powtec.2020.09.046

Zheng, Q. J., Z. Y. Zhou, and A. B. Yu (2013). "Contact forces between viscoelastic ellipsoidal particles". *Powder Technology* 248, pp. 25-33. DOI: 10.1016/j.powtec.2013.03.020