A phase field model of pressure-assisted sintering

The incorporation of an efficient contact mechanics algorithm into a phase field sintering model is presented. Contact stresses on the surface of arbitrarily shaped interacting bodies are evaluated and built into the model as an elastic strain energy field. Energy relaxation through deformation is achieved by diffusive fluxes along stress gradients and rigid body motion of the deforming particles maintain contact between the particles. The proposed model is suitable for diffusion deformation mechanisms occurring at stresses below the yield strength of a defectfree material; this includes Nabarro-Herring creep, Coble creep and pressure-solution. The effect of applied pressure on the high pressure-high temperature (HPHT) liquid phase sintering of diamond particles was investigated. Changes in neck size, particle coordination and contact flattening were observed. Densification rates due to the externally applied loads were found to be in good agreement with a new theory which implicitly incorporates the effect of applied external pressure.


Introduction
Sintering is the thermal consolidation of discrete powder particles to form a coherent compact component and is an important fabrication process in the manufacturing of many engineering products. Certain materials that are used in demanding high temperature applications are inherently resistant to sintering by high temperature alone [1]. This includes materials such as certain ceramics (alumina, carbides), refractory metals (titanium, aluminium) and even some alloys (steels, Nimonic, TiAl) [2]. Other materials such as diamond and cubic boronitride (cBN) also require high pressures in order to maintain phase stability during sintering [3].
The classical driving force for bonding between the particles during sintering is mediated by surface energy reduction. During this process, a gradient in the chemical potential between the free particle surface and the inter-particle contacts results in atomic diffusion and bonding. During pressure-assisted sintering, the particles attain a complex stress state with significant elastic energy stored in the particle contacts. Differences in stress between the centre of the contact and free surface provide additional gradients in the chemical potential. Material transfer along the stress gradients result in deformation and particle centres moving closer to each other. Mechanical stresses thus serve to enhance the driving force and lead to an increase in the densification kinetics during sintering [1].
There are three main densification processes active during pressure assisted sintering, namely: cold compaction, hot compaction and diffusion. HPHT diamond sintering, shown schematically in Fig. 1, posseses all three. At relatively low temperatures, cold compaction dominates; stresses above the yield stress of the material result in densification through fracture and particle rearrangement. At higher temperatures, hot compaction becomes dominant; stresses above the yield strength result in dislocations and plastic deformation. At moderate to high temperatures, mechanical stresses below the yield strength still result in densification via creep deformation or diffusion [4]. This paper will focus on the latter of the three processes. For an elastic defect-free material, the diffusive process can be further subdivided into three distinct mechanisms: Nabarro-Herring creep [5], Coble creep [6] and pressure-solution [7]. The type of deformation mechanism experienced by the material is dependent on the sintering conditions and material properties. The dominant mechanism is determined by the fastest diffusional pathway. The pathways for Nabarro-Herring creep, Coble creep and pressure-solution are lattice, grain-boundary and surface diffusion respectively. In reality, the deformations operate simultaneously with the densification rate summed over the active mechanisms [4]. In the case of diamond sintering, a metal catalyst, typically cobalt or nickel, becomes liquid at the sintering temperature and provides a much faster diffusional pathway for material transfer. Many different computational models describing the microstructure evolution during sintering exist; notable papers and reviews of these models are available [8,9,1,10,11] stresses in these models is, however, sparsely reported in the literature; this is because tracking of the particle interfaces during deformation is not trivial, especially when other processes such as surface energy reduction due to inter-particle bonding are considered. Bjørk [12] incorporated stresses into a sintering model by coupling the vacancy annihilation phenomenon in a kinetic Monte Carlo algorithm with the finite element method. This multi-scale approach required that the microstructure be re-meshed at each time step as it evolved. Recent advances in rigid body motion [8] and model dimensionalisation [13] have made the phase field simulation of sintering an extremely attractive option. One of the biggest strengths of phase field models is that they allow the evolution of complex arbitrary interfaces to be resolved in time without having to explicitly track their exact position. Furthermore, the versatility of phase field models means that new physics can easily be incorporated. A review of the phase field method for microstructure evolution is available [14]. To this end, an elastic phase field model for pressure assisted sintering is presented.
The development of micro-elasticity theory and its effect on microstructure evolution in phase field models has previously been reported in the literature; however, these models typically concern themselves with examples such as the strain energy resulting from elastic mismatch at the interfaces [15] and phase transformation strain [16]. Other examples include the effect of elastic strain energy on precipitate growth [17] and a dissolution front [18]. The inclusion of elastic strain in phase field sintering models is, however, limited. The effect of strain energy on the grain growth kinetics for a pore-free fully dense polycrystalline material has been investigated [19].
The application of a general micro-elasticity theory to a phase field sintering model for a porous material is presented in [20] and [21]. In both cases presented, the particles were not initially in contact. Displacement boundary conditions were applied to the edges of the simulation domain and the void space surrounding the particles acted as a pressure transmitting medium. Elasticity throughout the domain was computed using a coupled fully implicit finite-element solver. While the incorporation of micro-elasticity in the aforementioned sintering model would allow the study of the evolution of sintering particles subjected to external loads and strong contact interaction between individual grains, the implementation of a contact algorithm to allow the study of pressure-assisted densification in these models is missing. Incorporation of a contact algorithm in such a system would require an advanced finite-element/discrete-element approach and would increase the computational time significantly, especially in the presence of a large number of particles. In light of this limitation and in the absence of a model that explicitly captures the interactions between grains due to contact and related elastic deformations, we propose the implementation of a modeling strategy that allows incorporation of contact mechanics within a phase-field sintering model to capture pressure-assisted densification. The presented model allows the calculation of stresses at the inter-particle contact regions to be calculated semi-analytically and thus scales to systems involving very many particles with minimal computational overhead.
In Section 2, the theoretical framework is presented. The model is then applied to the high pressure liquid phase sintering of diamond particles and details of the computational implementation are provided in Section 3. The effect of external pressure on a single particle pressing onto a non-interacting rigid wall, two sintering particles experiencing a uniaxial load and lastly, the high-pressure sintering of a multi-particle compact are investigated in Section 4. To end with, the paper is concluded by providing and discussing the simulation results and evaluating opportunities for further research within the field in Section 5.

