Unravelling densification during sintering by multiscale modelling of grain motion

The resulting microstructure after the sintering process determines many materials properties of interest. In order to understand the microstructural evolution, simulations are often employed. One such simulation method is the phase-field method, which has garnered much interest in recent decades. However, the method lacks a complete model for sintering, as previous works could show unphysical effects and the inability to reach representative volume elements. Thus the present paper aims to close this gap by employing molecular dynamics and determining rules of motion which can be translated to a phase-field model. The key realization is that vacancy absorption induced motion of grains travels through a grain structure without resistance. Hence the total displacement field of a green body is simply the superposition of all grains reacting in isolation to local vacancy absorption events. The resulting phase-field model is shown to be representative starting from particle counts between 97 and 262 and contains the qualitative correct dependence of sintering rate on particle size.


Introduction
The sintering process is an important step in many materials processing routes, from the humble coffee cup to complex applications such as solide oxide fuel cells.Especially for the latter, the properties of the material are critical for the application.The materials properties depend on the microstructure of the material and thus predicting the microstructure is of large importance.While analytical models provide a first avenue to microstructural predictions, the necessary simplifications in geometry and other aspects often prevent quantitative predictions.Simulations offer another avenue which is less restrictive in the necessary simplifications.In recent decades, field-resolved simulation methods such as the Monte Carlo method [1,2] and the phase-field method [3][4][5][6][7] have attracted much attention for the simulation of sintering.It was recently shown [7] that the most popular phase-field model of sintering showed a significant influence of the green body size on the densification behaviour.Since this runs counter to physical intuition and experimental evidence, the present paper aims to eliminate this problem.
The goal of this paper is twofold: First, qualitative rules of motion during sintering are determined, as to allow the prediction of grain motion and length changes.This will be achieved by conducting molecular dynamics simulations in geometries which allow the isolation of the relevant processes.Second, these rules are translated into a phase-field model in order to allow for simulations of arbitrary geometry and scale.Following this translation, the phase-field model is compared against other phase-field models of sintering before being used to determine representative volume elements.

