## Abstract

We conduct experiments to investigate the sintering of high-viscosity liquid droplets. Free-standing cylinders of spherical glass beads are heated above their glass transition temperature, causing them to densify under surface tension. We determine the evolving volume of the bead pack at high spatial and temporal resolution. We use these data to test a range of existing models. We extend the models to account for the time-dependent droplet viscosity that results from non-isothermal conditions, and to account for non-zero final porosity. We also present a method to account for the initial distribution of radii of the pores interstitial to the liquid spheres, which allows the models to be used with no fitting parameters. We find a good agreement between the models and the data for times less than the capillary relaxation timescale. For longer times, we find an increasing discrepancy between the data and the model as the Darcy outgassing time-scale approaches the sintering timescale. We conclude that the decreasing permeability of the sintering system inhibits late-stage densification. Finally, we determine the residual, trapped gas volume fraction at equilibrium using X-ray computed tomography and compare this with theoretical values for the critical gas volume fraction in systems of overlapping spheres.

## 1. Introduction

The sintering of high-viscosity droplets to form a denser, connected mass is important in a range of industrial and natural scenarios, including the fabrication of ceramics [1], metals and glass, the welding of volcanic ash [2] and the vitrification of Iron Age fortification walls [3,4]. In each case, the dynamics may differ because the physical origins of the stresses that drive and oppose sintering may vary, and the materials are variably heterogeneous. We focus on what is commonly called ‘viscous sintering’—the sintering of two or more viscous droplets in the regime where interfacial tension drives fluid flow—which constitutes a viscous end-member of droplet coalescence problems. The viscous sintering problem has been studied extensively since the early theoretical works of Frenkel [5] and Mackenzie & Shuttleworth [6]. More recent studies have built on those works using both experimental [2,7–9] and theoretical constraints [10,11]. Implicit in models of surface tension-driven viscous sintering is that the liquid droplets are in the low Eötvös number and high Ohnesorge number regimes. The Eötvös number is given by
*ρ* is the liquid density, *g* is gravitational acceleration, *R* is the radius of the liquid droplet and *Γ* is the surface tension. For Eo≪1, the surface tension stress dominates the gravitational stress acting on the droplet. The Ohnesorge number is given by
*μ* is the liquid viscosity. For Oh≫1, viscosity is sufficiently high that inertial effects resulting from the surface tension-driven motion can be neglected.

In theoretical studies of viscous sintering [2,6,7,9–12], the starting geometry is usually approximated as a packing of liquid spheres (figure 1). The progress of sintering is characterized by the evolution of the gas volume fraction *ϕ* as a function of time *t*, typically non-dimensionalized by a characteristic timescale [2,9,13], where the gas volume fraction is defined as the ratio of the volume of interstitial gas to the total volume of the sample or, equivalently, the ratio of the bulk sample density to the liquid density. Previous experimental work has been limited by the resolution of measurement of sample volumes such that uncertainties on *ϕ* preclude the validation of one model over another. An example of this resolution deficit is in the poor constraint of the final gas volume fraction: the critical porosity [14,15] at which the system is in volume equilibrium over timescales much longer than the experiments.

Silicate glass is well suited to the experimental constraint of viscous sintering because, when heated to temperatures in excess of the glass transition, the resultant undercooled liquid has a Newtonian viscosity that is sufficiently high [16–18] to permit the observation of the viscous processes over timescales that are amenable to laboratory investigation. It is also the material most pertinent to ceramic applications for which the viscous sintering process has been of particular interest [10]. High temperature optical dilatometry is suitable for the study of viscous sintering of glasses because it permits quantification of sample geometry *in situ* during heating [8,19] (electronic supplementary material, figure S1). Here we provide methods to convert two-dimensional observed sample geometries into high-resolution datasets for *ϕ*(*t*) and use these to test various theoretical models, described and developed in the next section.

## 2. Theoretical framework

### (a) Viscous sintering by interfacial tension

