1. Introduction

Solid-state sintering is one of the most widely used material processing technologies in modern manufacturing ^[1]. Although sintering has been utilized for decades and novel sintering techniques are continuously being developed, it is still a nontrivial task to establish rigorous analytical or numerical models for industrial scale problems, due to the complex multi-physical processes involved in modern sintering technologies ^[2]. During sintering, material densification can be achieved either with or without externally applied pressure. In the latter case (i.e., free sintering), densification is driven by surface energy reduction and mainly achieved via diffusion of atoms or vacancies. Therefore, a fully densified sintered material is generally difficult to achieve while the final porosity is a key parameter to determine the properties and performance of sintered materials ^[3,4,5]. Due to its great importance, a significant amount of research work has been devoted to understanding and modeling microstructure evolution in free sintering ^[6].

The focus of this work is on sintering with externally applied pressure, which is also referred to as hot pressing sintering. During a hot pressing process, the overall densification can be achieved via many processes, including particle rearrangement ^[7], particle deformation and fracture, plastic flow ^[8], diffusion ^[9,10], etc. Under different circumstances, a single or a combination of these mechanisms may dominate a sintering process. Analytical models on pressure assisted diffusion and creep have been established, which lead to accurate knowledge of the effective stress on the grain boundaries and power-law creep behavior ^[9,11]. It has been shown that diffusion and creep last throughout the entire sintering process, while particle rearrangement and particle fracture occur in very early sintering stages and complete within seconds to minutes in an typical hot pressing sintering process ^[9]. However, the influence on the overall densification of this short initial period is significant. During the initial sintering stages, particles are squeezed by external pressure into free space and plastically deform their shape to accommodate the external forces, while particle fracture may also occur. Subsequently, interfacial diffusion, vacancy annihilation and other mechanisms gradually take over as the major driving forces, leading to neck formation and growth as well as further densification ^[12]. Therefore, in present work, we consider that a typical hot pressing sintering process consists of two distinct stages: initial powder compaction and subsequent diffusion controlled morphology evolution and densification.

To date, several compaction models have been proposed for particle systems at different length scales and under different processing conditions ^{[13,14,15,16]}. These models can be generally classified into two categories: discrete-particle-based models (e.g., discrete element model (DEM)) and continuum models (e.g., finite element method) ^[6]. For example, Martin et al. ^[17] used DEM simulations to investigate the effects of contact law, the relative density and the type of stress on plastic deformation and mutual rearrangement of particles during compaction. Chtourou et al. ^[18] employed finite element method and a "cap" constitutive model to simulate compaction processes, suggesting the choice of stress and strain measures. In general, the compaction process driven by the external pressure may result in grain boundary sliding, grain rotation, particle deformation and even particle rupture ^[19]. To the best of our knowledge, no existing hot pressing sintering models have incorporated the aforementioned compaction mechanisms.

In this work, we develop a novel hot pressing sintering model which incorporates powder compaction mechanisms including simultaneously rigid grain motions (sliding and rotation) and plastic grain deformations in the initial stages as well as interfacial diffusion and coarsening in the later stages. In particular, we combine the first-order compaction model proposed by Janssen and Walker (the J-W model) ^[20] and a recently developed kinetic Monte Carlo sintering method ^[21] to simulate the sintering behaviors of pixel-based discrete particle systems under external pressure. Specifically, uniaxial compression is considered, in which uniform pressures are applied at two opposite boundaries of the powder compact. The stress within the sintering system decays exponentially as one moves away from the boundaries and reaches a minimal value at the center of the system. Rigid-grain motions including translations and rotations are directed by the local stress state. A novel "geometric force" is introduced to realistically model geometrically necessary plastic particle deformation to accommodate the fast densification due to external pressure. The subsequent evolution of overall morphology and inter-particle connection is dominated by interfacial energy minimization and coarsening is considered in the later sintering stage. The utility of our model is illustrated by sintering 2D compacts of poly-dispersed circular particles and equal-sized elliptical particles, both from an initial loosely packed configuration to a dense sintered microstructure. Interestingly, we find that significant coarsening occurs in the poly-dispersed circle system, while only slight grain growth is observed in the equal-sized ellipse system. Although illustrated with 2D examples, our model is readily applicable to 3D systems.

The rest of the paper is organized as follows: In Sec. 2, we describe powder compaction model. In Sec. 3, we provide algorithmic details of our model for simulating sintering with hot pressing. In Sec. 4, we apply the model to two distinct 2D systems. In Sec. 5, we make concluding remarks.