Molecular dynamics
Molecular dynamics (MD) is a method in which the dynamics of individual atoms under the influence of an interaction potential can be simulated.The individual atoms are assumed to follow Newton's laws of motion, with the interaction potential determining what kind of material is being simulated.Sintering has previously been investigated with MD by various authors [8][9][10][11], but with a general focus on identifying the sintering mechanisms at the nano-scale rather than determining rules of motion for coarser spatial methods.The present study is somewhat similar to Hawa and Zachariah's work [8,9], in which they investigated how a chain of amorphous Si sintered and considered the influence of chain length and particle size and how these affect the velocity distribution and sintering time.
For the present investigation LAMMPS [12] is used to conduct MD simulations.As a model material copper is employed by using EAM potential developed by Foiles et al. [13].A comparison with more recent copper EAM potentials was conducted.While the quantitative results did change, the qualitative trends did not and the employed potential was much faster to calculate.The timestep employed within the MD simulations is generally 0.004 ps.
The primary goal of the following simulations is to predict the sample length change ∆L during sintering, parametrized by variables accessible on the mesoscopic phase-field scale.For this purpose the geometry depicted in Fig. 1 is developed.It contains a chain of rectangular cuboid grains of alternating, different orientations, with both ends of the chain being free as to allow movement.On each of the grain boundaries a pore may be placed, which then vanishes during the sintering process, which in turn induces movement of the grains.By employing cuboid grains extending to the periodic boundary, the grain rotation as commonly found in MD simulations of sintering [11] is largely suppressed.This makes the tracking of the pore region much simpler.It also makes the calculation of rigid displacements simpler because no rotational displacement needs to be removed.
The following results will be mainly based on placing spherical pores of radius ≈ 1.2 nm, but the qualitative trends of the results do not change when the size is varied or when a cylindrical pore is placed.The visualization of the results is done with OVITO [14] and matplotlib [15].The view of the simulations will be from the positive z direction as indicated in Fig. 1, unless mentioned otherwise.
The system is prepared as follows: First, a bicrystal of size (nL0, mA0) is prepared with a base length L0 = 141 Å and base area A0 = 1306 Å2 with free surfaces along the [100] axis (x-axis) of the simulation cell.n, m are positive integers, which axis.These GBs are chosen because they were easy to directly construct within LAMMPS.After the atoms are set, a conjugated gradient minimization at 0 K is conducted, followed by an constant number of atoms N , pressure p = 0 and temperature T (NPT) ensemble heating run from 1 K to the target temperature T = 700 K over 200 ps, followed by another 320 ps of equilibration at constant target temperature T = 700 K.This condition of p = 0 and T = 700 K will also be held for the rest of the paper.The system is then copied and shifted until the desired chain length is reached and equilibrated for a period of 400 ps to allow the system to relax the newly constructed grain boundaries.Once the system is ready, regions on the grain boundaries are defined and the atoms removed in order to place pores.By counting the number of atoms within these regions during the simulation, it is possible to obtain an estimate of the number of absorbed vacancies.Furthermore, the center of mass (COM) of the individual grains is tracked as to allow calculation of grain displacements and the total length change.The grain displacements are calculated directly by subtracting the center of mass xi(t) of grain i at time t from that at time t = 0.The total length change in a direction is assumed to be the length change of the vector connecting the center of mass of the first and last grain in the chain.
It is taken to be positive for a shortening of the vector.The sintering simulations on structures with pores are generally run in increments of 8 ns.If a pore is observed to have vanished, the simulation run may be terminated early.Not all simulations are run to complete pore elimination as the goal of the study is to find rules which are also applicable during the process and not only after it.A typical simulation result, starting with a cylindrical pore, is shown in Fig. 2 with the displacement per atom on the bottom.As expected, the pore vanishes over time and its vanishing is correlated to an inwards movement of the free surfaces.The displacement per atom is observed to be largely homogeneous within the grains, with large deviations only found close to interfaces, lending credence to the common assumption of rigid-body motion during sintering.The displacement in x, the shrinkage direction, generally directly scales with the pore size and is influenced by the area of the grains.While there were differences for the displacement in y, they generally did not show any monotonic relationship to pore size or grain geometry.In the shown simulation a gradient for the y displacement seems to exist, though this did not always manifest in the other simulations.Interestingly, the displacement in z did not show a sign difference between the grains.This might be due to the tilt axis being the z axis, though further research in this direction should be conducted.For the present paper we shall focus on the displacement in the shrinkage direction and assume that it is described by a rigid-body motion.
Next, the influence of grain geometry on the length change is investigated in a bicrystal.In Fig. 3 the results are shown for various grain lengths and areas.In this and the following plots, only every 25th point is marked unless mentioned otherwise; the line always indicates all data points.It is easily observed that the results are clustered by the grain area, whereas the grain length has no consistent influence on the results.The small influence of the grain length is likely due to excess stress caused by finite size effects.This also causes the pore to vanish at different times.Note that there is a linear relationship between vacancy absorption and displacement.By only plotting the system response over the number of absorbed vacancies, differences in kinetics can be removed from the problem.
As the present setup does not allow the differentiation of the GB area (GBA) from the grain cross section (GCS), a second series of simulations is performed.In these the area around the grain boundary is largely removed, effectively decoupling the GB area from the grain cross section.This is achieved by leaving only a rectangular area of size GBA in a region of length 50.61Å around the grain boundary.Thus new surfaces are introduced to a thin region between the grains, which need to be relaxed.This relaxation is done for 16 ns before a spherical pore is placed on the GB, followed by another run for up to 16 ns.It should be noted that the newly generated surfaces can also act as vacancy sink/source, whose contribution to the supposed absorbed vacancy count is not easily accounted for.Thus the relationship between displacement and vacancies shown in Fig. 4 will differ from that of the previous results.More specifically, more vacancies are being absorbed than plotted which also causes a larger displacement 1 .The squares and triangles show results for which GCS = GBA, whereas those in which GCS ̸ = GBA are marked with crosses and dots.The length change seems to be mainly influenced by GBA instead of GCS.
Hence we may formulate the first two rules of motion: The length change of a bicrystal sample due to vacancy absorption is antiproportional to the grain boundary area of the sample.The length change of a sample due to vacancy absorption is proportional to the number of absorbed vacancies.
Next we shall consider the influence of adding more grains, and hence grain boundaries and pores, to the system.The results of investigating chains with up to 8 grains are depicted in Fig. 5.For consistency, two grain boundary areas were considered and the first rule of motion is confirmed again.If one follows the line described by a smaller simulation, one can then reasonably predict the length change observed in the larger chain.Hence the displacement induced by the vacancy absorption on each grain boundary tends to be transported along the whole chain without any resistance.This can be interpreted as a kind of superposition property inherent in the solution, i.e. for a system containing multiple vacancy absorption sites, the total solution is the sum of all single vacancy absorption site problems.This is verified by running seven simulations for a chain with 8 grains, but with only a single pore being placed on a different grain boundary each time.The seven individual solutions are then added together and compared against the full solution, shown in Fig. 6.The calculated solution and measured solution match well, hence we may define the third rule of motion: The total length change of a system is determined by the superposition of all length 1 Using the later model eq.( 3) and presuming that the volume change is caused by both area and length changes ∆V = A∆L + L∆A yields a bit more than twice as many vacancies being absorbed, which fits quite well with these results.for two grain cross sections (GCS) and four grain boundary areas (GBA).Once GBA is independent of GCS, it is the determining factor of ∆L.Due to the quick pore removal for GBA = 0.5A 0 every point is marked in this plot.
changes due to vacancy absorption sites.
Recapping these rules we have the following properties: 1. the length change is antiproportional to the grain boundary area and independent of the grain length 2. the length change is proportional to the number of absorbed vacancies within the system 3. the length change of a system with multiple vacancy absorption sites is the superposition of the individual length changes In the following we shall shortly derive a model which contains all these properties.In spirit it is rather similar to DeHoff's theoretical developments in the 1980s [17,18], but without requiring the grain structure to be decomposed into a space-filling cell structure.First, based on item 2 we assume that each vacancy contributes a certain volume change dV proportional to the atomic volume Ω, with dN being the number of vacancies which have just been absorbed.Second, we assume that this volume change is due to a rigid movement of the entire crystal lattice of magnitude dL, i.e. dV = AdL with the grain boundary area A, motivated by item 1 Thus one may write