Formulation
The model is constructed as a modification of the advection-diffusion sintering phase field model initially proposed by Wang [8] and later extended to include elastic energy by Biswas et al. [20]. The proposed model thus shares the same free energy functional and phase field equations as in [20], however, as will be shown later, differs significantly in that we proposed an original way to deal with contact stresses and to implement them within the phase field framework. Within this formalism, the solid particles are initiated by a scalar field c (r, t), called the phase field variable. The phase field variable represents the mass density, or in the case of diamond sintering, the fractional solute concentration of carbon in the liquid metal phase. The solute concentration is a function of position, r, and time, t. Within a solid particle, the solute concentration is considered fully saturated and c takes a value of 1. Conversely, in the pure liquid phase, c takes a value of 0. At the interface between the solid and liquid there exists a smooth diffuse region where c ∈ (0, 1). The particles are further described by another scalar field, known as the order parameter. Each particle in the system possesses its own order parameter field, η i (r, t), where i = 1, 2, …, p, is the index of the particle and p is the total number of unique particles in the simulation domain. The order parameter is thus used to describe the unique identity of each particle with its own crystallographic orientation. A schematic of two contacting particles and their phase field and order parameters is shown in Fig. 2.
According to the principles of gradient thermodynamics [22], the free energy, F, of the system within the domain, V, is constructed as follows: (1) where f o is the non-equilibrium bulk free energy, ∇ is the gradient operator, β c and β η are the gradient energy coefficients and f el is the elastic strain energy density resulting from inter-particle contact. The first term of Eq. (1) describes the energy associated with a homogenous mixture of the different phases and takes the form of the following multi-well Landau polynomial: A and B are material specific constants related to the height of the energy barriers for phase mixing. The bulk free energy term thus acts to separate out a mixture of phases. The second term in Eq. (1) describes the energy associated with the gradient of the free surfaces between the solid and liquid phase. The third term describes the energy associated with gradients in the order parameter and is related to the grain boundary energy. The conservation of phases is enforced by the continuity equation: where j is the total solute flux density through a particular point in space and is comprised of a diffusional flux, j dif , and an advective flux, j adv : In order to minimise the energy of the system, a variational derivative of the free energy with respect to the phase field variable gives rise to the generalised chemical potential as follows: The driving force for diffusion is the gradient of the chemical potential: where j dif is the diffusional flux density. Assuming a constant molar volume, V m , the mobility, M, can be estimated from the diffusion coefficient, D, according to the equation: where R is the universal gas constant and T the temperature. The diffusion coefficient can be decomposed into its surface, lattice and grain boundary diffusivities as demonstrated in [8]. This form of the diffusivity would enable the model to study pressure-solution, Nabarro-Herring and Coble creep deformation mechanisms by taking into account the respective diffusional pathways. By combining Eqs. (1), (3), (4) and (6), we obtain the modified form of the Cahn-Hilliard non-linear advection diffusion equation for the temporal evolution of the phase field variable: The temporal evolution of the non-conserved order parameter field for the ith grain is similarly derived by the variational derivative with respect to the order parameter and is given by the modified Allen-Cahn structural relaxation equation: where j i,adv is the advection flux for the ith grain and L is the Allen-Cahn mobility which is related to the grain boundary mobility M b by the relation: β η L = γ ss M b [13]. The diffuse interfaces are allowed to develop for a small finite amount of time, t init , and the contact patches between interacting grains Fig. 2. Schematic illustrating two contacting diffuse interface particles of radius r. Γ indicates the particle interior, p the centre of mass and Ω the diffuse particle boundary which has a width δ. A plot of the spatial variation of the phase field and order parameter variables along line H-H are displayed on top. The overlap of the diffuse interfaces Ω 12 along line K-K is plotted on the side with Ω 12,max the centre of the contact and a the diameter of the contact area. The blue and yellow arrows indicate the advection direction, N , and normal force vector, F n , acting on particle 1 respectively.
i and j are defined by the overlap of their order parameters fields according to the geometric operation: where q is the overlap threshold for identifying contact between adjacent particles. The operation Ω ij is typically used to identify grain boundaries, but its use has been extended to identify particle contact in the proposed model. The location of the maximum value of the operation, Ω ij,max is thus henceforth called the centre of the contact between particles i and j. By normalising the values of the order parameter fields in the overlap Ω ij with respect to the maximum value within the contact in question, the function Ω ij /Ω ij,max is used as the discrete approximation of the normalised stress distribution within the contact.
The peak stress in the centre of a spherical contact, σ 0 , is calculated as follows: where F n is the normal force acting on the particle surface and, a, is the contact area measured from the result of the function Ω ij . This form of the peak stress is applicable to spheroidal particles at the start of sintering as presented in this study. It's form would be different for nonspheroidal contacts and is the subject of further investigation. The elastic strain energy density in the contact is thus calculated as follows: where E is the elastic modulus of the material.
As material diffuses along elastic energy gradients away from the centre of the contact, rigid body motion is required to maintain contact between the two particle surfaces. The restorative force for any particle in contact with its partner is thus calculated as: N is direction vector between the centre of mass of the particle, q i , and the contact point centre, Ω ij,max and has a magnitude of unity. The second term is the magnitude of the advective force, and is calculated as the deviation of the phase field variable at the centre of the contact, c| ϕ ij,max , from the value obtained at t init which corresponds to the time when the vacancies are assumed to be saturated at the initial contact point and the elastic energy in the system is switched on, c| Ω t , ij,max init . As in Wang's model [8], the driving force behind the advective force is the vacancy annihilation mechanism which, in our case, is amplified by the over-saturation of vacancies resulting from elastic energy in the contact. The force thus ensures that surfaces stay in contact to maintain mechanical equilibrium from the applied pressure. Due to the fact that the majority of confined particles possess several contact partners, the contribution of torque acting on the particles is likely to be minimal and therefore rotational motion was omitted from the model. The advection velocity field acting on the order parameter of the ith particle is thus a sum of the pairwise forces resulting from all interacting particles: where m t is the translational mobility and V i is the volume of the ith particle. Strictly speaking, a force on the particle should result in its acceleration; however, since the particle motion due to the advective force is a highly dissipative process [8], the velocity of the particle is used instead. The advection flux density acting on the order parameters is then given by: The overall advection velocity field acting on the solid phase is thus a sum of the velocity fields acting on the respective order parameters:

Implementation
In order to demonstrate the model, the initial implementation and all simulations were performed in two dimensions. The phase field equations were dimensionalised and parametrised to a real material by relating the constants M, A, B, β c and β η to the diffusion coefficient, D, grain boundary energy, γ ss , Nickel-Diamond interface energy, γ sl , and grain boundary width, δ, according to the equations of Ahmed et al. [13]. In the present example, only two phases were considered at the chosen temperature of 1800 K, namely, the solid phase which is diamond and the liquid phase which is molten nickel. To make the model applicable to liquid phase rather than solid state sintering, the solid-vapour interfacial energy in Ahmed's work was replaced with the solid-liquid interfacial energy. The model parameters used for the Nickel-Diamond system are listed in Table 1.
Several assumptions were made to simplify the model for demonstration purposes: elasticity and diffusion were assumed to be isotropic and constant at the given temperature, only limited plastic deformation is expected and was therefore not accounted for and the grain-boundary energy for diamond was taken as the average value over many possible mismatch angles. It is also assumed that the inter-particle forces resulting from the surface tension of the liquid are negligible compared to those arising from the application of external pressure.
Pressure gradients in the liquid arising from rigid body motion and migration of the interfaces are likely to be small and negligible; convection in the liquid phase was thus not considered.
In order to strike a balance between computational efficiency, sufficient diffuse interface resolution and minimised diffuse interface effects at the chosen grain-size scale, careful consideration of the mesh spacing and diffuse interface width was necessary. General guidelines based on a quantitative analysis [23] are that the diffuse interface width be at least three mesh points, δ ≥ 3Δx with the radius of the smallest particle be at least the diffuse interface width R ≥ δ. Based on mesh convergence tests and the aforementioned guidelines, mesh 3.71 × 10 −6 J/m Calculated according to [13] Phenomenological constants Value Source not available, a typical ceramics grain boundary mobility selected from [28] was used instead. b Phenomenological constants were used for q, m t and δ. An intentionally large value for the diffuse interface width, δ, was chosen for computational expediency (it has been shown that this does not affect the kinetic or thermodynamic driving forces [31]).
spacing and interface width were fixed at Δx = 0.5 μm and δ = 5Δx respectively for all the presented simulations. Further details of each set of simulations are provided in the subsections that follow.

