Revealing process and material parameter effects on densification via phase-field studies

Sintering is an important processing step in both ceramics and metals processing. The microstructure resulting from this process determines many materials properties of interest. Hence the accurate prediction of the microstructure, depending on processing and materials parameters, is of great importance. The phase-field method offers a way of predicting this microstructural evolution on a mesoscopic scale. The present paper employs this method to investigate concurrent densification and grain growth and the influence of stress on densification. Furthermore, the method is applied to simulate the entire freeze-casting process chain for the first time ever by simulating the freezing and sintering processes separately and passing the frozen microstructure to the present sintering model.


Introduction
Many products being used today, from the humble coffee cup to solar cells, undergo the sintering process at some point in their processing chain.While the precise materials properties of a coffee cup might not be particularly important barring liquid leakage, the electrical contacting of solar cells is important for their efficiency [1][2][3].Furthermore, since different materials are being co-sintered, solar cells also face the problem of warping [4], similar to solide oxide fuel cells.Finally, sintering is also usually used in the freeze-casting process chain to improve the mechanical properties of the structure.The porosity produced during the freeze-casting process is of two types: Macropores, aligned to the freezing direction, and micropores located in the walls of particles formed during the process.The goal of sintering this structure is to close the micropores, giving mechanical strength to the body, while keeping the macropores and thus achieving a highly porous, low-density but still sufficiently strong material.For example, this kind of material is of interest for tissue engineering scaffolds, in which the lack of small, interconnected porosity avoids the slowing down of cell proliferation into these regions [5].Furthermore, since the porosity can be made gradiated and is highly anisotropic, this also favors penetration of tissue into the freeze-cast scaffold [6].This requirement for controlling the microstructure, resulting materials parameters and dimensions is shared across many sintered products.Hence predicting the microstructure and its properties as well as the sintered dimensions is of great importance.Key properties of the microstructure are its density, grain and pore size as well as shape distributions, based on which further properties such as resistivity can be predicted.
In the past two decades, the phase-field method has been used by many researchers [7][8][9][10][11][12][13][14][15] to investigate the sintering process.More recent works have pointed out several weaknesses of these models [16][17][18], which were rectified subsequently [17][18][19].Based on these works, the present paper aims to elucidate several factors of interest in the solid-state sintering process: External pressure is often applied to speed up densification: Capillary pressure acts as a driving force for sintering and applying external pressure enhances this.Furthermore, the entire sintering process is fundamentally one of mass transport, with multiple concurrently active transport pathways such as grain boundaries (GBs) and surfaces.The ratio of their contribution plays a significant role, since densification is due to vacancy absorption on grain boundaries, with surface diffusion possibly acting counter to this mechanism.As elevated temperatures are usually employed for sintering, grain growth is another matter of concern.This not only increases the mean grain size, reducing e.g.mechanical strength, but also allows for porosity to detach from boundaries and be trapped within grains.Once this occurs, these pores cannot be eliminated without pressure on processing timescales.Hence the combination of grain mobility, as well as species mobility along interfaces, is another factor of interest.Finally, the mass distribution and pore anisotropy of a green body can cause anisotropic densification [20,21].
The effect of pressure will be investigated first, as it is a comparatively simple investigation.Since sintering and creep are closely related processes, a short excursion to creep is done as well to show the versatility of the model.Next, the materials properties of diffusivity and grain mobility will be varied, effectively representing sintering at different holding temperatures or with appropriate dopings.Special attention is paid to how different transport and growth mechanisms affect each other's influence on the microstructural evolution.In order to investigate the effect of green body geometry, a previously generated particle structure of the freeze-cast process is sintered, showing the experimentally observed anisotropic shrinkage of freeze-cast structures.The results of these investigations are presented first, with the employed phase-field model being succinctly described in Phase-field model and in further references therein.
The materials parameters are those of Table 1 unless mentioned otherwise.These parameters approximate copper at around 700 K. Copper is chosen as an example material as much experimental data of its sintering is available.Furthermore, since it is a metal, many of the complexities introduced when sintering ceramics, such as charged defects and space charge zones across GBs, can be safely ignored.Applications using sintered copper include microelectronic applications as well as foams used for heat exchangers, which can be produced via freeze-casting followed by sintering [22].Simple scalings of the diffusion and mobility values are taken to approximate different holding temperatures, such as to isolate the effects of the individual parameters.

Results & Discussion
Figure 1: The initial three-dimensional green body used for simulations, except for simulations concerning creep and freeze-casting, is shown.The overall shape is a cube, with spheres packed into the cube based on a larger packing from a discrete element method simulation; for details, see [19].The grains are shown as a copper-like material.
Based on prior work [19], the density evolution of 3D phase-field simulations of solid state sintering does not change substantially after more than 262 particles are contained within the green body, given that it starts with a uniform grain size and holds this size over the simulation run.If not mentioned otherwise, we shall employ a packing with 3445 spheres of radius r = 12 nm in a 400 nm 3 box for the initial conditions and observe how its evolution changes when various processing and materials parameters are changed.A small grain size is employed for computational efficiency, but as long as Herring's scaling law holds as shown in the supplementary material, the present results can be mapped onto coarser grain sizes.The initial green body of this packing is depicted in Fig. 1 representing the grain structure.The obtained densification should be independent of domain size based on [19], but the same need not hold for the grain size evolution since the number of remaining grains can become quite small.Additionally, the influence of green body geometry is analyzed by employing a freeze-cast structure obtained with the model of [23] and sintering it.