ΩdN = AdL
(1) with the length change ∆L taken to be positive for a shortening of the sample.This model is verified by testing it not only against the already presented data, but also several grain boundary types, grain lengths, grain boundary areas, pore shapes, pore sizes, and numbers of pores.The grain boundary area A of each simulation is estimated based on the equilibrated area before the pore is placed.The atomic volume Ω is determined by observing the volume of an fcc lattice of copper atoms in a periodic box at T = 700 K employing an NPT ensemble with p = 0, resulting in Ω = 1.22×10 −29 m 3 which is close to the value given by [19].The number of absorbed vacancies is known     The presented model so far seems to only describe the total length change ∆L and not the motion of individual grains described by their displacements uα, or equivalently velocities, as required in a field-resolved method as the phase-field method.In order to resolve this, consider the implication of a length change ∆L in an effectively onedimensional bicrystal αβ: Since both grains move as rigid-bodies, the length change is the sum of the individual displacements and thus is given by which is also the displacement jump ∆u αβ across the grain boundary.The superposition property item 3 is now exploited to enforce this simultaneously for all grain boundaries, leading to which is a usually overdetermined linear system of equations.The matrix C consists of rows with zeros and only one +1 and −1 each and acts on the unknown N grain displacements u = (u1, . . ., uN ) T .The sign of the entries is determined by the onesided grain boundary plane vector nα = −n β in the laboratory frame, with nα being normal to the αβ grain boundary and pointing out of the α grain.The right-hand side vector is determined by eq. ( 3) for each of the B grain boundaries.This system may be solved e.g. in a least-squares sense, with the conservation of momentum accounted for afterwards by subtracting the mass-weighted average displacement from the solution u.
For the special case of a linear chain of grains, there are B = N − 1 grain boundaries.Adding conservation of momentum to the system of equations makes the matrix C square and since the individual rows are linearly independent2 , it also is of full rank  6) and the observed data for two simulation states.A good match is observed for both.
and thus uniquely solvable.Furthermore, this formulation only accounts for motion due to vacancy absorption.If other processes inducing displacement occur, e.g. grain boundary sliding, then these will require a separate treatment.As a first test of this, we use the time-dependent data of one MD simulation and input these as the right-hand sides of the system of equations.The contacts between grains which fill the matrix C are also determined from these.Since the chain is linear, including conservation of momentum makes it uniquely solvable.The system is solved by direct matrix inversion since it is rather small, with a comparison of the calculated grain movement and the observed grain movement shown in Fig. 8 for two simulation states.As the figure shows, there is a close match between the calculation and measurement, giving a measure of confidence in this approach.

Phase-field
In this section a new phase-field model will be described, following by a small-scale validation to ensure that the green body size effect is no longer present before large green bodies are calculated.