Single particle test
The first test case was a single particle pressing onto a non-interacting diffuse interface wall. A 30 μm particle was initiated as a discrete object with a sharp-interface boundary. The particle was then allowed to develop a diffuse interface for a short finite amount of time, t init , such that the position of the middle of the diffuse interface at c = 0.5, η 1 = 0.5 coincides with the position of its respective sharp interface. A similar procedure was done with the wall. The wall was then introduced to the particle as a constant η wall field only, with a null concentration field. A schematic of the setup is shown in Fig. 3. Thus, in order to make the wall non-interacting, there was no time integration of the Cahn-Hilliard or Allen-Cahn equations for η wall . By initiating the wall in this manner, it served only to provide an imaginary rigid surface upon which the particle was pressed. A downward body force of 0.4 N was exerted on the particle after t init = 2.94 s and the simulation was run for 260 s. This force is typical of what a particle might experience during the HPHT sintering of diamond.
Dirichlet boundary conditions for the particle were imposed at the non-interacting diffuse wall such that there was a zero energy gradient across the wall: Furthermore, the particle was allowed to advect to but not across the wall boundary: In order to impose the boundary conditions above, a second-order central finite-difference scheme and an explicit Euler scheme was used as the discrete approximation of all spatial and time derivatives respectively. The time-step for the single particle simulation was 0.00588 s.

Two particle test
The effect of contact pressure on the sintering kinetics of two 40 μm particles with 0 N, 0.4 N and 0.8 N uniaxial body forces was studied. The total simulation time was 16 s. Since a uniform mobility field was used for demonstration in this example, it was possible to solve the phase field equations numerically using the semi-implicit Fourier spectral method [32]. For simulations where a spatially variable mobility is desired, a finite difference algorithm is typically employed [33]. The timestep size used for the simulation was 0.0588 s with t init = 2.94 s.