In liquid-phase sintering, droplets that share contacts undergo time-dependent coalescence driven by the interfacial tension between the liquid and the ambient fluid in the interstitial pore space. In the high Ohnesorge number regime (equation (1.2)), this process is dominated by the viscosity of the liquid droplets, and in the low Eötvös number regime (equation (1.1)), the stress driving flow arises from the excess surface pressure *P*, which is proportional to the local radii of curvature. For spherical liquid droplets, the two principle radii of curvature are equal to one another and to the radius of the droplet (figure 1) so this excess pressure is given by the Laplace general spherical solution *P*=2*Γ*/*R*.

The characteristic timescale associated with viscous flow, driven by interfacial tension, and neglecting inertia and buoyancy, is the capillary relaxation time λ
*l* is a characteristic lengthscale. Normalizing the time of observation *t* by this timescale provides a non-dimensional capillary time

Frenkel [5] proposes a model for the growth of necks between particles that share an initial contact, in which the initial radius of the droplet *R*_{i} is the characteristic lengthscale, yielding a dimensionless neck-formation time (denoted by subscript n)
*a*_{i} is the characteristic lengthscale, yielding a dimensionless bubble relaxation time (denoted by subscript b)
*μ*(*T*) and *T*(*t*) are known. For isothermal conditions, equations (2.5) and (2.6) reduce to equations (2.2) and (2.3), respectively. These equations allow us to develop dimensionless forms of the Frenkel and the Mackenzie and Shuttleworth models in §2b,c, respectively, that are suitable for both isothermal and non-isothermal conditions.

### (b) The neck-formation model

Frenkel [5] derives a solution for the growth of the radius of a neck *R*_{n} forming between two liquid droplets of equal initial radius, as a function of time
*h* can be related to *R*_{n} and *R*_{i} by *h*≈*R*^{2}_{n}/(4*R*_{i}) (figure 1*b*). Combining this approximation with equation (2.7), Kang [21] derives a linear shrinkage equation for spheres in series, cast as the length of the system *L* relative to the initial length *L*_{i}
*et al.* [22] and Soares *et al.* [8] extend this analysis to the volumetric isotropic strain in an array of cubically packed monodisperse spheres, deriving a model for the gas volume fraction *ϕ* as a function of time
*ϕ*_{i} is the gas volume fraction at time *t*=0. Introducing the normalization

### (c) The vented bubble model

Mackenzie & Shuttleworth [6] present an idealized model of sintering in which the interstitial, gas-filled pore space surrounding the droplets is represented as an array of spherical bubbles, evenly distributed throughout the liquid. Each bubble of radius *a* sits in a spherical liquid shell of radius *S*, and shrinks under the action of the surface tension between bubble and liquid (figure 1*c*). They derive an expression for the evolution of the bulk density of the bubble–shell unit as a function of time. In this scenario, the gas is assumed to be able to escape freely (despite the lack of physical escape routes) so we term this the ‘vented bubble’ model. Conceptually, the formulation is very similar to that used in studies of the growth of bubbles in magma [23] and, in the electronic supplementary material, we provide an alternative derivation of the vented bubble model, based on a bubble-growth model.

The Mackenzie & Shuttleworth [6] solution can be cast as a rate of change of gas volume fraction to give
*N*_{b} is the bubble number density in the system. In this model *ϕ*→0 at *N*_{b} in terms of the initial gas volume fraction, which is more easily measured in practice than *N*_{b}, via the equivalence *N*_{b} is constant throughout the sintering process. As with the Frenkel [5] model, we can normalize *ϕ* by *ϕ*_{i} and use equation (2.3) to obtain a dimensionless form of equation (2.12)

The differential equations above cannot be recast to give gas volume fraction as a simple function of time or temperature. However, if we make the simplifying assumption that *ϕ*≪1, so that equation (2.13) becomes

### (d) The exponential approximation

Chiang *et al.* [24] make the observation that the relationship *N*_{b}4*πa*^{3}/3=*ϕ*/(1−*ϕ*) allows the Mackenzie & Shuttleworth [6] model (equation (2.11)) to be simplified to give
*a* and *ϕ*) rather than the initial radius and gas volume fraction (*a*_{i} and *ϕ*_{i}) that yield equation (2.12). Consequently, non-dimensionalization requires the use of a modified bubble relaxation time, in which the bubble radius is a function of time. Following the approach outlined in §2a, we couch

where the superscript asterisk indicates that *a* is a function of time. Note that this form also accounts for non-isothermal conditions if *μ* is also treated as a function of time. The dimensionless form of equation (2.16) is then
*a*(*t*) is not known *a priori*. Nonetheless, equation (2.19) and variations on a non-isothermal formulation have been widely used to describe sintering of liquid droplets [2,7,8,19,25], but without acknowledgement of the implicit approximation that *a*=*a*_{i}, for all *t*. We explicitly adopt this approximation and will assess its value against our experimental data later, in §4.

### (e) Extension to account for non-zero final gas volume fraction

It is a common observation that the gas volume fraction of a sintered mass does not reach zero, but approaches a final gas volume fraction *ϕ*_{f} [2,7,9,25–27]. This is discussed further in §2f, but we find it convenient here to accommodate this observation by substituting *ϕ*−*ϕ*_{f} in place of *ϕ* and *ϕ*_{i}−*ϕ*_{f} in place of *ϕ*_{i} in our system of equations. The normalization of *ϕ* then becomes *ϕ*_{f}, and propose that any consequent loss of fidelity will be outweighed by the advantage gained from capturing the non-zero final porosity. We leave the models to be tested against data later, in §4.

### (f) The initial bubble radius a_{i}