Phase-field model with advection
The model in the following is based on [7,20], with the advection velocity being calculated with a model based on the MD results.This new model in its entirety is dubbed MDi, as it is inspired by MD.The evolution equations for the fields are the same as in [7]: cannot be constructed from prior rows.
This represents the evolution of the N phase-fields ϕα and the chemical potential µ for one independent component, taken to be copper in the present paper.The phasefield tuple ϕ distinguishes the surrounding vapour (ϕV, V = 0) from copper grains of arbitrary orientation (ϕa, a > 0).The evolution of the chemical potential µ accounts for the species conservation via the concentration c and takes into account the effect of phase changes due to changes in ϕ.For further particulars of the terms the interested reader is referred to [7,20].The calculation of the grain velocities follows the ideas outlined in the previous section.This will be formulated in terms of instantaneous displacement u and number density of absorbed vacancies ∆n to be consistent with the MD section.The instantaneous velocity v is linearly related to u as v = u ∆t with the time interval ∆t.The concentration c is related to the number density n as n = Na Vm c with Avogadro's constant Na and the molar volume Vm.Each grain boundary αβ absorbs an amount ∆N αβ = GB ∆n αβ dV of vacancies in a time interval ∆t, with the density of absorbed vacancies ∆n αβ .This is related to the vectorial displacement jump in which the orientation of the grain boundary plane was taken into account by a similar approach as Wang [3], but employing normalized phase-field gradients representing the normal vector, i.e.
|∇(ϕα−ϕ β )| instead of just the phase-field gradients.Note that ⃗ n αβ needs to be chosen consistently with the later input to the matrix equation.It is always taken to be from the lower α index to the larger β index, defining a unique orientation for each αβ pair.The grain boundary region GB is defined to be the region in which g = ϕαϕ β ≥ gT holds, i.e.only the region where both grain phases have significant volume fractions.In the following the value gT = 0.14 is chosen.The grain boundary area A αβ is resolved by dividing the grain boundary volume V αβ by the equilibrium grain boundary width l0 = for the employed obstacle potential.The remaining 4ϕαϕ β and V αβ terms act as a weighted average to assign higher weight to regions which contain more grain boundary phase 4ϕαϕ β .
Equation (11) still contains an unknown, namely the number of absorbed vacancies.For this we assume that a grain boundary has a certain equilibrium number density of atoms n gb eq and that it relaxes towards this number density: which allows the calculation of the number density of absorbed vacancies ∆n = − ∂n ∂t ∆t.Since this describes a relaxation process, one can consider the term n − n gb eq as the driving force for this process -once it vanishes, densification via advection stops.Note that this assumes atoms and vacancies are both conserved quantities.Only the number of atoms is actually conserved, with vacancies and lattice sites being destroyed and generated during absorption to accommodate the volume change.However, this rough treatment suffices to show the capabilities of the model and is in fact quite standard in phase-field modelling of sintering [3][4][5][6][7].In spirit this is similar to Wang's model [3] but the relaxation time tr is identified explicitly here, which can in turn be determined via molecular dynamics.In the MD simulations tr was observed to depend strongly on the grain boundary orientation relationship; it is strongly related to the efficiency of a grain boundary at absorbing vacancies such as described in [21].In the following simulations, all grain boundaries are assumed to be of a (210)/[001] STGB (53.1 • ) type.Furthermore, n gb eq does not have the meaning of a "grain boundary equilibrium density" in this context.Specifically, if it is below the bulk density as one would expect based on physics, it will push apart grains instead of attracting them.This phenomenon has been investigated in-depth in [7,22].Based on the suggestions therein, n gb eq is calculated based on the observed average chemical potential on the particle surface μ, which should approximate the capillary pressure.This ensures that the resulting velocities are consistent with the free energy functional and the theoretical dihedral angle is recovered [7].Note that this allows a rather trivial addition of an external isotropic pressure to the system, as the capillary pressure can simply be shifted by the external pressure.
It should also be noted that properties such as ⃗ ∆u αβ and V αβ need to be tracked for each grain boundary individually and thus their memory and communication requirements scale as O(nN 2 ) if implemented naïvely, with the number of parallel processes n and number of phase-fields N .While this is not a problem for a few hundreds of grains, once thousands or tens of thousands grains are resolved this will dominate the memory and communication costs.This is resolved by only storing the actual contacts (thus being a sparse representation) and distributing it across all parallel processes.The message passing interface (MPI) is employed for the parallelization and updates to this distributed matrix are realized via one-sided communication.The details of this scheme will be published in a separate paper.
The displacement jumps are used to solve for the particle displacements u, for each direction separately, by building a system of equations Cu = ∆u (13) in which the matrix C is filled according to the connectivity determined during the simulation.The structure of the matrix is clarified by the following example: Consider a 2x2 grid of grains, depicted in Fig. 9.Each particle has two contacts, one along each dimension.These are always taken to be from the lower grain index a of ϕa to the higher one, i.e. we have the contacts described by the ordered set {C1,2, C1,3, C2,4, C3,4}.The size of this set is equal to the number of grain boundaries B in the system.A matrix C d of dimension B ×N is constructed per dimenions d, with a corresponding RHS ∆u d of size B. For each contact a row is added to the matrix, with only non-zero entries for those grains which are connected by this contact.The magnitude of the entries is always 1, but the sign is determined consistently with ⃗ n αβ in this direction, going consistently from the lower to the higher index.The right-hand side displacement jump is already in a vector form and therefore can be easily split into its components ∆u d .Thus eq. ( 13) can be written, for the x dimension, as with the sign on the left-hand side determined based on the grain boundary normal.
If the grain boundary normal has no component in a dimension, the sign within the matrix plays no role as the right-hand side will be zero.This effectively says that no relative motion occurs.Note that this matrix, although square, is singular, since e.g. the fourth row can be constructed by adding the first and second row and subtracting the result from the third row.Conservation of momentum is accounted for afterwards by subtracting the mass-weighted average displacement from each grain displacement.
The resulting system of equations is usually overdetermined and only in the special case of a linear chain of particles can be solved exactly if conservation of momentum is included in the system.However, the system can be solved in a least-square sense.This could be done with e.g. an Alternating Direction Method of Multipliers approach [23], which requires a collective reduction of N scalars per iteration.Alternatively, the matrix C and its transpose could be partially distributed and then be solved by using some Krylov subspace method (e.g. a block Conjugated Gradient Least Squares method [24]) which would need a collective reduction of a single scalar per iteration.However, both approaches lead to excessive communication time (on the order of entire field sweeps) and thus an approach suited to the present problem is developed.Typical solutions of the system were sought by generating packings, from which eq. ( 13) was determined while assuming the right-hand side is given by a normal distribution with its mean and standard deviation σ.A mean of µ = 1 is generally employed with variable standard deviation σ.A spatial dependence can be included by simply adding a function of particle position to the random sample generated by eq. ( 14).
The system is solved employing the Least Squares via QR factorization (LSQR) [25] algorithm.
It was observed that the spatial distribution of the right-hand side tends to determine the shape of the solution.If it is assigned randomly without any spatial dependence, a mostly linear function of position is observed, with local inhomogeneities.If a linear dependence on the position is added, the displacement field becomes a quadratic function of position.This implies that the solution of the system basically integrates the spatial distribution of right-hand sides.This can also be seen from the structure of the matrix C: Each row effectively represents a finite difference formula, with the right-hand side giving the slope, i.e.C is a differentiation operator.Thus its generalized inverse is an integration operator.
Given that the same grain boundary type is assumed for all contacts, with similar initial neighbourhoods, it seems reasonable to assume the displacement field is given by a linear function.Hence one may approximate the full problem by replacing the particle displacements u by the relation with the known center of mass of each grain x, the known total center of mass xm and an unknown slope m.Hence only m remains to be determined, which can be done exactly in a least-square sense by employing the normal equations: which only requires a parallel reduction operation for p and q.This approach is compared against the full solution via LSQR.The error is evaluated with the root-mean-square error (RMSE) and relative RMSE (RRMSE) defined as for the solution vectors of the LSQR method and the linear fit ansatz.It is generally observed that the RRMSE is unaffected by the choice of mean µ, while the RMSE changes due to the scale in displacement.For regular packings the (R)RMSE is observed to scale mostly linearly in the standard deviation σ of the random distribution, with errors on the order of machine precision for σ = 0.The irregular packings employed later for the PF simulations generally show some non-zero error even for σ = 0, though starting from about RRMSE ≈ 0.16%.Even if σ is comparable to the mean of the normal distribution, RRMSE ≈ 1.6% and thus still quite acceptable.A visual comparison for the effect of σ on the displacements in a 3D packing containing 3445 particles is shown in Fig. 10.As can be seen, the solution shape and scale are always well-preserved.What the fit ansatz obviously cannnot match however is the local variation of absorption activity, modelled by the random distribution of displacement jumps.The interested reader is referred to the Supplementary Material wherein the code employed for this test is published in full, along with the irregular packings employed later for the large-scale sintering simulations.As a final note, the specification of the vacancy absorption rate is the main weakness of the presented model, since it cannot be completely linked to quantities obtainable via MD.Hence a direct comparison to MD simulations would likely yield a significant mismatch in the temporal evolution.However, if an improved model for the vacancy absorption rate is developed, it can be included easily into the current approach.This is due to the model effectively being split into a kinematic part eq. ( 13) and a dynamic part eq. ( 12) which can be changed independently.