Sintering under stress
The first foray will consider the influence of stress on sintering, as it is quite simple but also important for densification.Pressure is often applied during sintering in order to either speed up densification or to remove isolated/detached porosity by forcefully dissolving the contained gas into the material.In this section the former effect is investigated by considering the effect of the applied stress on the equilibrium concentration of vacancies on grain boundaries.This effect can be easily included in the model of [19] as shown in the modelling section.There is a roughly linear relationship between pressure and densification rate as also experimentally observed by [24].
The pressures p ∈ {0, 50, 100, 200}MPa will be considered.The scale of the pressure is based on the magnitude of the capillary pressure O(pc) ∼ O( γs r ) ≈ 100 MPa with the surface energy γs and the particle radius r.The resulting density over time is shown in Fig. 2(a), with increasing densification rate with higher pressures.In order to determine the qualitative relationship between pressure and densification, the densification rate of the individual simulations is calculated at test densities of ρ ∈ {0.8, 0.85, 0.9}.For this, the numerical densification rate, i.e. the slope in Fig. 2(a) is calculated, followed by an interpolation across its density space, which then allows the calculation of densification rates at any density within the valid density range.This yields Fig. 2(b) which shows a roughly linear relationship between densification rate and pressure, as also observed by experiments [24].Classical theory e.g.due to Coble [25] also predicts a linear relationship between external pressure and densification rate.Hence this simple addition of pressure compares well against both experiment and theory.As sintering and creep are intimately related processes, any model of solid-state sintering should also be able to approximate creep.The present model can achieve this by setting the stress to be uniaxial, which will cause vacancy generation on grain boundaries aligned with the tensile stress direction.Four tensile tensile stresses σ ∈ {−100, −150, −200, −300}MPa in the X direction are considered, with σ = −200 MPa being run longer than the others to show porosity increase.The simulations are conducted by starting from an almost dense (97%) sintered body produced by prior simulations.The evolution of normal strain in the coordinate directions as well as porosity are shown in Fig. 3: The samples lengthen (negative strain) in the X direction and shrink in the Y and Z directions, with the Y and Z strains effectively overlapping.After an initial transient, the strain is observed to be a linear function of time as expected of secondary creep.The steady strain rate is observed to be linearly dependent on the applied stress, as is expected of Coble creep [26].Porosity is reduced at first, but with sufficient time, the porosity rises and would eventually lead to failure.This suggests that the model extension successfully models the stress dependence of Coble creep.