The radii of initially spherical glass spheres are trivial to constrain using a variety of techniques (see §3) providing constraint of the lengthscale *R*_{i} for use with the neck-formation model [5]. However, the lengthscale *a*_{i} that appears in the vented bubble model [6] and our extensions thereof in §2c,d is a less easy-to-constrain parameter. However, Torquato [28] and Torquato & Avellaneda [29] provide a rigorous expression for a mean pore size *a* occurring between particles in arbitrary packing. Their scheme can be cast for a packing of completely impenetrable ‘hard’ spheres: an arrangement identical to our initial case of packed glass beads (figure 1*a*). This is given in the form of a cumulative probability density *F*(*ζ*) of the pore size distribution for which *ζ*=*a*/*R*
*E*_{V}(*ζ*) is a pore nearest-neighbour exclusion probability function. In our system, this is a conceptual tool akin to finding the expected fraction of space available to a pore of radius *a*. To solve for *E*_{V} is a non-trivial problem that has received significant attention [28,29]. A validated expression for *E*_{V} as a function of *R* is given by Torquato [28] based on Torquato & Avellaneda [29] and reproduced here for completeness, where we cast it in terms of the gas volume fraction *ϕ*
*ζ*≥0 and contains coefficients *y*_{n} which are given by
*n*th moment of the probability density function of *ζ*, termed 〈*ζ*^{n}〉, is then related to the cumulative probability density function *F*(*ζ*) in equation (2.23) by integrating as follows:
*n*=1) value of 〈*a*〉 is *a*=〈*ζ*〉/〈*R*〉.

Equation (2.23)–(2.26) can be used to find *a* in the monodisperse limit of *R*. Torquato [28] describes the solution of Lu & Torquato [30], which further constrains a polydisperse solution which is again validated against data and reproduced here for completeness. In this form, the pore nearest-neighbour exclusion probability function is the polydisperse *e*_{V}(*ζ*) and is
*S* is the ratio of the specific surface of the polydisperse system to that of the monodisperse system at the same value of *ϕ*. *S* is given by
*R*^{n}〉 is the *n*th moment of the probability density function of the particle radius distribution. As before, the coefficients *z*_{n} are given by specific solutions, here with a dependence on *S* and 〈*R*^{n}〉
*a* as a function of *R*, the same method as the monodisperse limit is applied but where *E*_{V}(*ζ*) in equation (2.23) is replaced by *e*_{V}(*ζ*) and equation (2.26) remains unchanged.

An example of a polydisperse limit is the Schulz distribution [31] for a polydispersivity factor *m*=0. This is found by relating 〈*R*^{n}〉 to *m* by 〈*R*^{n}〉=〈*R*〉^{n}(*m*+*n*)!/[*m*!(*m*+1)^{n}] such that when *m*=0 the particle size distribution is heavily weighted to small particle sizes and with a broad tail at the high particle size limit. See §4 for application of these constraints to our system.

### (g) The concept of volume equilibrium

In the above sections, we have explored and developed non-dimensional solutions to the main sintering models for the low Eötvös, high Ohnesorge number regime in §2b–d, and presented constraint of the pore size as a function of the initial droplet size distribution in §2f. Finally, we find it useful to constrain the final gas volume fraction at which the system is in volume equilibrium. That is, the percolation threshold *ϕ*_{c} at which the gas volume fraction becomes disconnected and forms isolated bubbles suspended in the liquid. Sintering must halt at *ϕ*=*ϕ*_{c} because at this point the bubbles are no longer ‘vented’ such that the gas permeability *k*→0 and pressure equilibrium between the now-isolated bubble and the liquid is established. A period of rounding of the bubble over the capillary timescale, equation (2.1), will occur but the gas volume fraction *ϕ* should remain in equilibrium (although the bubbles may rise out of the system buoyantly over longer timescales). This *ϕ*_{c} arises from percolation theory and is the same parameter we account for empirically via §2e [14,15,28]. Both Elam *et al.* [14] and Kertész [15] constrained values of *ϕ*_{c} for systems of randomly located monodisperse overlapping spheres that are in mutual agreement regardless of the sphere size used. This is close to, albeit statistically different from, the polydisperse *ϕ*_{c} found for the two distributions of droplets investigated by Rintoul [32]. The values for *ϕ*_{c} found by Elam *et al.* [14], Kertész [15] and Rintoul [32] are collated in table 2 for comparison with our experimental data. All these studies use numerical models that randomly place spheres that are permitted to fully overlap. In our system, this can be viewed as the initially spherical liquid droplets encroaching into one another with time until *ϕ*=*ϕ*_{c}. We will compare the predicted values of *ϕ*_{c} with our data in §5b.

## 3. Experimental materials and methods

### (a) Material properties and experimental method

To assess the viscous sintering of droplets we use populations of glass spheres which, when heated above their glass transition temperature, are metastable undercooled liquids. As an experimental starting material, we use soda-lime silica glass spheres (Spheriglass^{®} A-glass microspheres product number 2530, Potters Industries Inc.). This material has been shown to have a reproducible glass transition onset *T*_{g} at a heating rate of 10 K min^{−1} of approximately 824 K and a stable mass over the temperature range 270–1670 K [2]. The temperature dependence of the Newtonian liquid viscosity above *T*_{g} can be fitted using a Vogel–Fulcher–Tammann equation of the form *μ*=*A*+*B*/(*T*−*C*) for which the fitted *A*, *B* and *C* parameters are −2.63, 4303.36 and 530.60, respectively [2], and *T* is in units of Kelvin (figure 2*a*). We note that when shear stresses acting on silicate liquids are large, non-Newtonian effects have been measured [33]; however, the shear stresses imposed by surface tension are sufficiently small that a Newtonian viscosity is sufficient to describe the rheology. Surface tension is negligibly dependent on temperature; we use the value *Γ*=0.3 N m^{−1} for dry silicate liquids [34,35].