Parameters and data evaluation
The scales employed are listed in Table 1 and the materials parameters in Table 2.These are the same as in [7] except for the newly introduced atomic volume Ω and the relaxation time tr.The relaxation time is determined by running MD simulations for a (210)/[001] STGB (53.1 • ) in which atoms are randomly removed from the grain boundary and observing the time it takes until the atom count within the grain boundary has stabilized.Several such simulations were run and the order of magnitude for the relaxation time then used for the value of tr.This is more of qualitative approach, but tr behaves similarly to the stiffness κ in the classical rigid-body motion (RBM) model of Wang [3]: In the classical model, the RBM velocity scales as v ∝ κ(n − n eq gb ) (22) whereas in the present model it is rather i.e. the velocity is proportional to κ, but inversely proportional to tr.Shi et al. [26] could show that once κ is sufficiently large, the resulting advection velocity does not change upon a further increase in κ.Hence the proportionality is only applicable to a certain limit, after which other processes control the velocity.Since tr divides the driving force n − n eq gb whereas κ multiplies it, the same behaviour applies, but reversed: Once tr is small enough, making it smaller will not change the velocity.One can think of this as saying that the problem ought not to be controlled by the absorption rate and hence the absorption rate should be appreciably faster than the slowest other process.The time step is determined by calculating stable time step widths within the explicit scheme and employing the minimum to ensure stable time integration as described in [7], with the table only listing the largest stable timestep with zero advection velocity.
The evaluation of the data is the same as in [7] in terms of strain and density.The contact number later employed in the three-dimensional simulations is calculated with the package cc3d's [27] function contacts on the phase-field label field and the simply counting the pairwise occurrences of labels.The surrounding vapour phase is treated as background, with the phase-field label field being defined pointwise as the label which has the highest phase-field value.Since surface particles always have some missing neighbours, including these would induce a particle size and simulation size bias in the coordination number.These are excluded by determining the bounding box of the green body, then shrinking it by 32 nm in each direction and only averaging over particles contained in this shrunken box.In the 3D simulation this excludes at least two particle layers, removing the surface effect.
The simulations in section 3.3 are conducted locally while employing GNU Parallel [28] for efficient job management.The three-dimensional simulations in section 3.4 are calculated on the Hawk supercomputer at the High Performance Computing Center in Stuttgart.Hence the processor employed for the three-dimensional simulations is the AMD EPYC 7742, with 64 cores running at 2.25 GHz.A single core performance of 9.4 GFLOP/s, i.e. 26.1 % of the theoretical peak, is achieved.