---

2. Powder Compaction under Hot Pressing

---

2.1. Mechanisms for Compaction under Hot Pressing

The Janssen-Walker (J-W) model was first developed for 3D powder compact in a cylindrical container under uniaxial pressure ^[22,23,24]. This model can be readily applied to 2D systems in rectangular containers which can be regarded as lateral cut of a cylindrical container. In particular, the J-W model assumes that stress transmission ratio (the ratio of transmitted stress over the applied punch pressure) is an exponential function of compact height, material properties, container shape, Coulombic friction coefficient between powder and wall and applied punch pressure ^[25,26,27]. This ratio decays exponentially as the thickness of the compact increases. For uniaxial pressure applied at opposite boundaries of the container, the largest stress occurs at the boundaries, which reaches a minimal value at the center (see Figure 1).

When external pressure is applied on both ends of the powder compact, the particles will be forced to move towards the central region forming a densified green body. If the particles are hindered by their neighbors, they may slide, rotate, deform and even break into parts to get around the obstacles. Although complicated, the particle-particle interactions and responses can be generally categorized into two types: rigid body movement and plastic deformation. Rigid body movement is dominated by total mass, moment of inertia, the direction and magnitude of total contact force and torque with respect to mass of center, etc. Plastic deformation is controlled by material stiffness, elastic strain energy, and plasticity, etc. We note that grain rupture, which involves detailed calculation of stress distribution within each sintering particle and fracture mechanics analysis that is beyond the scope of our sintering model.

In general, the coupling of different physical mechanisms and complex nature of hot pressing sintering make it extremely difficult, if not impossible, to rigorously treat each individual mechanism involved both analytically and numerically. Therefore, we have devised a functional, called "cohesive energy", which quantitatively characterizes particle geometry and can be used to calibrate plastic deformation during particle-particle interaction ^[21]. For a proper definition of the cohesive energy, we need to accurately quantify particle geometry and understand how this geometric quantity varies due to local stress and interface evolution. Furthermore, the definition of cohesive energy should lead to a relatively simple way to calculate energy change before and after rigid body movement and plastic flow.

---

2.2. Geometric Moments and Cohesive Energy

Here we use the central geometric moment ^[28] to quantify the shape of a particle, which also affects the center-of-mass (COM) movement due to the total contact force exerted on a particle. Note that the geometric moments have been extensively used in geometry recognition ^[29,30] and it has been shown that a finite number of lower-order geometric moments (e.g., first and second order) can be sufficient to successfully reconstruct a relatively complex morphology ^[31], which indicates that these moments can be employed to quantify a geometrical shape of a complex object.

For a 2D image (or field) $f\left( x,~y \right)$, the geometric moment of ${{\left( p+q \right)}^{th}}$ order is defined as ^[32]:

${{m}_{pq}}=\iint_{\xi }{{{x}^{p}}}{{y}^{q}}f\left( x,y \right)dxdy$

(1)

where denotes integration domain. For a discrete image I(x, y), its ${{\left( i+j \right)}^{th}}$ moment can be calculated as:

${{m}_{ij}}=\underset{x}{\mathop \sum }\,\underset{y}{\mathop \sum }\,{{x}^{i}}{{y}^{j}}I\left( x,y \right)$

(2)

The geometric center is then defined as:

$\bar x = {{{m_{10}}} \over {{m_{00}}}},\bar y = {{{m_{01}}} \over {{m_{00}}}}$

(3)

The central moment of order of $\left( p+q \right)$ is defined as:

${{\rm{u}}_{{\rm{pq}}}} = \mathop \sum \limits_{\rm{x}} {\mkern 1mu} \mathop \sum \limits_{\rm{y}} {\mkern 1mu} {\left( {{\rm{x}} - {\rm{\bar x}}} \right)^{\rm{p}}}{\left( {{\rm{y}} - {\rm{\bar y}}} \right)^{\rm{q}}}{\rm{I}}\left( {{\rm{x}},{\rm{y}}} \right)$

(4)

The cohesive energy is defined as:

$\text{E}\left( \text{I}\left( \text{x},\text{y} \right) \right)={{\text{ }\!\!\alpha\!\!\text{ }}_{0}}{{\text{u}}_{00}}+{{\text{ }\!\!\alpha\!\!\text{ }}_{10}}{{\text{u}}_{10}}+{{\text{ }\!\!\alpha\!\!\text{ }}_{01}}{{\text{u}}_{01}}+{{\text{ }\!\!\alpha\!\!\text{ }}_{20}}{{\text{u}}_{20}}+{{\text{ }\!\!\alpha\!\!\text{ }}_{02}}{{\text{u}}_{02}}+{{\text{ }\!\!\alpha\!\!\text{ }}_{11}}{{\text{u}}_{11}}$