Coupling grain growth and densification
Densification and grain growth are intimately related and thus the concurrent simulation of both is of paramount importance for quantitative simulation of sintering.In this section, these coupled processes are investigated by parameter variations of the relevant mobilities.The mobility of an interface is controlled with the inverse mobility τ , for which high values imply low mobility.For the grain-vapor interface a small value of τgv is chosen such that the evolution is diffusion-controlled.For the graingrain interfaces, two different values of the inverse grain mobility τ gb are chosen, with their values also serving as simulation labels: Kinetically suppressed grain growth is achieved with τ gb = 100τgv and unsuppressed grain growth with τ gb = τgv.
For each of these mobility variations, different diffusion coefficients will be tested.The base case is that of Table 1, with variations on the surface diffusivity Ds and the GB diffusivity D gb by simple factors, with the factors being chosen arbitrarily.The base, τ gb = 100τ gv base, τ gb = τ gv D gb = 0.1D gb,0 , full set of investigated diffusivity combinations is {(D gb,0 , Ds,0), (D gb,0 , 0.1Ds,0), (D gb,0 , 0.01Ds,0), (0.1D gb,0 , Ds,0), (0.1D gb,0 , 0.01Ds,0)} with the color scheme being used consistently to identify the simulations.The base case (D gb,0 , Ds,0) is identified explicitly, with departures from it being used as labels with the appropriate color.The surface to GB diffusion ratio in the base case of is about 3, with this being varied from about 1/30 to 30 with the employed factors.
A qualitative, visual comparison of some of these results is given in Fig. 4 showing the GB network (blue) as well as isolated (grey) and detached porosity (red) for selected densities.It can easily be seen that in the base case without grain growth (left), no pores are detached.For otherwise the same parameters (middle), allowing grain growth leads to a significant fraction (≈ 80%) of detached porosity.Finally, if the GB diffusivity is reduced (right), then densification slows down significantly and within the allotted simulation time the maximum reached density is about 82%, whereas both simulations with quick GB diffusion achieved at least 99% density.As an exemplary quantification, the time to reach 80% density is increased by a factor of 33 even though the GB diffusion was only reduced by a factor of 10.From these images it should also be clear that even at the same density, microstructures need not be comparable.Rather, these are the complex product of the temporal interplay of densification, pore destabilization and grain growth, with a single variable being unable to capture the full picture.Videos of the time evolution, also showing a kind of Plateau-Rayleigh instability of isolated pore channels, are deposited with the supplementary material (https://doi.org/10.5281/zenodo.8263532).
Let us start the quantitative investigation by considering the evolution of surface area with density, as the reduction of surface area is the driving force for densification.It is often observed [27][28][29] that surface area and density are linearly related during sintering.This relationship is demonstrated in Fig. 5 for simulations with and without suppressed grain growth separately.The legend employed will be the same for all following plots, with colors indicating variations in diffusivity and symbols indicating suppressed grain growth (diamond) or unsuppressed grain growth (circle).The main point however is that regardless of the employed parameters, a linear relationship is observed.The linear fits' (dashed lines) coefficient of determination R 2 > 0.98 also shows this in a quantitative manner.Therefore, the model's densification behaviour is in qualitative accordance with experiments.Local deviations, such as for the blue circles, are mostly due to either a change of sintering stage or the onset of grain growth.Hence next the density and grain size evolution will be considered.
The simulation measurements in terms of density ρ, grain size G and GB area A gb are collected in Fig. 6.Focus first on the diamonds in Fig. 6(a): These represent the density of simulations with suppressed grain growth.The base case (teal), as well as the simulations in which only surface diffusion Ds is reduced (purple, orange), have essentially the same densification behaviour.This is due to vacancy absorption happening fast enough that surface diffusion contributes relatively little to neck growth compared to vacancy absorption.There is a small effect of reduced densification rate with reduced surface diffusion up to ≈ 90 % density, after which the densification rate with reduced surface diffusion is larger.This density-dependent influence of the surface diffusion is likely due to the inversion of the surface mass flux observed by Luo et al. [30].In contrast, decreasing the GB diffusion D gb (green) significantly reduces the densification rate throughout the process.In this case a significant amount of neck growth is due to surface diffusion filling the neck without concomitant vacancy absorption, which also causes the neck curvature to be lower at the same density, shown in the supplementary material.Hence decreasing the surface diffusion (blue) for this case does significantly speed up densification (≈ factor 3 less time to 80% density), though not to the original levels.The two simulations with reduced grain boundary diffusion (green, blue) will also henceforth be called slowly densifying, in contrast to the quick densification exhibited by the remaining simulations.Focus now on the effect of grain growth on densification by comparing the circular symbols (unsuppressed grain growth) and the diamonds (suppressed grain growth) in Fig. 6(a): For the quickly densifying simulations there is hardly any influence on the density evolution.This is due to densification happening so quickly that grain growth only really starts past about 90% relative density.Grain growth however does lead to pores detaching from the GBs, which at the end of the simulations causes about 1% porosity to remain in a detached state.As before, decreasing the GB diffusion causes a significant slowdown of densification, with grain growth further limiting the achievable density for a fixed time.Since grain growth has a significant effect on these, it will be discussed next.
The grain size evolution is depicted in Fig. 6(b).Comparing the base case (teal) with only reduced surface diffusion (purple, orange) shows a slight reduction of grain  growth via pore drag, since pore motion is limited by surface diffusion.The effect is not particularly pronounced due to the quick densification removing pores quickly as well.Pore drag becomes much more evident when comparing the quickly densifying simulations to the slowly densifying ones: As densification proceeds more slowly, more pores will be present at the same time and hence grain growth is slowed down.Unexpectedly, when in addition to reducing the GB diffusion the surface diffusion is reduced as well, grain growth starts much later in term of time and density.This is likely due to the sintering neck being formed only slowly in this case (D gb = 0.1D gb,0 Ds = 0.01Ds,0, blue), which also suppresses grain growth: The flux of atoms leading to grain growth is the contact area times the flux density.The flux density is due to curvature differences between grains and will be roughly comparable between simulations prior to the start of grain growth, as the initial conditions are the same.The contact area evolution differs significantly however, as is shown in Fig. 6(c) by identifying it with the GB area A gb .The simulation D gb = 0.1D gb,0 Ds = 0.01Ds,0 (blue) indeed shows much slower GB area growth.For this particular simulation, once a few grains have achieved a sufficient neck size, these grow rapidly (up to 6 times the mean grain size) until they are slowed down by porosity again, leading to stagnant grain growth for a short interval.This heavier porosity loading on grain boundaries is shown via the insets in Fig. 6(b), which show the GB network for similar grain sizes: Both simulations with reduced GB diffusivity show increased porosity on grain boundaries, with their grain growth hence being more significantly affected by pore drag.
In order to quantify grain growth with respect to experiments, the grain growth law is evaluated by fitting the grain size G data for G > 13 nm to power laws of the form G = At n , with the plots showing the results being deposited with the supplementary material.The grain size filter is employed to fit only the regime where grain growth is taking place.Exponents ranging from about 1  3 to 4 5 are observed; the experimentally observed range is 1  4 to 1 2 [31].The simulations with D gb = 0.1D gb,0 show values close to 1  3 which is often experimentally observed and can be due to pore drag.The simulations showing n > 1  2 are likely due to the start of grain growth significantly deviating from the rest of the curve.For example, if the data is filtered to G > 20 nm, then these exponents move appreciably closer to 1  2 , whereas D gb = 0.1D gb,0 achieves n = 1  3 .Hence filters for larger grain sizes will cut off more of the initial regime, but also include more data in the regime of low grain numbers.Ideally a larger system containing more grains and hence capable of reaching steady state would have been simulated.But since the primary goal was investigating densification, a too small number of grains was chosen for the initial conditions.
Based on these observations, one may conclude that below a certain ratio of Ds D gb , densification will proceed unhampered by surface diffusion.The coarse spacing of the ratios of the present investigation does not allow a particularly accurate estimate of the critical ratio, but it should be in the interval [1,10], i.e.GB diffusion should at worst be only a magnitude slower than surface diffusion.Otherwise surface diffusion will inevitably account for major parts of neck growth and hence reduce the achievable density.Given that surface diffusion and GB diffusion are similar in their grain size dependence following Herring's scaling law, this should also extend to larger particle sizes and when grain growth occurs.At some point however volume diffusion becomes, in terms of total mass transported, relevant w.r.t. the interfacial fluxes as it is less affected by grain size.
Let us close the classical consideration by verifying the experimentally observed linear relationship between grain size and the inverse square root of porosity G ∝ 1 √ P [29] with P = 1 − ρ.This relationship is plotted for the simulations which exhibited grain growth (G > 13 nm at any point) in Fig. 7 and shows a quite linear character once grain growth has started.This is quantified with linear fits, whose coefficient of determination is generally R 2 > 0.97, except for D gb = 0.1D gb,0 Ds = 0.01Ds,0 (R 2 ≈ 0.75) with its stagnating grain growth.The scatter at high densities ρ = 98% ∼ 1 √ P ≈ 7 is due to the density measurement being affected by surface roughness of the green body, which is effectively counted as open porosity.However, at this stage there are no pore channels left on the green body's surface, hence any "porosity" counted near the surface is a measurement error.Furthermore, by normalizing experimental data available in the literature with the initial grain size, it is possible to compare the present results directly with experimental results.This is shown in Fig. 7(b) and a quite good agreement between the experiments for copper and the simulation with both reduced grain boundary diffusion and mobile GBs is observed.The remaining experimental data points also show that, in terms of the sintering trajectory, different material systems can be approximated by the present model.Besides the classical quantities of density and grain size, the microstructure can also be characterized via its shape.An integral approach to this characterization is the Euler characteristic χ of the pore space, which can be used to distinguish stages of sintering [27,33]: The initial stage is characterized by a negative and decreasing value, with the intermediate stage continuously increasing the value up to a positive maximum value.In the final stage, the characteristic decreases towards zero for infinite time, representing a fully dense structure.The Euler characteristic is plotted over density in Fig. 8 for the present simulation set and good accordance to the expectations is observed.Note that grain growth has little effect on the evolution of the characteristic, since the symbols generally overlap.It is however influenced by GB diffusion via densification speed and by surface diffusion because it determines the time scale of pore evolution.In order to verify the suggestion of [33] that the maximum of χ indicates the final sintering stage, the bottom part also shows open and closed porosity separately over density.Assuming that the start of the final stage is given by the equality of isolated and open porosity, the suggestion is confirmed for the present simulations.Furthermore, it can be observed that isolated porosity is formed more slowly with lower surface diffusion.This is simply due to the pore instability being mostly limited by surface transport.Fast densification also has an apparent pore instability suppressing effect, since Ds = 0.01Ds,0 (orange) has generally less isolated porosity than D gb = 0.1D gb,0 Ds = 0.01Ds,0 (blue) at the same density.The likely origin of this is that densification progresses rapidly enough that the time to destabilization isn't reached before the pores are eliminated.
An approach which allows a more granular shape description is the calculation of shape factors characterizing individual objects.The factors employed are those of MacSleyne et al. [34] employing second moments of the mass distribution of the object to characterize its shape.These have three invariants Ωi, i ∈ {1, 2, 3} , which are characteristic for a fixed shape e.g. a sphere, a tetrahedron or the truncated octahedron (tetrakaidecahedron) suggested by Coble [35] for modelling the intermediate and final stage grain shape.The invariants can be normalized to a specific shape s, Ωi = , and thus the distance from this shape is related to how far from the value of 1 the normalized invariant Ωi is.The sphere is chosen as a reference shape and the normalized variants are collected in a vector Ω.In order to summarize the invariants, we introduce the lumped invariant employing the Euclidean norm ΩE = ||Ω|| √ 3 of the vector Ω.The factor of √ 3 stems from the use of the Euclidean norm, such that ΩE = 1 still represents a sphere; the value itself can be interpreted in the same way as the components Ωi.The invariants are separately calculated for the grains, isolated pores and detached pores, with the detached pores forming a subset of the isolated pores.Hence if only detached pores remain, their invariants will be the same as for the isolated pores.
The mean of ΩE is plotted over density in Fig. 9 for the separate structures.The standard deviation was also calculated, but in general it is quite large and hence would obscure the evolution.It was ensured that there are at least 5 objects of which the mean is calculated, as otherwise the initial trajectory is dominated by small, short-lived pores for the pores.
Focus first on the evolution of the isolated and detached pores: At high densities, reached by simulations with high GB diffusivity, both tend to ΩE = 1, which is consistent with curvature minimization.The transitory period shows a non-monotonic behaviour, likely attributable to isolated networks of pores first forming, then splitting into single pores via Plateau-Rayleigh-like instabilities.These instabilities are likely induced via grain boundaries since this is typically faster than instability growth from perturbations.[36].Since the timescale for this instability is dependent on the surface diffusion, the convergence to a spherical shape is also slower for reduced surface diffusion.High grain mobility allows pores to detach and afterwards spheroidize within the grains.This causes the invariant for the simulations with high grain mobility (circles) to generally be closer to that of the sphere than for those without, since pores on grain boundaries deform to reach the proper dihedral angle.In the slowly densifying cases (low GB diffusivity), the isolated pore networks still remain and tend to dominate the shape factor.Note that an entire pore network can become detached from the GB network, which significantly slows down its destabilization to isolated spheres.Hence apparently isolated porosity in 2D micrographs might very well be the same pore.
Focus next on the evolution of the grain structures: From the initial somewhat spherical shape, a mostly monotonic decrease is observed, with the speed being related to both interfacial diffusivities and the grain mobility.The dashed black line indicates the cube's invariant, with the solid black line indicating the truncated octahedron's invariant.All simulations pass the truncated octahedron's invariant, but do not converge towards it.One might argue that pore-laden structures have a different invariant compared to the pore-free structure.But even for the base case without grain growth, in which only few pores remain at the end, a value significantly different, i.e. ≈ 1 3 of the distance between a sphere and a cube, from that of the truncated octahedron is observed.Experimental observation of the annealing of dense (≈ 99%) strontium titanate [37] also indicate that real grains do not show the invariants of the truncated octahedron, though these of course also have effects from crystalline anisotropy.Together with the evolution of grain coordination number before the final sintering stage [19,38], this suggests that the assumption of truncated octahedra being representative of the grain structure during intermediate stage sintering should be dropped.However, the present results do not allow a constructive suggestion for a replacement shape.The generation of shapes from invariants, while possible, usually requires invariants of higher order moments as shown in [39]; furthermore, these shapes are likely not as analytically tractable as the truncated octahedron.Hence while the truncated octahedron is not a satisfactory approximation for the grain shape during intermediate stage sintering, it can still be used in geometric modelling while being aware that prefactors derived from it will likely be wrong.As a final note, the invariants for simulations with unsuppressed grain growth will  3 of the normalized invariant vector Ω is plotted over density for all simulations, with the structures of isolated and detached porosity as well as the grains being shown separately.Given sufficient densification, the porosity tends to approximate spheres, i.e.Ω E = 1.The steady invariants for the grains depend on densification rate and grain mobility, with the latter influence likely due to finite size effects.
be affected by the initial green body shape in their later stages.This is due to the few leftover grains approximating the original green body shape.Therefore, the grain invariants' values in the later stages with unsuppressed grain growth should not be taken at face value.The plots of the unlumped invariants, generally showing similar trends, are deposited with the supplementary material.

Sintering of freeze-cast structures
Freeze-casting [40] is a novel process for the production of porous materials as well as near net shape casting.A suspension of a liquid, usually water and hence assumed thus, and a target material is mixed and frozen, resulting in a microstructure of the target material and ice.The ice is sublimated next, leaving a porous structure which is then usually sintered, completing the freeze-casting process chain.A previously developed model for the phase-field simulation of the freezing part of freeze-casting [23] is employed to generate a three-dimensional freeze-cast structure.The volume fraction of target material is known for each cell of the domain and is employed to generate a sphere packing approximate the freeze-cast body.Where a sufficient volume fraction is reached, spheres of r = 8 nm are placed into a domain, approximating the continuous volume fraction field with a discrete packing of spheres.The freeze-casting simulation itself employed particles of radius 250 nm, with the change in sphere radius being for computational efficiency, as smaller particles will sinter faster; this does mean however that a smaller version of the actual structure is sintered next.The space of low volume fraction is left empty, producing macropores, and thus a rough approximation of the sublimation process is achieved.Finally, the resulting sphere packing is computationally sintered with the present model to completely simulate the process chain of freeze-casting for the first time ever.
The resulting microstructural evolution over time is shown in Fig. 10.On the left the entire structure is shown, with grains being rendered as a copper-like material, since the materials parameters of Table 1 approximate copper, and the dark lines within the structure represent grain boundaries and higher order junctions.As can be seen, even a complex, inhomogeneous structure can be sintered with the present model.Furthermore, some grain growth has occurred as would be expected from the low grain size.At the end of the simulation run, the particulate region is completely dense, with no isolated porosity remaining and no pore channels penetrating through the structure.On the right a fracture surface is shown together with a rough approximation of the porosity, rendered with a water-like apperance.This porosity is continuously reduced, with many grains already having eliminated their surrounding porosity at t = 0.09 ms, which limits their further contribution to densification.Note that almost no vertical motion is evident between t = 0.09 ms and t = 1.2 ms.Densification in this direction has stopped at about t = 0.09 ms and only the directions normal to it continue to densify.The surface in this direction also tends to have fewer grains, allowing these to bulge out more and thereby influence the strain measurement.
The shrinkage is characterized with the strain in each spatial dimension, which is found to be anisotropic: In the freezing direction, upwards and parallel to the macropore, a shrinkage of about 10 to 13% is observed, but in the directions normal to it, a shrinkage of 16 to 19% is observed.The strain range is due to inhomogeneous strain measurements.Experimental evidence of anisotropic shrinkage after freezecasting exists [20,21], but the experiments disagreed on which direction shrinks less.Farhangdoust et al. [21] observed that the freezing direction shrunk less, comparable to the present results.However, as Lichtner et al. [20] note, Farhangdoust et al. may not have removed the initial, isotropic structure generated by freeze-casting and neither the continuous skin formed on the outside of a freeze-cast cylinder.Lichtner et al. investigated the anisotropic shrinkage with experiments and the discrete element method, experimentally and simulatively observing less shrinkage in the plane normal   (a-c) show the full body, whereas panels (d-f) show a close-up for a fracture surface normal to the Y direction, also indicated by the blue plane at t = 0 ms, including a rough approximation of the porosity rendered with a water-like appearance.It can easily be seen that the porosity within the walls changes substantially between t = 0 ms and t = 0.09 ms, with many grains already having eliminated their surrounding porosity.
to the freezing direction.Following Olevsky [41] it was argued that the anisotropic shape of the pores induces anisotropic shrinkage.The reason for this is that the sintering potential is proportional to curvature, and if the curvature is anisotropic, so is the sintering potential.Hence directions which exhibit a larger curvature will densify faster than others with smaller curvature.It should be said however that the macropores' macroscopic curvature is low and it is not clear whether it can cause the observed magnitude of shrinkage anisotropy.Another reason stated in the discussion of [20] was that particles within the walls will sinter isotropically, as their contacts are isotropically distributed.Particles on the walls' surface have missing contacts and hence will tend to move anisotropically.This disregards however that the sintering potential can also be distributed inhomogeneously and anisotropically.
The present results can be explained by analyzing the anisotropy of the sintering potential, represented by how much GBs deviate from their equilibrium concentration within the simulation.Define the anisotropy factor of a property Π as f (Π) = Π Y Z Π X , i.e. simply the ratio of the property in the Y Z (average over both Y and Z directions) plane, normal to the freezing direction, to the property in the freezing direction X.The anisotropy factor of both the strain and sintering potential are plotted against each other in Fig. 11.Two regimes are evident here: Initially, there is a slight anisotropy in the sintering potential, which is sufficient to cause the strain to become increasingly anisotropic.Although the anisotropy is small (≈ 1.05), this small anisotropy causes about half of the observed anisotropy in the strain.This is due to the high driving forces for densification at this early stage.The anisotropy in the sintering potential increases eventually, due to the sintering potential being connected with the chemical potential (stress) on a particle's surface within the phase-field model.As the structure densifies, more and more grains completely lose their connection to the surface and become embedded in a surrounding grain matrix; this can easily be seen on the right of Fig. 10.Once this happens, the model assumes that there is no additional stress due to curvature and correspondingly reduces the sintering potential of the associated grain boundaries.While this occurs isotropically, the average sintering potential in a direction is dependent on both such embedded grains and surface grains.A simple model of this effect could read as e. the average sintering potential σ d in a direction d is the addition of both a volume contribution Nvσiso, acting isotropically, and a surface contribution N d s σ d s , potentially inducing anisotropy.If Nv ≫ N d s , isotropic behavior is observed, which was the case in the previous sections; though if the anisotropic contribution is the same in each direction, this would also result in net isotropic behavior.In the present geometry however, Nv ≈ 2N d s in the Y Z plane since the walls are on average about 6 particles thin.In contrast to this, Nv ≈ 5N d s in the X direction (12 particles) and hence this direction will be affected less by the surface grains, especially compared to the Y and Z directions.Thus, the observed anisotropy of the sintering potential might be due to the phase-field model giving an unwarranted extra weight to particles with exposed surfaces by virtue of their increased sintering potential σ d s .Finally, the increase in the anisotropy of the sintering potential at roughly constant strain anisotropy is due to little to no porosity remaining.This correspondingly magnifies the aforementioned effects, though given that little to no densification is possible anymore, no further increase in strain anisotropy is observed.The non-monotonic behaviour in this final region is mostly due to measurement errors via grain growth, since densification has effectively stopped.
Three other confounding factors exist: First, the present geometry only contains a single macropore open to the surrounding vacuum, raising questions of represen- tativeness.Second, within the present phase-field model it was observed that linear chains of particles only achieve particle count independent strain evolution starting from about 16 particles in the chain [19].Up to that point, increasing the number of particles in the chain decreased the strain rate monotonically.The strain rate ratio between a four and sixteen particle chain, roughly comparable to how thick the freezecast structures is in the Y Z and X directions, is on average about 1.17.This is in rough agreement with the initial increase of strain anisotropy up to about a factor of ≈ 1.2, after which the increase in sintering potential anisotropy becomes the dominant influence on the strain anisotropy.The change in strain rate between chains of different particle counts was explained via the differing sintering potentials of the chain's end particles.Since their shape evolution is only restricted by one grain boundary, they will generally have an average surface chemical potential which differs from the inner particles.Since there are many thin sections in the freeze-cast structure, this effect is present in freeze-casting sintering simulation as well, and is the likely origin of the sintering potential anisotropy.The final confounding factor is that some grain growth (≈ 10% average grain size increase) occurs.While for the present results it does not result in desintering [42], the green body resulting from freeze-casting is generally liable to experience this phenomenon -the bridges spanning the macropores are similar to the bridges spanning cracks in the classic work of Sudre and Lange [42].Since desintering would generally reduce the connectivity normal to the freezing direction, it would support Lichtner et al's [20] experimental results.In total, while the present results are quite encouraging for relatively thick structures, further research and modelling work should be done on the calculation of the sintering potential within the model.

Conclusion
It could be shown that pressure-assisted sintering and Coble creep can be effectively modelled by the presented phase-field model of sintering.This was shown by comparing the dependence of densification and strain rate on the applied stress magnitude, with the same scaling being found as predicted by theory [25,26] and found in experiments [24].Furthermore, the coupled phenomena of grain growth and densification were investigated by systematic parameter variations.A high grain boundary diffusion and a low surface diffusion are generally found to increase the densification rate.As expected, once grain growth occurs formerly isolated porosity can detach from grain boundaries and stays stable within the grains.This causes an eventual stagnation of densification prior to full density.The extent of the difference between this stagnated density and full density depends on the densification speed relative to the speed of grain growth: If grain growth occurs similarly fast as densification, 10 to 20 percent of porosity is observed to remain in the green body, whereas if densification occurs much faster than grain growth, only only about one percent of porosity remains.These simulations also allow comparing sintering trajectories with experimental data, showing good agreement for some experimental systems even though the parameters were chosen without prior reference to the experimental data.Furthermore, the grain morphology during the sintering process was analyzed and found to generally not conform to the classically assumed shape of a truncated octahedron.Finally, the sintering process was applied to a more complicated geometry than a cube, namely one which was generated by the freeze-casting process.Anisotropic shrinkage was found, with a comparison against literature data on anisotropic shrinkage of freeze-cast structures suggesting further possible model improvements.
In total, the results show the versatility and applicability of the phase-field method to the sintering process.Given accurate information about diffusivities and mobilities, it can be employed for quantitative predictions, since all the necessary qualitative features of the sintering process are present.Especially for novel high heating rate sintering processes, the control of the temperature-time curve is critical to avoid grain growth following densification, with experiments being costly and time intensive.Simulations could partially replace these and allow the determination of desirable temperature profiles.However, finding a consistent data basis in terms of diffusivities and mobilities for a single material is still challenging, as these are not only dependent on temperature but also on the grain boundary and surface characters, with impurities as well as dopants also playing significant roles.
Future modelling work will investigate variations in how the sintering potential is calculated.Separate from this, a non-conservative evolution of vacancies coupled with other species will be considered.This will enable including gas pressure, and hence pore stabilization, via another conserved species.The rigid-body assumption for individual grains can be relaxed by solving a momentum transport equation while accounting for local strain rates induced by e.g.vacancy absorption.This would allow both volume sinks of vacancies, e.g.dislocations, as well as using an orientation field for resolving grains of different orientation, which in turn could allow the inclusion of atomistically informed interface properties.

Methods
Phase-field model The model employed in this paper is based on the recent works [11,17,19].In the supplementary material a derivation of the model is given together with the necessary assumptions.The final results are summarized here.A coupled system of partial differential equations is solved, with the N values ϕα representing the surrounding vapor ϕ0 = ϕV as well as copper grains of arbitrary orientation ϕα, α > 0. The evolution of the copper concentration c is indirectly calculated with the chemical potential µ, employing the grand potential formalism [43] to decouple bulk free energies from interfacial terms.
The concentration 1 − c can be interpreted to be the vacancy concentration, which is low in the solid and high in the surrounding vapor.The meanings of the remaining terms are described in detail in [11,17].On all fields no-flux conditions are employed as boundary conditions.The equations include advective terms which are calculated based on a model inspired by molecular dynamics simulations: [19] ∆⃗ Cu = ∆u (6) in which the instantaneous velocity follows from the instantaneous displacement.The displacement jump ∆u across a GB of area A αβ is based on how many vacancies ∝ GB ∆c gb αβ dV were absorbed on a single αβ GB; this is scaled by the atomic volume Ω.The expression c − c gb eq (S) can also be thought of as the sintering potential; once it vanishes, densification via advection stops.The factor Na Vm with Avogadro's constant Na and the molar volume Vm accounts for the conversion between units of concentration c (mole fraction) and number density n.The absorption rate is based on a relaxation to an equilibrium vacancy concentration c gb eq (S) on a GB.This depends on the simulation state S in order to account for the Gibbs-Thomson effect.It is approximated via the concentration of a bulk α grain cα at its average surface chemical potential µs [17].Finally, the displacement jumps are connected to individual grain displacements u in a matrix equation [19].For efficient parallelization, a linear ansatz function for the displacements u is employed.This restricts the simulation geometry to those in which the vacancy absorption rate eq.( 5) does not have an explicit dependence on position, but can still vary from GB to GB.However, the model itself is capable of handling arbitrary variations in grain size or vacancy absorption rates as shown in the supplementary material.
By applying a bias in µs the effect of pressure on the equilibrium vacancy concentration on grain boundaries c gb eq (S) can be modelled.Furthermore, the effect of uniaxial stress can be approximated by a projection onto the local grain boundary plane ⃗ n αβ : c gb eq (S) = cα(µs + σ(⃗ n αβ )) ( 7) with the unit vector ⃗ vσ describing the direction of stress application.This effectively presumes that grain boundaries with normals aligned to the stress direction experience the full stress, whereas those perpendicular to it experience no stress.This modifies the vacancy absorption rate per grain boundary and depending on the sign of σ can also cause vacancy generation.In the case of an isotropic pressure, this projection can be skipped and the pressure p = σ can be added directly.This approach of adding stress is inspired by Coble's [25] classical model in which the capillary pressure and external pressure are simply added together.A more general approach would include the stress and strain fields in the energy functional (cf. the derivation within the supplementary material), coupling them with the chemical potential and thus also change the local equilibria.Furthermore, if the stress fields are locally resolved, the effects of stress gradients on densification would also be included.However, this requires solving for the stress evolution as well, which adds considerable additional computational effort, especially if a realistic material contrast between solid and vapor is considered.The effect of stress is not limited to densification, but also influences grain growth: Pressure modifies the diffusivity [44] and therefore the grain mobility, with higher pressures typically decreasing the diffusivity.However, at the same time, stress concentrations can lead to higher driving forces for grain growth.
Let us contrast this to the work of Dzepina et al. [13]: The model of Wang [8] is extended by adding an approximate elastic energy term into the functional.Its energy contribution is based on normal forces homogeneously acting on GBs, with the origin of the force being a deviation from the equilibrium vacancy density on the GB.This construction allows the addition of an external force (pressure) by shifting the force generated from the vacancy disequilibrium by the external force.This is similar to the present approach, except that by working directly in pressure (chemical potential) space there is no need to assume a spherial contact within the present model.However, since the motion of a particle is the result of a direct-neighbours resultant force calculation, not a global one, the densification obtained by [13] is unlikely to be homogeneous.The origin of this effect is described in the supplementary material's model derivation.Note that the stress-field is not explicitly solved for in [13].
The parameters of the simulations are given in Table 1, with more details on the choices being given in [17,19].The diffusion and interfacial energy values approximate copper at 700 K, with the choice of inverse mobility τ being such that diffusion-controlled growth is achieved for all simulations.

Data evaluation methods
The density of a green body is calculated as described in [17] with a Monte Carlo approach of sampling inside the convex hull of the green body.The strain is evaluated by dividing the domain into one-dimensional slabs per evaluated dimension d.The slabs are positioned based on the current minimum and maximum position of the particles in that dimension, with a fixed number of slabs employed per simulation.The grain-wise displacement u d,α (0) − u d,α (t) of each grain α is calculated and averaged over for each grain within the slab.With the average displacement per slab and the slab position, finite differences yield the normal strain ϵ dd per slab position, which can be averaged over if homogeneous strain per dimension is assumed.This definition will yield positive strains for a shrinking compact, consistent with the usual definition of strain in sintering.
The grain volume Vα and the equivalent grain size Gα = ( 3Vα 4π ) 1/3 easily follow from volume integrals over the phase-field domain.Integral interfacial areas (surface s, grain boundary gb) are defined by pairwise clamped sums over the appropriate phase ϕ gb = min([ A gb = 1 l ϕ gb dV (12) with the clamp being necessary as the sum of paired phase-field products is not constrained to be ≥1.The field described by ϕ gb is also used to visualize the grain boundary network.The sintering potential in the axis directions is calculated by a ϕ gb -weighted average of (c−c gb eq (S)) with the components in each direction being based on the local grain boundary normal vector.
The package cc3d's [45] connected component labelling function is used to identify isolated porosity by applying it to the vapour phase-field.This is followed by a sweep of the field: If any cell of an isolated porosity object contains a non-zero value of ϕ gb , then it is assumed to still be attached, but detached otherwise.After these labelling steps, the diffuse vapour field is split into an isolated and detached pore field which are used for analysis and visualization.
The central moments of second order µpqr required for the calculation of the moment invariants Ω [34] are calculated by integrating µ A pqr = f (A)(x − xA) p (y − yA) q (z − zA) r dV (13) over the domain with p+q+r = 2, for each object A with its center of mass (xA, yA, zA) and f (A) = 1 inside of the object and zero outside.The invariants follow directly as described in [34].The truncated octahedron's second moments are calculated with LattE integrale [46], yielding numerical values of the invariants Ω = (12.732, 162.10, 2063.8).
The Euler characteristic χ is evaluated on a triangle mesh of the vapour phase-field using its definition with the number of vertices V , edges E and faces F of the triangle mesh.Paraview [47] is used for three-dimensional visualization and plots are drawn by employing matplotlib [48].

8 ,Figure 2 :
Figure2: The influence of pressure on the density evolution is depicted in (a), with its influence on densification rate depicted in (b).With rising pressure, quicker densification is achieved.The dashed lines in (b) indicate linear fits with coefficient of determination R 2 given in the legend.There is a roughly linear relationship between pressure and densification rate as also experimentally observed by[24].

Figure 3 :
Figure 3: Strain and porosity evolution for several creep simulations, with tensile stress being applied in the X direction.After an initial transient, a constant strain rate is observed, characteristic of secondary creep.

Figure 4 :
Figure 4: Comparison of simulations with suppressed grain growth (left, ac), with grain growth (middle, d-f) and a simulation with grain growth and reduced GB diffusivity (right, g-i).The GB network is shown in transparent blue, isolated porosity in grey and detached porosity in red.The depicted states were chosen based on their density and can represent different times.

Figure 5 :
Figure 5:  The relationship between surface area and density is observed to be roughly linear, with the dashed lines indicating best fit linear functions.

Figure 6 :
Figure 6: The influence of grain growth and variable interfacial diffusivity on densification (a), grain size (b) and total GB area A gb (c).The GB networks (blue) for three simulations are shown as insets in (b) for similar grain sizes, with white space within the image indicating porosity.

Figure 7 :Figure 8 :
Figure7: The inverse square root of porosity P is plotted against the grain size G.Once grain growth starts, a rather linear relationship is observed; the dashed lines indicate fits to simulation data with G > 13 nm to account for this.In (b), only the fits together with normalized experimental data is depicted.The data are taken from the overview article[29] figs.8 (Cu) and 9 (Fe, ZrO 2 , steel, Al 2 O 3 ) and[32] (Y 2 O 3 ).The grain size is normalized to the grain size of the first data point.

Figure 9 :
Figure 9: The lumped invariant Ω E = ||Ω|| √3 of the normalized invariant vector Ω is plotted over density for all simulations, with the structures of isolated and detached porosity as well as the grains being shown separately.Given sufficient densification, the porosity tends to approximate spheres, i.e.Ω E = 1.The steady invariants for the grains depend on densification rate and grain mobility, with the latter influence likely due to finite size effects.

Figure 10 :
Figure10: A freeze-cast structure with a macropore undergoing sintering.The freezing direction was parallel to the X direction shown in the figures.The structure itself shrinks while keeping the macropore intact.Panels (a-c) show the full body, whereas panels (d-f) show a close-up for a fracture surface normal to the Y direction, also indicated by the blue plane at t = 0 ms, including a rough approximation of the porosity rendered with a water-like appearance.It can easily be seen that the porosity within the walls changes substantially between t = 0 ms and t = 0.09 ms, with many grains already having eliminated their surrounding porosity.

Figure 11 :
Figure 11: Anisotropy factor f (Π) = Π Y ZΠ X for the properties Π of strain and sintering potential.