Validation
In this section the system size convergence of the present model will be investigated.In effect, this tests whether the superposition rule of motion from the MD simulations has been transferred successfully.As will be shown below, the results for a two-particle model are virtually identical between the present model MDi and a grand potential model including advection, i.e. the model (ADV-µ) of [7].Hence the accordance with classical theory in terms of neck growth and approach of centers, as well as Herring's scaling laws for a two-particle model as shown in [7], are transferable to the present model.Thus the strain in a system of increasing size will be investigated, first in a particle chain as suggested by [22], then in a rectangular grain geometry.Simulations are run for both models, with the number of particles in the chain given by n ∈ {2, 4, 8, 16, 32}.The simulations for the same geometry are all run to the same simulation time te.The strain is evaluated based on the movement of the barycenters of the first and last particle, i.e. e = L(t)−L(0) with L(t) = x m,last − x m,f irst and xm, * being the the barycenter of the first/last particle.The geometries considered in this section can also in general be used as a relatively cheap benchmark geometry for determining whether size-independent densification is captured by the model.
The strain over time for a chain of circular particles is shown in Fig. 11(a).While model MDi seems to converge at around 16 particles in a chain, model ADV-µ still shows a large change at this particle count.Let us first consider how model ADV-µ fails to converge by looking at the velocity distribution in Fig. 11(b).The velocity is plotted over the relative barycenter Xi,r = x m,last −x m,f irst which allows the compact viewing of chains of arbitrary physical length.The absolute scale of velocity reached is similar for both models, but while model MDi always yields a linear function by design, model ADV-µ tends to produce curved velocity profiles which directly imply inhomogeneous densification.While a small amount of inhomogeneity is to be expected due to the discussion below, this should drop off rapidly from the outermost particle.
One might ask why model MDi does not instantly converge in this case, given that  The dependence on the strain becomes negligible at around 16 particles for the MDi model but is substantial across all investigated particle counts for model ADV-µ.The velocity distribution shows the inhomogeneous densification (non-constant velocity gradient between the particles) from which this lack of convergence originates.the model without the phase-field information did so?This is due to the chain not actually being homogeneous and only becoming so at a sufficient number of particles.The simulation geometry at the end of the simulation is depicted in Fig. 12, showing that the chain ends take up a different shape than inner chain particles.This shape difference is simply due to the inner particles being restricted from free movement, whereas the outer particles can easily adjust.This difference will also eventually lead to grain growth.
The chain ends being different can be somewhat mitigated by employing the rectangular grain geometry from the MD simulations and placing pores on the grain boundaries.It will not be fully mitigated, as the absorption rate still depends on the average surface chemical potential, which is different for end grains and for inner grains.However, as Fig. 13 shows the convergence is sped up with this geometry for the MDi model, but the lack of convergence for ADV-µ becomes even more obvious.It should be noted that the MDi model eliminates the pores on the GBs at similar times, whereas model ADV-µ eliminates them step by step from the outer parts of the chain, with videos showing this being deposited with the Supplementary Material.This is also the reason for the different velocity magnitudes between the two models: Up to the point of pore elimination, the size of the pores in the simulations with model MDi is roughly comparable to that of the outermost pore of ADV-µ at the same time.Thus the pores for simulations with model MDi are on average smaller than for model ADV-µ which implies a larger driving force for vacancy absorption.
Based on these results, one can conclude that at about 16 particles in a onedimensional chain the model MDi becomes representative.Presuming that this result extends to three-dimensions, one would expect representative simulations to start from about 16 3 = 4096 particles.However, it is likely that the end geometry problem makes this a significant overestimation, and thus we proceed to three-dimensional simulations to test this.