The particle radii were measured, after sieving to size fractions below approximately 63 *μ*m, using a Beckman Coulter LS^{TM} 230 laser refraction particle size analyser with a measuring range 0.375–2000 *μ*m. The mean radius 〈*R*〉 was then calculated to be 24.7±1.6 *μ*m (figure 2*b*). Using the particle size distribution, we can calculate a predicted mean of the bubble radii interstitial to the polydisperse particles using equations (2.23) and (2.26)–(2.29), yielding 〈*a*_{i}〉 of 5.9±0.53 *μ*m (figure 2*c*).

The starting samples for our experiments were free-standing cylinders of beads. The cylinders, approximately 3 mm tall with radius approximately 1.465 mm, were formed by filling a die with glass beads, compacting with a pressure-gauged push-rod, and extracting the sample onto a ceramic plate. These samples were loaded into a Hesse Instruments EM-201 optical dilatometer, which consists of a halogen lamp, tube furnace and camera in series such that the camera is able to record a cross-sectional image of the sample (electronic supplementary material, figures S1 and S2). Samples were heated to different experimental temperatures *T*_{0} in the range 841 to 1164 K. Software developed by Hesse Instruments processes silhouettes of the sample, which are then converted to binary images collected at 1 Hz sampling rate. The cross-sectional area, height and width of the sample silhouette is computed in real time and continuously recorded in units of pixels. Samples were held at *T*_{0} until volume equilibrium was attained.

A type-S thermocouple is embedded in the ceramic sample holder within 1.5 mm of the base of the sample and is used to monitor the experimental temperature, calibrated to within 1.6 K. Low heating rates of 10 K min^{−1} were used to ensure nominal thermal equilibrium of the sample [2] during heating. Following heating, samples were held within ±2 K of *T*_{0}.

### (b) Isotropic shrinkage of a cylinder