(5)

where ${{\alpha }_{0}}=10$, ${{\alpha }_{10}}={{\alpha }_{01}}=5$, ${{\alpha }_{20}}={{\alpha }_{02}}={{\alpha }_{11}}=2$. Note that these values are simply chosen based on computational convenience and subsequent calibrations.

For simplicity, in our model, we use a vector to store the chosen geometric moments to represent the geometry of individual particles, under the assumption that the associated cohesive energy is linearly related to weighted sum of the stored moments. This assumption is validated by the fact that it leads to realistic particle morphology evolutions during sintering. A more detailed description for this form of cohesive energy was provided in ref. ^[21]. The difference between geometric moments' values before and after particle deformation is then a natural measurement of the "response" of a geometric deformation, which is defined as the "geometric force" here, i.e.,

$F=\partial {{E}_{2}}-\partial {{E}_{1}}$

(6)

Detailed procedures for modeling force-directed rigid-boy movements and cohesive-energy-driven particle deformation via kinetic Monte Carlo method are described in subsequent section.

---

3. Algorithmic Details for Hot Pressing Sintering Model

As discussed in Sec. 1, hot pressing sintering process involves an initial fast densification stage ^[33] similar to a typical compaction process and a subsequent diffusion controlled morphology evolution stage due to interfacial energy minimization. Although the interface-energy-driven diffusion takes place during the entire sintering process, rigid-body motion and the resulting plastic deformation are the dominant mechanisms for microstructure evolution in the fast densification stage typically within a few seconds. Therefore, in our model, we consider three distinct stages, for which the major mechanisms for structural evolution are distinct. In the initial stages, only rigid body movement and particle deformation due to particle-particle and particle-wall collision will be explicitly modelled; in the intermediate stages, neck formation and growth due to vacancy annihilation at interface will be treated by interfacial energy minimization; finally, coarsening and grain growth will be considered in the later sintering stages. The sintering time is simply measured by the Monte Carlo steps and we consider the initial fast densification stage involving rigid particle motion virtually does not consume sintering time.

---

3.1. Modeling Powder Compaction

---

3.1.1. Rigid-body Translation during Hot Pressing

Our compaction model is established for powder compact under uniaxial pressure applied at both ends of the container ^[34]. The process of compaction is simulated here by selecting a particle and moving it along the direction of total contact force with a distance linearly related to the magnitude of the total contact force. The magnitude of the contact force at a specific contact point is simply the magnitude of the transmitted stress at that position, and its direction is along the propagation direction of the transmitted stress. The probability of a particle being chosen is calculated according to the magnitude of transmitted stress at the center of mass of that particle.

As discussed in Sec. 2.1, the stress transmission ratio as given by the J-W model ^[35,36] is:

$T\left( h \right)=\frac{{{\sigma }_{t}}}{{{\sigma }_{a}}}=\text{exp}\left( -C\frac{h}{D}~ \right)$

(5)

where ${{\sigma }_{t}}$ and ${{\sigma }_{a}}$ are respectively the transmitted and applied stresses, $C$ is a material dependent scaling parameter, $D$ is a characteristic length scale (e.g., the average particle diameter) and $h$ denotes the position in the compact. In the present model, we consider the probability of displacing a particle at position $h$ possesses the same mathematical form as stress transmission ratio, i.e.,

$\text{P}\left( \text{h} \right)=\exp \left( -\text{C}\frac{\text{h}}{\text{D}}\text{ }\!\!~\!\!\text{ } \right)=\text{exp}\left( -\text{ }\!\!\alpha\!\!\text{ h} \right)$

(6)

Due to the symmetry of our system, if we set $h$ to be 0 at one end and equal to $L$ at the other end, the probability $P$ can be re-written as:

$\text{P}\left( \text{h} \right)=\exp \left( -\text{C}\frac{\text{h}}{\text{D}}\text{ }\!\!~\!\!\text{ } \right)+\exp \left( -\text{C}\frac{\text{L}-\text{h}}{\text{D}}\text{ }\!\!~\!\!\text{ } \right)=\exp \left( -\text{ }\!\!\alpha\!\!\text{ h} \right)+\text{exp}\left( -\text{ }\!\!\alpha\!\!\text{ }\left( \text{L}-\text{h} \right) \right)$