Multi-particle tests
For the multi-particle tests, it was necessary to employ the discrete element method (DEM) to obtain the initial starting arrangement and normal forces between the particles. This was achieved by using the Granular package of the open-source DEM engine LAMMPS [34]. A bimodal particle size distribution of 15 μm and 25 μm particles with a number ratio of 1:1 was used. 50 particles were introduced randomly and the simulation domain was gradually reduced in size until the desired pressure on the domain walls was achieved. The DEM simulation was performed in two dimensions with periodic boundaries in the x and y directions. Particle interaction was done according to a frictionless Hertzian formulation with zero tangential stiffness. The particle positions, radii and pair-wise normal forces were then exported and introduced into the phase-field model. Due to the large number of order parameters, a bounding box algorithm [35] was used to improve computational efficiency. As in the two particle case, the Cahn-Hilliard equations were solved numerically using the semi-implicit Fourier spectral method with t init = 1 s for a total simulation time of 10.65 s with a time-step of 0.0588 s. A set of simulations at different external pressures, from 1 to 8 GPa, were performed and compared with the pressure-less case. Three sets at the respective pressures, each with different random starting configurations, were used to gather statistics. Average values and standard deviations between the three respective sets were obtained. The density of the sintering particles was calculated within a square region positioned inside the compacted specimen; to reduce edge effects, the sides of the square region were positioned such that their distance from the outermost particle edge of each side was a distance of two times the radius of the largest particle.

Results and discussion
The classical Hertzian contact theory for two spherical bodies was used to find contact pressure distributions from a geometrical consideration of their elastic displacements within the contact patch between the two mating surfaces [36]. Similarly, within the proposed diffuse interface formalism, the pressure distribution within the contact is approximated by the geometric overlap of the order parameters of the contacting bodies. Fig. 4 compares the contact function Ω ij /Ω ij,max for two overlapping circular diffuse interfaces to the normalised pressure distribution along the contact patch as calculated by the Hertzian stress equations [37]. In reality, many ceramic powders, diamond included, are not spheroidal. The particles possess an irregular and faceted morphology and the contact points are likely to be angular with some small degree of rounding at the tip due to plasticity. In a similar manner to the spheres, the contact function is compared to the analytical solution for the pressure distribution at the contact of a rounded-tip square indenter [38] in Fig. 5. There is good agreement between the diffuse contact function and the analytical formulations; deviations are attributed to the discretisation of the solid particles. This significant observation means that the proposed method is capable of calculating the contact stresses of arbitrarily shaped bodies fairly accurately in a very simple and efficient manner without the use of complicated analytical formulas. Furthermore, transient pressure distributions as the particle changes shape over time are implicitly handled by the model. This strengthens the advantages of using a diffuse interface approach by eliminating the need for tracking of particle surfaces.

Single particle test
The extent of pressure solution deformation on the contour of a single 30 μm particle experiencing a 0.4 N downward body force is shown in Fig. 6. The contour of the particle was taken at the location in the diffuse interface where c = 0.5. Bulging of the particle around the periphery of the contact was caused by diffusion of material from the centre of the contact, where the elastic strain energy was highest, to the Fig. 3. Schematic illustrating the diffuse interface of a single particle, Ω 1 , interacting with the diffuse interface of a non-interacting wall, Ω wall . The red square, Ω 1wall,max , indicates the position of the maximum overlap of the respective diffuse interfaces. The blue area, Ω wall* , indicates the area where Dirichlet boundary conditions are applied.
surrounding unstressed particle surface. The stressed particle experienced 0.55 μm vertical displacement.
The bulk, c-gradient, η-gradient and elastic energies were calculated during the simulations. The respective energy densities within the contact as a function of time are plotted in Fig. 7(a) for the 0.4 Nsingle particle simulation. The elastic energy within the contact dropped exponentially as the contact patch grew in size. Eventually an asymptote was reached where the contribution of the elastic energy to the overall energy state within the contact was dwarfed by the other energy terms.
A global energy analysis at the end of the single particle simulations with 0 N and 0.4 N force was also performed. Since the phase field method is used to describe an irreversible thermodynamic process, it is expected that the overall total energy of the system decreases as the system attempts to find a minimum energy configuration along the space-energy hypersurface; this was found to be true. The difference in respective energies between the 0 N and 0.4 N simulations is plotted in Fig. 7(b). From the graph it is evident that the incorporation of elastic energy resulted in a net positive difference in the bulk energy and a net negative difference in the respective gradient energies. This is due to the smearing action of the elastic energy and is attributed to the generation of new surface which is to be expected as the particle surface changes shape.