Our samples are initially cylindrical and shrink with time. If we first assume that the geometry remains cylindrical during shrinkage, we can determine the volume as a function of time as follows. A cylinder of height *L* and radius *r* that has shrunk isotropically by a factor *α*=*L*/*L*_{i}=*r*/*r*_{i} has volume *V* =*α*^{3}*πr*^{2}_{i}*L*_{i} and maximum cross-sectional area in the *r*–*L* plane *A*=2*α*^{2}*r*_{i}*L*_{i}, where the subscript *i* indicates the initial dimensions. It follows that *V* =*V* _{i}(*A*/*A*_{i})^{3/2}, hence continuous measurement of *A* permits the cylindrical volume to be computed for all times. Volume can be converted to gas volume fraction by *ϕ*=1−*m*_{i}/(*ρV*), where *ρ* is the liquid density and *m*_{i} is the mass of the system assuming that the only contribution to *m*_{i} is from the liquid and that no mass changes occur with time. Errors associated with this measurement arise from pixel resolution and the method of pixel size calibration.

### (c) Axisymmetric volumes by the solid of rotation

The volume of the sample can also be calculated from its silhouette by treating it as a solid of rotation (electronic supplementary material, figure S2*a*), and integrating the radial distance from the axis of symmetry to the sample edge as a function of vertical position, to give *b*) and, rastering up the image in 1-pixel horizontal slices, the width of the cylinder is measured for each slice; the radius of the cylinder *r* is taken as half the width of the slice in pixels. The volume of the cylinder calculated is in voxels, which we calibrate against the final volume of the sample (see §4a). Knowing that the glass density is *ρ* and sample mass is *m*_{i}, we can convert final volume to final gas volume fraction *ϕ*_{f}. Once *ϕ*_{f} is determined, we can obtain the time-dependent gas volume fraction using *ϕ*=1−(*V* _{f}/*V*)(1−*ϕ*_{f}), where *V* _{f} and *V* can remain in voxel units.

## 4. Results and interpretation

### (a) Calibrating time-dependent porosity curves

We compare the time-dependent porosity data determined via the methods described in §3b,c (figure 3*a*). The cylindrical assumption (§3b) consistently overestimates the porosity by an average of 4.7% compared with the solid of rotation method (§3c). Throughout the analysis in this study, we use the data derived from the more general method §3c but point out here that the cylindrical assumption of §3b, which is commonly used in other studies [2,8,19], systematically overestimates the majority of the sintering process, albeit with a small error.

To calibrate the volumes computed we used X-ray computed microtomography to measure the final gas volume fraction in selected post-experimental samples that had attained volume equilibrium. Samples were mounted onto an alumina rod and clamped to the rotation rig. Images were captured using a Phoenix Nanotom E system operating at 80 kV using a 0.1 mm Cu filter to reduce beam-hardening. Images were reconstructed from 1440 projections using standard proprietary filtered back projection algorithms and the resolution (pixel sizes) range from 1.42 to 1.59 μm. Image visualization and analysis was performed using Avizo^{TM}. Pore volumes were segmented from the central region of each sample to avoid edge effects. Segmentation was performed using a standard gradient-based algorithm using the moments of the intensity distribution. All pores with volumes less than 125 voxels were discarded from the analysis as below this value the error on absolute volumes exceeds 5% [36] and these objects only comprise approximately 0.06–0.18 vol.% of the sample. Pore volumes were calculated for remaining segmented features. We note that there is a negligible average difference of approximately 0.5 vol.% in the computed pore volumes if we take the three-dimensional data from the whole sample rather than from a cropped sub-volume.

We expect that the actual value of *ϕ*_{f} measured at room temperature will be slightly lower than that measured at *T*_{0} because the bubbles may shrink in response to thermal contraction of the gas on cooling from *T*_{0} to *T*_{g}. We can estimate the magnitude of this contraction by using the ratio of the ideal equilibrium gas volumes at temperatures *T*_{0} and *T*_{g} which, assuming the number of moles of gas are constant in these isolated pores, reduces to *T*_{0}/*T*_{g}. For our experimental temperatures this equates to a maximum fractional error in *ϕ* on cooling of approximately 0.01. We do not correct the curves and rather consider this a minor source of error in our determination of the final porosity. We additionally assume that the contraction of the liquid and glass is negligible on heating or cooling [37] compared with the contraction of the gas, or equivalently, the density of the liquid and glass is taken as constant. This assumption is justified because the thermal expansion coefficient for silicate glass is approximately 10^{−5} K^{−1} [38]; hence the variation in glass density between room temperature and *T*_{g} on heating or cooling is less than 1%, which is small compared with density changes resulting from sintering.

### (b) Testing models using best-fit droplet and pore radii

After applying the solid of rotation to obtain time-dependent volumes of high temperature experimental samples *in situ* in voxels and subsequently converting these to *ϕ*(*t*) (§3), we obtain the data presented in figure 3*b*. The curves all show a rapid onset of *ϕ* decay followed by a long tail at high values of *t*. While the shape of the curves is similar across all *T*_{0}, the absolute rate of this process is systematically dependent on *T*_{0}.

We use our experimental results to test the theoretical models presented in §2. In all cases, the models are tested in dimensionless form, which necessitates transforming the raw datasets—i.e. *ϕ*(*t*) for each experimental run—into *ϕ* is trivially non-dimensionalized as *t* as *Γ* is constant at *Γ*=0.3 N m^{−1}. If equation (2.2) is used, the initial droplet radius *R*_{i} is either taken from the measured particle size distribution, or treated as a fitting parameter. If equation (2.3) is used, the initial pore radius *a*_{i} is either calculated from *R*_{i} following the approach outlined in §2f, or treated as a fitting parameter.

Where non-isothermal behaviour is important, non-dimensionalizing *t* is slightly more complex. The temperature–time data for the run is used to calculate *Γ* is constant. If equation (2.5) is used, the initial droplet radius *R*_{i} is either taken from the measured particle size distribution, or treated as a fitting parameter. If equation (2.6) is used, the initial pore radius *a*_{i} is either calculated from *R*_{i} following the approach outlined in §2f, or treated as a fitting parameter.

In §4b, we allow the initial droplet radius *R*_{i} and initial pore radius *a*_{i} to vary freely, as fitting parameters; in §4c, we constrain these radii based on measured particle size distributions. This two-step analysis allows us to assess the consistency of each model across the large range of liquid viscosities investigated before generalizing the models without any fitting procedure.

#### (i) The neck-formation model

For each dataset, the best fit of the neck-formation model (equation (2.10)) is found with *R*_{i} as a free fitting parameter using a least-squares regression procedure. When the experimental time *t* is normalized using the best-fit value of *R*_{i} obtained, the data collapse to close to a single curve. Compared with the model itself (equation (2.10)), this produces a moderate fit for all temperatures (figure 4*a*). Using equation (2.2) and the definitions of the isothermal viscosity μ and the surface tension *Γ* discussed, a least-squares regression method can be used to fit the linear relationship across all temperatures in figure 4*b*, which determines a best-fit characteristic radius *R*_{i}. For all values of μ, this yields an estimated average *R*_{i}=11.4±1.2 μm (black curve figure 4*b*) via equation (2.2). This *R*_{i} compares favourably with the mean radius from the measured particle size distribution of 〈*R*_{i}〉=24.7 μm. When a single best-fit characteristic *R*_{i} is found for all isothermal temperatures, the coefficient of determination is 0.984 (figure 4*b*; table 1).

The model of Frenkel [5] tested here is based on the formation of necks between liquid droplets and as such is expected to describe the early part of the sintering process better than the later part. This is confirmed by figure 4*a*, in which the model curve decays to

In order to explore this further, we repeat the fitting process multiple times for all datasets, each time fitting a slightly greater fraction of the data. In each case, we start the fit at *t*=0) and fit to a volume fraction *R*_{i} for that fraction of the data. Figure 4*c* plots *R*_{i} against *R*_{i} normalized by 〈*R*_{i}〉, such that a value of *R*_{i}/〈*R*_{i}〉 closer to 1 indicates a good fit between computed and measured particle radius. The plot demonstrates that the model and data are in closest agreement when we fit only the early sintering data, and that the fit worsens as more data are included in the fit. The data presented in figure 4*a*,*b* were calculated using *ϕ*_{i}, corresponds to *ϕ*=0.2, which is the value above which Prado *et al.* [7] claim the Frenkel model is applicable. However, here

#### (ii) The vented bubble model and exponential approximation

Next, we test the vented bubble model, which was modified after Mackenzie & Shuttleworth [6] in §2c. This model uses the dimensionless time *t* is normalized to the bubble capillary timescale λ_{b}. The model, which is solved numerically, can either go to *a*) or to *b*). In the former case, the best-fit timescale λ_{b} for all experiments yields a best-fit bubble radius *a*_{i} of 11.0±1.2 μm and in the latter case, a best-fit *a*_{i} of 9.6±0.9 μm, both of which compare favourably with the 〈*a*_{i}〉 value of 5.9 μm computed following the approach described in §2f (equations (2.23) and (2.26)–(2.29)). To fit for the characteristic *a*_{i}, we use a least-squares approach and equation (2.3). Finding a common *a*_{i} for all experimental temperatures yields a coefficient of determination of 0.993 (figure 5*c*; table 1). The agreement between the fit *a*_{i} and the measured 〈*a*_{i}〉 is manifest in the success of the collapse to a single curve of *a* and 6*b*). The final gas volume fraction *ϕ*_{f} refers to the minimum observed value that is confirmed by X-ray micro computed tomography (§4a) and is measured to be approximately 0.03 (table 2).

The Mackenzie & Shuttleworth [6] model yields an analytical approximation when *ϕ*_{i}≪1 by equation (2.15) and (2.21). We test this against the experimental data in figure 5*d*,*e*. Whether *d*) or to the empirically observed *ϕ*_{f} (figure 5*e*), the best-fit *a*_{i} is within error of the estimated 〈*a*_{i}〉 for all experimental values of μ (figure 5*f*). For our samples, for which the average *ϕ*_{i}=0.45±0.02, the small *ϕ* approximation provides an excellent collapse of the data to a single curve for a common characteristic radius *a*_{i}, with a coefficient of determination of 0.994 (table 1).

Finally, we test the commonly used [2,7,8,19,24] exponential approximation of the vented bubble model (equations (2.19) and (2.22)) as described in §2d, where we note the implicit assumption that bubble radius is independent of time. Despite this assumption, which must, in reality, be violated, the results of fitting for the timescale λ_{b} for both the *g*,*h*) and almost indistinguishable from those of the small *ϕ* approximation, resulting in average best-fit radii in excellent agreement with 〈*a*_{i}〉 and a coefficient of determination of 0.995 (table 1).

All best-fit radii and coefficients of determination are summarized in table 1.

### (c) Testing models without fitting

In figure 6, we show all data from figure 3 normalized by the capillary timescale λ_{b} in which the lengthscale is the radius of the bubbles interstitial to the mean of the particles 〈*a*_{i}〉 estimated using Equations (2.23) and (2.26)–(2.29), which we have shown provides a reasonable approximation to the lengthscale controlling the best-fit timescales across all experiments (figure 5). This permits all data to be collapsed to a single description of *ϕ* approximation of the vented bubble model; and (iii) the exponential approximation [24]. We show both the solutions when *a*) and when *b*). The vented bubble model modified from Mackenzie & Shuttleworth [6] and the small *ϕ* approximation thereof are almost indistinguishable from one another for the values of *ϕ*_{i} represented by our samples, and both provide a good agreement with the data. However, we note that there is a systematic deviation from the predicted behaviour at values of

## 5. Discussion

### (a) The permeability problem

Figure 6 shows that the high temperature data deviate from all models at high values of _{D} when the gas is driven through the network by the capillary pressure in the pores, such that
_{g} is the gas viscosity, *L* is the sample lengthscale (in our case the sample length itself; figure 2), and *k* is the gas permeability. We can estimate *k* for our initial system by using the universal scaling for impenetrable (hard) sphere packings proposed by Martys *et al.* [39], which was validated for glass sphere packings, sandstones and sintered materials [40]:
*s* is the specific surface area of the pores (i.e. the surface area of the pores normalized to the sample volume). This yields values of *k* for our initial system of approximately 10^{−11} m^{2} using values of 〈*a*_{i}〉 estimated following the approach in §2f. It should be noted that Martys *et al.* [39] offer a second, very similar solution which is dependent on the *n*th moments of the particle size distribution; however, this is less easy to implement as a function of time as we are concerned with the evolution of pores and not particles. Using values appropriate for our initial system (*L*=2.98 mm and μ_{g}=10^{−5.5} Pa s for argon gas), equation (5.1) predicts λ_{D}≈10^{−2.22} s.

To assess how λ_{D} evolves with time, we need to know not just *k* and *a* for the initial system, but how these parameters evolve with time. The length *L* is extracted, as a function of time, from the image analysis described in §3 and is continuously recorded by the optical dilatometer. The radius of the pores as a function of the porosity is given by
*a* can be converted to a bulk specific surface area by *s*=3*ϕ*/*a* assuming, for these illustrative purposes, a monodisperse value of *a*. Equation (5.3) allows us to convert *ϕ*(*t*) (see §3) into *a*(*t*) and *s*(*t*). Now equation (5.1) can be assessed as a function of time, or more usefully, of *a*). We can see that the Darcy timescale is a strongly nonlinear function of time and goes toward infinite values at high values of _{D}/λ_{b} which yields a Capillary Darcy number