Large scale three-dimensional simulations
The model will now be employed to simulate three-dimensional green bodies in order to determine representative volume elements (RVEs) by verifying that the densification is independent of the green body size.For this, the packings as described in [7] are employed: A voxel domain of size N3 v , Nv ∈ {200, 400, 800} is filled based on a packing generated with the discrete element method.In order to minimize boundary effects, it is ensured that there are at least 15 voxels between the outermost edge of a particle and the global boundary.All fields employ zero-flux conditions on the global boundary.The voxelization happens with a fixed number of voxels R ∈ {8, 12, 16} used to resolve the particle radius, allowing the investigation of particle size effects with different simulations.Given the employed non-dimensionalization scales and discretization, these correspond to particles of size R ∈ {8, 12, 16}nm.Table 3 lists the number of grains Ng within each combination of (N 3 v , R), as well as on how many cores the simulations were run and for how long.Figure 14 shows one of the structures at different simulation times; videos of the entire process are deposited with the Supplementary Material.The left side shows a simple visualization of the entire green body, with the right side showing a fracture surface generated with a (011) plane and removing any grains beyond this plane.The mesh visualized here is based on a cellwise maximal value of the phase-field vector excepting ϕV.Contour levels l > 0.5 of this field cause an etching-like effect to appear starting from the highest order junctions 3 .A contour level of l = 0.6 is employed which entirely reveals the triple lines in the three-dimensional structure.In any case, both visualizations show the macroscopic densification of the body, with the fracture surface view also showing that the grains transform from spheres to polyhedra.
The evolution of the density is shown in Fig. 15, with the present model as well as the results of [7] for a particle size of 12 nm.It can easily be seen that the present model shows highly similar density evolution over all the considered green body sizes.Only for R = 16 there is a slight effect of green body size from 200 3 to 400 3 , with the following simulations being quite similar again.It is likely that rather than the green body size, a certain minimum number of particles should be contained within a three-dimensional packing, with that threshold lying between 97 particles (200 3 , R = 16 nm) and 262 particles (200 3 , R = 12 nm).The particle count likely acts as a proxy variable for whether the geometry is sufficiently homogeneous.Thus green bodies with a polydispserse grain size distribution will likely need a larger number of Table 3: Initial grain counts N g for the employed packings as well as the number of cores C employed and the total runtime T .The longer runtime of larger particles is due to these being run for longer to achieve comparable densities.particles to be representative, but this will be left to future work.Furthermore, it is easily seen that densification progresses more slowly with larger grains.Let us shortly revisit why the present model does not fail to reach a RVE: The necessary requirement for densification is for the divergence of velocity ∇•v to be negative between volume elements.In advection models only employing nearest-neighbour interactions such as [3] and models based on it, the velocity of particles only depends on their immediate neighbours.Since within the green body proper, a grain's immediate neighbours will be similar, the neighbouring volume elements will be similar, with the exception of those volume elements containing the green body boundary.In contrast, in the MDi model the velocity of a single particle depends on all particles via solving eq. ( 13).Thus there is nothing forcing neighbouring volume elements to be similar.The simplification of using a linear ansatz for the particle displacement of course forces a constant ∇ • v between neighbouring volume elements.Given that solving the complete system for spatially uncorrelated absorption does result in the particle displacement being a linear function of position, this is quite justified.24) are depicted.The coordination number rises monotonically with density, but shows a systematic deviation from German's relation.However, the slope of the curve is highly similar, as is shown by also plotting a vertically shifted version of German's relation.
The present results also allow to qualitatively test whether the geometries assumed in intermediate stage sintering are found.Coble's classical model [30] assumes a tetrakaidecahedron, i.e. a solid with 14 faces, for the grain shape.Thus the number of neighboring grains, also called the coordination number Nc, should tend towards 14.A more quantitative relation is given by German [31], based on a fit of literature data: with the fractional density f which is equivalent to the present usage of density.The average coordination number of all grains is plotted over density in Fig. 16, showing a monotonic increase of coordination number with density.After correcting for surface effects, a coordination number of 14 is reached at around 99 % density.This includes the effect of many small contacts, which might not be detected easily in experiments.Hence there is a systematic deviation from German's fit, but the slope is quite comparable.If the fit is shifted vertically by a value of 1.1, which is roughly the difference in starting coordination number at 60 % density, then a quite close match is observed.Conversely, the simulation data could also be filtered to exclude contacts with small surface area, also resulting in a reasonable match for higher densities.This is explored in the Supplementary Material, as choosing any one minimum surface area is quite arbitrary.In any case, based on both the experimental fit due to German and the present results, the tetrakaidecahedron shape assumption does not hold in intermediate stage sintering.It can however hold in the final stage.

Conclusion
In the present paper molecular dynamics (MD) is employed to investigate the densification behaviour of a chain of grains.The grains are observed to move largely as rigid bodies, i.e. the atomic displacement within a single grain is largely uncorrelated to atomic position.Three rules of motion for this displacement are found: The dis-placement is proportional to the number of absorbed vacancies, antiproportional to the grain boundary area and superimposable if multiple vacancy absorption sites are present.These rules were used to construct an analytical model which agreed well with the MD results.Following this, a previously published phase-field model (ADVµ) was extended with this new model for calculating velocities (MDi).The previously published phase-field model and the new model are then compared in terms of their strain evolution within a linear chain geometry.For the MDi model, the strain as a function of time is observed to become independent of the number of grains between 8 and 16 grains in the chain, depending on the particulars of the geometry.Model ADV-µ did not converge, as previously shown by [7].Finally, the model MDi is employed to sinter large-scale three-dimensional structures to determine representative volume elements.It is found that between 97 and 262 particles are necessary for densification to become independent of the green body size.Furthermore, the qualitative correct influence of particle size is included in the model, with green bodies consisting of larger particles sintering more slowly.Finally, reasonable agreement with a model linking the coordination number of grains to the density could be shown.
A future work exploring the applications of the model to pressure-assisted sintering and the investigation of concurrent densification and grain growth is planned.
(a) initial state showing all atoms (b) intermediate state showing HE atoms (c) final state showing HE atoms (d) displacement in x (e) displacement in y (f) displacement in z