(7)

where $C$ and $D$ are unified in the parameter $\alpha$, and $L$ is the linear length of the system along the pressing direction. Both $C$ and $\alpha$ are material-property dependent parameters. We note that according Eq. (7), the probability $P$ is equal to one at both ends where pressure is applied, and reaches its minimal value at the center, while its actual value depends on physical parameters of powder, friction coefficient of container wall. In the numerical examples presented in this work, the value of $\alpha$ is simply set to be $-1.0~\times ~{{10}^{-4}}$ to give relatively large value of $P$.

The distance $l$ (as measured via pixel length) that a chosen particle moves is determined by the magnitude of total contact forces (stresses) on the particle surface multiplied by a trivial scaling factor, i.e.,

$\text{l }\!\!~\!\!\text{ }=\text{ }\!\!~\!\!\text{ }\!\!\mu\!\!\text{ }\underset{\text{i}}{\mathop \sum }\,{{\text{ }\!\!\sigma\!\!\text{ }}_{{{\text{t}}_{\text{i}}}}}\text{ }\!\!~\!\!\text{ }=\text{ }\!\!~\!\!\text{ }\!\!\mu\!\!\text{ }\underset{\text{i}}{\mathop \sum }\,{{\text{ }\!\!\sigma\!\!\text{ }}_{\text{a}}}\text{exp}\left( -\text{ }\!\!\alpha\!\!\text{ }{{\text{h}}_{\text{i}}} \right)$

(8)

where $\mu$ is the length parameter, which depends on the sintering temperature and is chosen to be $0.05L$ in our simulations, and the summation is over the contact neighboring particles that contribute to densification. In general, this parameter can be calibrated from sintering experiments.

---

3.1.2. Shape Deformation due to Particle-Particle Interaction

As discussed in Sec. 2, in response to the contact forces, rigid body movement, including translation and rotation, take place frequently. However, if the rigid body movement is locally hindered, plastic deformation occurs if the magnitude local geometric force exceeds elastic limit. In our model, possible plastic flow ^[37] is realized by applying a simulated annealing scheme, i.e., randomly interchanging surface pixels to reduce the increment of cohesive energy as much as possible. In general, there are four major parameters in a typical simulated annealing process: total iteration number, error tolerance, cooling rate and initial temperature. These parameters are tested and optimized so that the selected parameter values lead to physically realistic sintering result.

Due to the fact that both rigid-body motion due to powder compaction and interfacial evolution take place simultaneously, formation and growth of necks between neighboring particles are controlled by multiple factors, including diffusion properties, local stress and temperature ^[38]. To simplify the scenario, in present model, the necks connecting neighboring particles can be artificially broken temporarily after the particle displacement directed by total contact force, which results in a slight increase of total surface energy. In the initial sintering stages, neck-breaking does not play an important role in microstructure evolution because neck area is relatively small and weak, which is easy to break. However, in intermediate and later stages, neck-breaking events may significantly increase the interfacial energy and affect subsequent microstructure evolution. Therefore, when implement the geometric-force-directed plastic flow, we not only need to minimize the geometric moment variations (i.e., cohesive energy) for the rigid-displaced particles but also the local change of interfacial energy arising from neck modification. The procedure for minimizing interfacial energy is discussed in Sec. 3.2, which allows us to statistically maintain the total interfacial energy of necks.

---

3.1.3. Rotation of Non-spherical Particles

As indicated in Sec. 1, during a powder compaction process, rotation takes place and must be explicitly considered for non-spherical grains. Similar to the procedure for force-directed grain translation described in Sec. 3.1.1, when the total torque on a specific particle due to contact forces is not zero, it will rotates to reduce the unbalanced part of the toque in order to achieve a new balance. Therefore, an appropriate definition of torque and physically reasonable rotating rule are necessary for modeling particle rotation during hot pressing sintering.

In classical mechanics, torque is generally defined as the cross product of displacement and contact force. In the present model, the total torque on a particle is computed according to the equation below:

$\text{N}=\overset{\text{n}}{\mathop \sum }\,{{\text{r}}_{\text{i}}}{{\text{f}}_{\text{i}}}{{\text{m}}_{\text{i}}}$

(9)