If equation (5.4) represents the competition of timescales that is responsible for the discrepancy between our experimental data and the vented bubble model, we would expect that, as sintering progresses, the system should transition from Da<1 (i.e. sintering is unimpeded by permeable gas escape) to Da>1 (i.e. sintering is impeded by permeable gas escape) prior to attainment of volume equilibrium at *ϕ*_{c}. In figure 7*b*, we show the residuals when we subtract the measured gas volume fraction *ϕ*_{meas} from that which is predicted by the vented bubble model *ϕ*_{VBM} as a function of the calculated Da. Here a residual of zero represents perfect agreement between *ϕ*_{meas} and *ϕ*_{VBM}. At very low values of Da the residuals are approximately 0, suggesting that, when the permeability of the system is high, the model well predicts the sintering. However, at *ϕ*. At Da∼10^{−6} the calculated permeabilities of the samples fall in the range ^{−4} after which the residuals decrease as the samples approach final porosity. Final gas volume fraction is associated with a residual of approximately 0.03, and Da in the range

If decreasing sample permeability were responsible for the mismatch between data and models as the sintering process nears completion, we might expect the mismatch to be greatest for Da≈1. In fact, both the onset of the deviation between *ϕ*_{meas} and *ϕ*_{VBM}, and the peak residual, occur while Da≪1 for which, as stated, Da is calculated for the bulk sample. Another alternative explanation for the mismatch arises from consideration that *a* is not monodisperse (§2f). The consequence is that isolation of pores may occur over a range of *ϕ*_{iso}, would be an increasing function of time at *ϕ*_{iso}=*ϕ*_{c}. These propositions certainly require further investigation, which could be accomplished via *in situ* 4D experimental datasets.