Figure 2 :
Figure 2: Simulation results exemplarily depicted for L = 2L 0 , A = 2A 0 and a cylindrical pore.The color in the top row indicates the local orientation in the viewing direction (Z), calculated via polyhedral template matching[16], allowing to distinguish the grains.In the last two two images, only high energy (HE) atoms (potential energy of > −3.2 eV) are shown, revealing the interfaces.The displacement of the atoms in the x (d), y (e) and z (f) directions after the pore has vanished are plotted over the original atomic position on the bottom.The displacement range is fixed to [-3, 3] Å since the large displacements on the surfaces would obscure the behaviour within the grain.The inner parts of the grains are generally homogeneously displaced, with the interfaces showing large deviations.Similar plots are obtained for all other simulations.

Figure 3 :
Figure3: Length change ∆L of a bicrystal containing a (210)[001] STGB for various geometries, over time and over the number of absorbed vacancies.The length change is observed to be strongly dependent on the area, with only a weak dependence on the total length.The pore elimination time, roughly given by when ∆L becomes constant, is slightly dependent on total length.A linear dependence of ∆L on the number of absorbed vacancies is also observed.

Figure 4 :
Figure4: The length change ∆L of a bicrystal containing a (210)[001] STGB for two grain cross sections (GCS) and four grain boundary areas (GBA).Once GBA is independent of GCS, it is the determining factor of ∆L.Due to the quick pore removal for GBA = 0.5A 0 every point is marked in this plot.

Figure 5 :
Figure 5: Total length change for two areas and up to 8 grains / 7 pores in the chain, with L = 2L 0 .A linear relationship between vacancies absorbed and the length change is observed.The results for longer chains roughly behave as if lying on the same line as for the smaller chains.

Figure 6 :
Figure 6: Comparison of a 8 grain chain with 7 pores and the solution via superposition.A good match for the total length change as well as the grainspecific displacements is observed.

Figure 7 :
Figure 7: The measured length change is plotted against the calculated length change based on eq.(3).A general match is observed, with the fitted line's slope (m) indicating a slight underprediction (5%) of the model.

Figure 8 :
Figure 8: Comparison of solving eq.(6) and the observed data for two simulation states.A good match is observed for both.

Figure 9 :
Figure 9: Example 2x2 setup of grains for clarifying the matrix structure.The number within the circles indicates the grain index.

Figure 10 :
Figure 10: Comparison of the full displacement solution calculated via LSQR and fit ansatz for two values of the standard deviation.Each marker indicates a particle's position and its resulting displacement in one dimension of the 3D packing.The shape and scale of the solution are always well-preserved by the fit ansatz.

Figure 11 :
Figure 11: Comparison of the MDi model (crosses) and the model of [7] (circles).The dependence on the strain becomes negligible at around 16 particles for the MDi model but is substantial across all investigated particle counts for model ADV-µ.The velocity distribution shows the inhomogeneous densification (non-constant velocity gradient between the particles) from which this lack of convergence originates.

Figure 12 :
Figure 12: The grain field, defined by 1 − ϕ V , at simulation end for 8 particles and model MDi.Yellow indicates the copper grains, dark indigo the surrounding vapour and green the interface between them.The phase transitions defined by the 0.5-isoline of the phases are represented with red lines.

Figure 13 :
Figure 13: Comparison of the MDi model (crosses) and the model of [7] (circles) when employing the rectangular grain geometry.While the MDi model converges at around 8 particles in the chain, model ADV-µ does not converge at all.The velocity distribution shows that model ADV-µ tends to produce nonlinear velocity profiles.

Figure 14 :
Figure14: Three-dimensional view of the 400 3 , R = 12 nm simulation for various times using Paraview[29].The left side shows the entire green body, with the right side showing a (011) fracture surface from the same angle.The dark lines delineating regions can be interpreted as grain boundaries and higher order junctions.The dark smudges on some grains are due to the dark areas being reflected via ray tracing and thus not actually part of the simulation data.Both macroscopic densification and the polyhedralization of the grains are evident.

Figure 15 :
Figure 15: Comparison of densification for the present MDi model and the results of[7] for a diffusion-only (DO) model and a model including advection (ADV-µ) based only on nearest-neighbour interactions.The results of[7] clearly have a strong dependence on domain size, which the MDi model eliminates once a RVE is reached.

Figure 16 :
Figure16: The Coordination number over density for all simulations as well as the relation of German eq.(24) are depicted.The coordination number rises monotonically with density, but shows a systematic deviation from German's relation.However, the slope of the curve is highly similar, as is shown by also plotting a vertically shifted version of German's relation.
Figure 1: Two-dimensional sketch of the considered geometry in the MD simulations.Grains of different orientation O 1 , O 2 are placed in a row, with pores on grain boundaries.The ends of the chain are free surfaces, with directions perpendicular to these being periodic.
will be varied.The main grain boundary orientation relationship which will be investigated in the present study is the (210)/[001] STGB.It is a symmetrical tilt grain boundary (STGB) with a misorientation angle of θ = 53.1 • about the [001] tilt axis, with the grain boundary plane being (210).Additional GBs for which the simulations were conducted are the (310)/[001] STGB and an asymmetrical tilt grain boundary with the left/right grain boundary planes being the (100) and (2 10) planes, with a 26.56 • rotation around the [001]

Table 2 :
Employed physical and numerical parameters for the simulations.