where $n$ is the total number of neighboring particles in contact, ${{r}_{i}}$ is the distance between the contact location of neighboring particle $i$ and the center of the mass of the interested particle, ${{m}_{i}}$ is mass of particle $i$, and ${{f}_{i}}$ is the transmitted stress at contact position, if the direction of ${{f}_{i}}$ is counterclockwise wither respect to the center of mass of the interested particle, $f$ is positive; otherwise, $f$ is negative.

Once the total torque is obtained, we need to determine the magnitude of rotation angle given the magnitude of unbalanced torque. Here we assume that the rotation angle linearly depends on the total contact torque, i.e.,

$\text{ }\!\!\theta\!\!\text{ }=\text{ }\!\!\gamma\!\!\text{ N}$

(10)

where $\text{ }\!\!\gamma\!\!\text{ }$ is an scaling constant and is chosen to be $\gamma =~\text{0}\text{.0001}$ to realistically produce the sintering dynamics. The rotation is implemented by transforming the position of each pixel of the rotating particle as follows:

$\text{{x}'}=\text{xcos }\!\!\theta\!\!\text{ }+\text{ }\!\!~\!\!\text{ ysin }\!\!\theta\!\!\text{ }$

(11)

$\text{{y}'}=-\text{xsin }\!\!\theta\!\!\text{ }+\text{ycos }\!\!\theta\!\!\text{ }$

(12)

and the central moment is then updated as in:

${{\rm{u}}_{{\rm{pq}}}} = \mathop \sum \limits_{\rm{x}} {\mkern 1mu} \mathop \sum \limits_{\rm{y}} {\mkern 1mu} {\left( {{\rm{x' - \bar x}}} \right)^{\rm{p}}}{\left( {{\rm{y' - \bar y}}} \right)^{\rm{q}}}{\rm{I}}\left( {{\rm{x'}},{\rm{y'}}} \right)$

(13)

note that we use the center of mass before rotation due to the fact that is quite small, and thus, the resulting center of mass displacement is negligible.

---

3.2. Interfacial-energy Driven Morphology Evolution

In next step, we focus on diffusion-driven microstructure evolution by minimizing interfacial energy which becomes a major driving force for densification ^[39,40]. To perform an overall diffusion-driven microstructure evolution to relax sintering stress ^[41], we first randomly pick a pixel at interface, and randomly exchange it with its neighbors within prescribed separation distance from that pixel. If the exchange leads to a lower interfacial energy, it is accepted; otherwise, it is accepted with a probability determined according to the Metropolis rule ^[42], with the initial temperature $\beta ~=~0.001$ and cooling rate $\gamma =0.99$. The system's total surface/interface energy is defined as:

${{E}_{S}}=\frac{1}{2}\underset{i}{\overset{N}{\mathop \sum }}\,\underset{j}{\overset{m}{\mathop \sum }}\,\frac{1}{{{d}_{ij}}}\left[ 1-\delta \left( {{q}_{i}},~{{q}_{j}} \right) \right]$

(14)

where $N$ is the total number of pixels in the system, $m$ is the number of neighboring pixels, ${{d}_{ij}}$ is the distance between two pixels and $\delta \left( {{q}_{i}},~~{{q}_{j}} \right)$ is Kronecker delta, which is one if ${{q}_{i}}$ and ${{q}_{j}}$ are of the same value, otherwise, it is zero. The effect of temperature is implicitly considered in the surface energy minimization, i.e., higher temperature results in faster surface diffusion and better minimized surface energy.

---

3.3. Coarsening and Grain Growth

During later stage of sintering, the final density and grain size distribution are all significantly influenced by coarsening process, which usually has negative effects on the mechanical properties of the sintered body. To incorporate the coarsening effect, we first track all of the interfacial pixels, then randomly pick one of them and change its value to be the same as one of its neighbors with equal probability. If this change decreases the total surface energy, it is accepted; if it increases the surface energy, it is accepted according to Metropolis rule with the initial temperature $\beta ~=~0.0001$ and cooling rate $\gamma =0.999$.

---

4. Results and Discussion

---

4.1. Sintering Compact of Monodisperse Elliptical Grains

We have postulated that geometric moment is an effective way to quantify the geometry of a particle and assist the implementation of rigid body movement, including translation and rotation of non-spherical grains. To validate this, we show that our model leads to realistic microstructure evolution and correctly capture the densification process for particles with non-spherical shapes. In particular, we simulate the hot pressing sintering process of the maximally random jammed (MRJ) packing of ~100 elliptical particles in 2D with aspect ratio 2:1 ^[43]. The simulation domain is discretized into pixels. The linear size of simulation domain is 200 pixels and the long semi-axis of the ellipse is 10 pixels.

Figure 2 shows the microstructure evolution of the ellipse packing during a hot pressing sintering. Clearly, it can be seen that the anisotropic nature of the particle shape is well preserved during the whole sintering process, resulting in elongated grains in the final sintered microstructure. However, the orientations of particles are evolved significantly from a random distribution to well-alignment along the direction perpendicular to the pressing direction due to rigid-body rotations. In addition, although coarsening is explicitly considered in the later stage, no significant coarsening or grain growth is observed. This is due to the fact that the grains are of the same size, and thus, no single grain is more energetically favorable than others. These results indicate the effectiveness of our methods for defining cohesive energy via geometric moments, handling rigid body movement in response to external pressure, as well as coupling interfacial diffusion and coarsening.

Figure 3 shows the density profile of the sintering body along y direction (i.e., the direction of applied pressure) at different stages. Figure 4 shows the simulation data and exponential fitting curve of the packing density in the central region (to avoid edge effect) during the sintering process. It can be easily seen that the major densification is achieved during the early sintering stages due to the force-directed rigid-body motions and the associated plastic deformation, consistent with experimental observations.

---

4.2. Sintering Compact of Polydisperse Circular Grains

We now apply our hot pressing model to simulate a powder compact composed of polydisperse circular grains in 2D. Specifically, a loose random packing of ~200 circular disks with a Gaussian distribution of radii with a mean of 10 pixels and a standard deviation of 5 pixels is generated using Monte Carlo simulation and subsequently employed as the initial configuration for the sintering simulation. The simulation domain is discretized in pixels and the linear size of simulation domain is 200 pixels. Figure 5 shows that particles located in the boundary regions of the compact undergo large movements and therefore exhibit faster densification due to the higher stress in those regions according to Janssen-Walker's model. Meanwhile, the particle shape is also evolved due to the plastic flow in response to the high stress level in such regions. Once the initial fast densification is complete, diffusion-driven interface evolution dominates and neck growth can be clearly observed. In the late sintering stage, significant system-wide coarsening and grain growth are observed.

Figure 6 demonstrates the density profile at different sintering stages, which is similar to the ellipse sintering described in Sec. 4.1. We note that the density is slightly smaller at the surface regions than in the central region, which arise from the edge effect. Figure 7 shows the density evolution as a function of sintering time (Monte Carlo stages) calculated in the central region of the system.

We note that different from hot pressing of equal-sized ellipses reported in Sec. 4.1, in which no coarsening was observed, in the sintering of polydisperse circular disks, significant coarsening of the grains are apparent. This is mainly due to the large particle size dispersity in the initial packing in the latter case. Finally, we find that in our model the final density is mainly determined by the following parameters: plastic deformation rate which depends on the material's yield strength, the total sintering time, and the magnitude of externally applied pressure. Calibration of these model parameters for real experiment requires a systematic investigation with combined experimental and computational efforts, which will be carried out in our future work. Here, we just select the parameter values to achieve realistic sintering dynamics for the model systems.

---

5. Conclusions

In summary, we have developed a new computational model for hot pressing sintering, which incorporates a variety of relevant densification mechanisms such as forced-directed rigid-boy motions, stress-induced plastic flow, interfacial diffusion as well as coarsening and grain growth. In particular, rigid-body movements of sintering particles including both translation and rotation are explicitly taken into account in the model and are directed by the local stress state of the sintering system. A cohesive energy and the resulting "geometric force" are defined and incorporated to stochastically model particle deformation due to plastic deformation induced by external pressure. Subsequent evolution of particle morphology and inter-particle connection are mainly driven by interfacial energy minimization. And coarsening is considered in the later stages. The utility of our model is illustrated by sintering 2D compacts of polydisperse circular grains and equal-sized elliptical grains.

For future work, we will apply the model to large-scale 3D powder compact to investigate the sintering behavior of industrial relevant sintering systems. For example, we will systematically investigate the effects of material properties, packing characteristics (e.g., particle shape, size distribution, initial packing density and contact distribution etc.), as well as sintering conditions (e.g., externally applied multi-axial pressure, temperature, etc.) on the microstructure evolution and mechanical properties of the final sintered material.

---

Acknowledgments

This work was supported by the Division of Materials Research at the National Science Foundation under award No. DMR-1305119. Y. J. also thanks Arizona State University for generous start-up funds.

---

Conflict of Interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.

---