We have shown that, although the exponential approximation includes the erroneous implicit assumption that *a*≠*f*(*t*), it nonetheless provides the best fit to data of all the models that we present, consistent with its widespread adoption in the literature [7,8,19]. Given that the exponential approximation always overestimates *ϕ* compared with the vented bubble model (figure 8), while the vented bubble model tends to underestimate *ϕ* compared with the data at moderate and high Da, we speculate that the good performance of the exponential model might result from a fortuitous cancelling of errors. Figure 8 shows that the poorer performance of the exponential model might be expected if *ϕ*_{i} were larger than that tested here and in other studies.

### (b) The critical gas volume fraction and volume equilibrium

The critical value of gas volume fraction at which pores become completely isolated has been the subject of extensive investigation, and the percolation threshold for monodisperse systems of fully penetrable spheres has been found to be in the range 0.033<*ϕ*_{c}<0.034 for any monodisperse system [14,15]. Investigation of the critical gas volume fraction for polydisperse fully penetrable spheres has suggested a small effect of polydispersivity, that results in a lower value of *ϕ*_{c}=0.029 [32]. These studies use stochastic methods to randomly locate spheres in volumes such that the number density of the overlapping spheres controls the bulk solid volume fraction. In our system, however, we have a constant mass of liquid (equivalent to their solid volume fraction) and it is the volume of the whole system that controls the bulk liquid fraction. We find a critical gas volume fraction of *ϕ*_{c}=0.03±0.008 for all temperatures investigated. This is in excellent agreement with these previous works [14,15]; however, the resolution of this value, manifest as the quoted uncertainty, is not sufficient to conclude whether our system of spheres is closer to the theoretical values of *ϕ*_{c} for the monodisperse [14,15] or the polydisperse [32] simulations. These results are summarized in table 2.