Two particle tests
A comparison of the resultant contour of two particles sintering at 0, 0.4 and 0.8 N uniaxial body force after 60 s is shown in Fig. 8. The calculated maximum Hertzian contact stress for these loads is 0, 64 and 102 GPa respectively which fall within the estimated range experienced by the diamond grains during cold-compaction.
The lengths of the necks between the particles, X, at the end of each of the simulations were measured. The results for each respective body force is shown in Table 2. In theory, the neck size grows in time at a rate proportional to the cube-root of sintering time ( t 1 3 ) [11]. At equilibrium, the classical asymptotic limit for the neck size between two sintering particles is determined by the particle diameters and the ratio of the grain boundary to solid-liquid surface energies [1]. As a check, extremely long simulations ensured the neck widths converged asymptotically to their theoretical values.
Since the sintering of the particles in this study was performed within a time frame much shorter than that required to achieve equilibrium, the neck widths were considerably shorter than their equilibrium values. Nevertheless, it is evident that the kinetics of neck growth were significantly affected in the simulations. At a load of 0.8 N, a 34% increase in grain boundary length was observed. The increased connectivity in the two grains would result in increased strength between the sintered particles.

Multi-particle tests
Obtaining the starting particle positions for a multi-particle simulation is possible using standard packing algorithms. However, computing the inter-particle forces for particles confined in a box, especially for non-monomodal distributions, is not trivial. For this reason, it was necessary to use the LAMMPS granular discrete-element solver. The particle's position, radius and inter-particle normal forces were exported. The results were then post-processed to show the inter-particle normal forces as connecting lines between their centres with the thickness of the line proportional to the magnitude of the force. The resulting force-chain for one of the particle packing arrangements is shown in Fig. 9. From the figure it is evident that the compact possesses a complicated force distribution. Most of the particles contribute to    Fig. 9. The resulting force-chain distribution for particles confined under pressure of 8 GPa. The thickness of the lines connecting the particle indicates the magnitude of the normal force.
supporting the confining pressure of the walls by the formation of a particle skeleton. Due to zero gravitational forces, some of the particles do not make contact with other particles and float about in free space while others, although touching, experience no significant inter-particle forces.
The effect of an 8 GPa confining pressure on the final microstructure after 10.65 s for the sintered particles is demonstrated in Fig. 10. Qualitatively, the liquid pores appear fewer and smaller with thicker necks between the particles for the high pressure simulation. An increase in particle coordination from 3.26 to 3.42 for this simulation run was also observed.
The total length of the grain boundaries at the end of each simulation were measured. The relative difference in total grain boundary length between the pressure-less case and each pressure is shown in Fig. 11. The fit to the graph is a second order polynomial indicating a square dependence of grain boundary length on pressure. Thicker necks and increased densities have been shown to give sintered compacts improved stiffness, hardness, fracture toughness and compressive strength [39][40][41].
The starting solid density for the three different initial configurations ranged between 80-82%. The changes to the solid density of the sintering compacts were computed as a function of time.
The densification curves for the typical diamond sintering pressures are plotted in Fig. 12 as the difference in density between the high pressure simulations and the pressure-less case. For the 8 GPa runs, the solid density increased by around 4% for the duration of the run. The results show that even in the absence of plastic deformation, elastic energy at the inter-particle contacts is responsible for significant densification observed during diamond sintering. The ever increasing error bars are indicative of how small changes in the initial particle arrangement can lead to different outcomes in the final microstructure. Kingery [7] provided a theoretical model for the shrinkage occurring during liquid phase sintering. Upon wetting of the particles, the surface tension of the liquid, γ lv , results in a negative pressure acting on each of the remaining gaseous pores. The pore pressure is equivalent to placing the entire system under a hydrostatic pressure which in turn results in appreciable forces acting on the particles. Derivation of the equation was modified (see Appendix) to include the external pressure as the driving force for densification rather than the surface tension of the liquid:   11. A plot of the relative increase in total grain boundary length between the pressure-less case and each confining pressure. The error bars indicate the standard deviation for measurements across the three different starting configurations for the multi-particle simulations. where ΔL is the change in length, L 0 is the original length and p 0 is the externally applied pressure. Q is a measured geometric coefficient equal to the relative difference in grain boundary length at the respective pressure, which has been derived from fitting the data shown in Fig. 11 and k is the constant of proportionality as reported in the appendix and was calculated to be 17.5, C the solubility of the solid in the liquid phase, V m is the molar volume, R the universal gas constant, T the temperature and r the particle radius. To make the model correspond to the phase field model, C was taken as 0.5, δ L is the liquid film thickness which was taken as the diffuse interface width and r was taken as the mean starting grain size of 20 μm. The analytical solution to Eq. (19) for each of the respective pressures is included as the dashed line in Fig. 12. Generally a good agreement in the shapes of the densification curves between the phase field and theoretical model was observed. At the beginning of the run, the theoretical model displays slightly increased densification rates. This observed difference is likely due to the delay, t init , used in the phase field models that resulted in partial bonding and larger contact patches initially. The delay, while unphysical, is necessary as it not only allows the diffuse interfaces to develop adequately, but also provides numerical stability for the model in the early stages of sintering where point-like contacts result in enormous elastic energies.
In reality, scanning electron micrographs reveal that the starting diamond particles are not spherical; grains are highly prismatic and angular. While the proposed phase field model is intrinsically capable of handling non-spherical particles, a more advanced discrete element model would be required to calculate inter-particle forces for the irregular surfaces. A complete model for the sintering of diamond would also require the implementation of crystal plasticity and quantification of plastic strain energy which is beyond the scope of this article. Further developments in these areas of high pressure sintering models are needed and more advancements will be presented in future contributions by the authors.