The *ϕ*_{c} discussed above is the equilibrium value, encapsulated by our model when we set *ϕ*_{f} to the empirically constrained *ϕ*_{c} (see §2e). However, there are processes that can modify *ϕ*_{c} from equilibrium by halting viscous flow prior to the completion of sintering. If sintering is coincident with crystallization of the liquid, or if sintering occurs during cooling through the glass transition, then sintering can stop prior to attainment of the equilibrium *ϕ*_{c}. In a crystallizing system of sintering droplets, high crystal volume fractions can result in a local yield stress, which may be larger than the surface tension stress driving flow at the droplet surfaces, or can result in a jammed state in which bulk flow is altogether inhibited [41]. However, when liquid droplets are the size of those studied here, rigid crystals tend to form at their surface rather than in the droplet interiors [42]. A model treatment of the effect of this surface phenomenon on sintering rates remains enigmatic, particularly as the viscous flow driven by surface tension is also dominantly local to the surface.

Eberstein *et al.* [26] investigated systems of glass fragments mixed with rigid crystal particles and found that the empirically observed maximum linear shrinkage of cylindrical samples (manifest as the minimum final sample height *L*_{f}) scales with a bulk crystal volume fraction *ϕ*_{x} such that *ϕ*_{x} exceeds the crystal fraction above which the crystal suspension is jammed at the surface. Therefore, sintering can halt before the crystal-free *ϕ*_{c} is attained. Beyond this observation, a model which satisfactorily scales *ϕ*_{f} with *ϕ*_{x} is lacking and would be a valuable future contribution. Parametrization of the processes that affect the viscous sintering rate, including crystallization, glass composition, non-isothermal trajectories and droplet initial geometries, have been investigated in detail due to their application to ceramic and glass fabrication [8,25,43] or to natural systems in which sintering takes place, such as volcanic interiors [2,9]. However, further constraint of the parameters affecting *ϕ*_{f} is needed (as discussed above). The concept of volume equilibrium is therefore one that requires knowledge of the time- and temperature dependence of the material properties and we suggest that attainment of the equilibrium *ϕ*_{c} is reserved for a special class of metastable undercooled liquids which do not readily crystallize.

## 6. Conclusion

We show that viscous sintering can be modelled using a modified version of the original Mackenzie & Shuttleworth [6] theory. We test this and other models and discuss their limitations for real systems of liquid droplets. We find an equilibrium final gas volume fraction preserved in the samples that is in agreement with theoretical values for overlapping penetrable spheres. Finally, we introduce a conceptual framework for the incorporation of gas permeability into a sintering Darcy timescale for gas escape, and show that decreasing gas permeability may have an important impact on late-stage sintering rates. We propose that the most pressing unresolved complexities in the sintering of high-viscosity droplets, namely the low-permeability final phase of the sintering process identified here, could be explored by future high temperature *in situ* X-ray techniques.

## Data accessibility

Processed data are available in the electronic supplementary file and raw data are available upon request from the authors.

## Authors' contributions

F.B.W. conceptualized the work and performed the experiments with the help of J.S. under the supervision of D.B.D. and B.S., and J.V., E.W.L. and F.B.W. performed data processing and analysis. All authors wrote the manuscript.

## Competing interests

All authors declare no competing interests.

## Funding

Funding was provided by the Fellowship Development Fund of the Department of Earth Sciences (Durham University), the European Union's seventh program for research under grant agreement VUELCO (282759) and the European Research Council's advanced investigator grant EVOKES (247076).

## Acknowledgements

We thank F.W. von Aulock, Y. Lavallée and L. Colli, for fruitful discussion and K.-U. Hess for facilitating the X-ray tomography.

- Received November 9, 2015.
- Accepted March 11, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.