Conclusion
The implementation of an efficient contact mechanics algorithm into an elastic phase field sintering models is presented. Contact detection and pressure distribution in the contact were approximated by the overlap of the order parameter fields. Advection fields were used to maintain contact between deforming surfaces. The model was demonstrated by applying it to three different scenarios. The model showed that significant contact flattening occurs when a single particle was pressed onto a non-interacting rigid wall at loads typically experienced during diamond sintering. An energy analysis for the single particle test revealed that the stored elastic energy was transformed into the generation of new particle surface. Typical inter-particle contact pressures for the LPS of diamond particles resulted in grain boundary lengths up to 34% larger than the pressure-less case. Such microstructural changes are expected to have significant changes in the mechanical properties of a sintered compact. For the multi-particle simulations, the initial packing configurations and inter-particle normal forces were obtained through discrete element methods. Densification kinetics for the elastic phase field model were in good agreement with theoretical predictions of shrinkage during LPS. At a confining pressure of 8 GPa, a 4% increase in solid density was experienced by the sintering compact for the duration of the run. A new formula was derived based on mechanistic considerations that, starting from equivalent theories developed for low-pressure liquid phase sintering, include the effect of contact stresses to estimate solid densification rates in multiparticle HPHT sintering processes; it was shown to provide very good predictive capabilities for the range of conditions and the system considered in this study.
Compressive stresses at the contact result in an increase in chemical potential, μ, and activity, a, as follows: where ΔP is the compressive stress at the contact. The increased activity in the contact results in a difference in solubility between the contact and surrounding particle surface such that the diffusive flux, J, away from the contact points equals The volume of material removed, V, is equal to: where h is the inter-penetration strain of the grains. For continuity, the volume of material removed must be equal to the flux away from the contact point: Combining Eqs. (21) and (22) we arrive at the relation: ( 1 ) 2 rhdh dt .
By approximating the exponential with the first term in its series expansion and rearranging we obtain: The contact pressure is related to the remote external pressure, p 0 by: where Q is the relative difference in contact path length as a function of external pressure and k is the constant of proportionality. The product kQ 2 is thus a geometric factor that relates the contact pressure by the external pressure. Substitution of Eq. (29) into equation (28)