Sintering of polydisperse viscous droplets

The full-text may be used and/or reproduced, and given to third parties in any format or medium, without prior permission or charge, for personal research or study, educational, or not-for-pro t purposes provided that: • a full bibliographic reference is made to the original source • a link is made to the metadata record in DRO • the full-text is not changed in any way The full-text must not be sold in any format or medium without the formal permission of the copyright holders. Please consult the full DRO policy for further details.


I. INTRODUCTION
Coalescence and sintering of many droplets is relevant to our understanding of early planetesimal formation [1], welding of volcanic ash in volcanic interiors [2][3][4][5] or volcanic supereruptions [6,7], to our efforts to fabricate partially sintered ceramics and glasses [8][9][10][11], and to recent advances in our ability to manufacture inexpensive building materials on other planets [12], among numerous other applications.In all of the above cases, the high temperature viscous droplets that sinter together are typically polydisperse, the temperature environments are nonisothermal, and the packing of droplets is random and heterogeneous.
The coalescence of viscous droplets has received attention for the case where two [13][14][15][16] or a few [17] droplets coalesce to form single larger droplets.Extensions of models for a few interacting droplets to the coalescence of many-droplet systems have been made [10,11,18,4,19,20].In this case the metric of coalescence becomes the interdroplet fluid fractionwhich is a gas volume fraction φ in the common case that the interstitial fluid is air.These extensions remain limited to arrays of droplets in which the extension of two-droplet kinetics to many-droplet systems relies on definition of a regular packing type for the droplets [11,18,21].For more realistic random and heterogeneous distributions of droplets, no general microphysical model for large droplet populations has been attempted.
In the regime where body and inertial forces are negligible, such that the Eötvös number Eo = ρgR 2 / is small and the Ohnesorge number Oh = μ/ √ ρ R is large, coalescence is driven only by the Laplace pressure P L = 2 /R of the interface between droplets and the interstitial fluid.Here ρ is the density of the droplet, the density of the interstitial fluid is taken to be negligible, g is the acceleration due to gravity, R is the droplet radius, is the interfacial tension, and μ is the droplet viscosity.
The Laplace pressure drives fluid flow to the point contacts between droplets such that the length of a two-droplet system decreases as the radius of the fluid collar, or neck, between the droplets grows [22].The characteristic time for this flow into the neck is the capillary time λ n , which depends on the initial droplet radius R i , according to λ n = μR i / [11,4,23] (where subscript i denotes an initial value).Simulation results agree very well with experimental data for droplet-droplet pair shortening, and show that the short-time dynamics can be scaled with this model in systems of up to ten droplets in chains or defined arrays [17]; however, this model for the coalescence of two droplets cannot easily be generalized to systems of many droplets because λ n scales with the coordination number of the droplets [17], which is difficult to predict, control, or measure.Experimental data for sintering of moderately polydisperse packs of many droplets in random heterogeneous coordination can be described well by this model up to times t ∼ 0.6λ n , but after this time the system deviates significantly [23].Therefore, the long-time dynamics of sintering of many droplets is not well captured by the neck-formation model.
An alternative class of models that account for coalescence of large systems of many droplets approximates the interstitial space between the droplets as an array of nestled bubbles which shrink under interfacial tension [24] (Fig. 1).Conceptually, these bubbles are vented to the outside of the system such that the gas pressure in the bubbles is in equilibrium with the ambient pressure and the droplet fluid pressure; we term this class of models the vented bubble model [23].For this alternative view, the initial droplet length scale R i is exchanged with the initial bubble length scale a i in both Eo and Oh, and the capillary time that is characteristic of the flow around the bubbles becomes λ b = μa i / [23].We have previously shown that, when a rigorous description is found for the radius distribution that best describes the array of bubbles among randomly packed droplets [25,26], the vented bubble model is highly effective for monodisperse populations of droplets [23].
So-called cluster models have been proposed [10,11,18,27,28], which combine the neck-formation model and the vented bubble model.These cluster models rely on a slope-matching method, in which the neck-formation model is adjusted empirically so that it transitions smoothly FIG. 1.The model geometry and experimental droplet size distribution.(a) A spherical pore or bubble of radius a in a liquid shell of radius L where the approximation of the pore spaces as bubbles yields L ≈ a + R. The gas escape vent is conceptual, and is not physically involved in the bubble shrinking.(b) A 3D approximation of a portion of a pack of droplets of different sizes; the curvature of the pore phase is approximated as a polydisperse distribution of bubbles nestled between the droplets.(c) The measured particle size distribution of the glass spheres used in the experimental work which, once heated to temperatures above the glass transition interval, are viscous liquid droplets.Inset: the initial gas volume fraction of the packed spherical droplets; note that the polydispersivity is defined as S = R 2 i R i / R 3 i , where R n i is the nth moment of the distribution of R i [25].Here all data are displayed and in subsequent figures only representative curves are given.(d-g) The 2D projections of the experimental cylinder packs of droplets on ceramic plates inside the furnace for an intermediate polydispersivity of droplets (S = 0.78) at four instances during the coalescence process; the corresponding porosity φ is calculated by integrating the edge of the sample about the axis of rotation (see Sec. III).
to the vented bubble model at a particular bulk sample density [11].Comparisons of the cluster model with the vented bubble model have shown that the latter is in fact a better description of continuous experimental data, which attenuates the universal utility of the cluster model adjustment [3].In this study we extend the vented bubble model to include polydisperse populations of droplets, which are more commonly found in experimental and natural coalescence scenarios.

A. Sintering of monodisperse viscous droplets: Vented bubble model
A first form of the vented bubble model was given by Mackenzie and Shuttleworth [24] based on a model for a single bubble shrinking in a liquid shell (Fig. 1) scaled to the system of n bubbles.The derivation of this model is available in both the original Mackenzie and Shuttleworth [24] and in a work that applied their full model to experimental data [23].By normalizing n by the volume of the liquid shells, they formulate the scaled vented bubble model in terms of the number density of bubbles N (number of bubbles per unit volume of liquid), yielding where t is time.Here we generalize and nondimensionalize the Mackenzie and Shuttleworth [24] vented bubble model for a monodisperse system of bubbles, before extending it to account for a polydisperse distribution.We note that, for a monodisperse system of bubbles, N is related to the initial gas volume fraction φ i by 4Nπa Nondimensionalizing time as tb = t/λ b and gas volume fraction as φ = φ/φ i gives Since μ in tb is dependent on temperature T for most liquids by some function μ(T ), Wadsworth et al. [23] develop a modified dimensionless time to extend the utility of Eq. ( 3) to nonisothermal conditions: This can be used when μ(T ) and T (t) are known; Eq. ( 4) reduces to the simple tb = t/λ b when the system is isothermal.We emphasize that the approximation of the real interstitial pore spaces as a population of polydisperse hypothetical bubbles has not received detailed attention.Below, we test this approximation and, on the basis of the discrepancy between the number density of the real pores and the hypothetical bubbles, we scale the model result to account for the inexactness of this approximation.

B. Polydispersivity of interdroplet pore sizes
To take account of the polydispersivity explicitly, we must define the distribution of model bubble sizes that occurs among a random pack of droplets, which may themselves have a polydisperse size distribution.The probability density function for pore sizes occurring among randomly distributed droplets is F (ζ ) where ζ = a i / R i and R i is the mean initial radius of the distribution of R i .We use the solution given by Lu and Torquato [25] in which F (ζ ) is dependent on φ and the polydispersivity of the droplets S. We additionally adopt their definition of S = R 2 i R i / R 3 i where R n i is the nth moment of the distribution of R i [25].S = 1 therefore represents the monodisperse limit.The probability density function F (ζ ) is given by where h V (ζ ) is a pore nearest-neighbor probability function.
Lu and Torquato [25] give a rigorous solution for h V (ζ ), which, when cast for our system, is where e V (ζ ) is the associated exclusion probability, which is given by for which parameters z n are related to the distribution of R: With the continuous probability density function F (ζ ), we can solve the dimensional form of the vented bubble model [Eq.( 2)] for every contribution of ζ , and therefore every contribution of a i , and integrate to find the evolution of the weighted average of the normalized gas volume fraction, which we term φ Procedurally we use Eq. ( 2) to give φ(t), which is then normalized by φ i to find φ(t) for each class of ζ before using Eq. ( 9) to find φ = (t).Equation ( 9) therefore represents a sintering model that explicitly accounts for droplet polydispersivity in random heterogeneous packing.Our method differs from that of Prado et al. [11], who present a bubble cluster model that can be used for polydisperse suspensions.That model is built on a geometrical argument for the total number of droplets of a given radius that can be accommodated around a droplet of some larger radius; by summing these droplet cluster arrangements over the droplet size distribution and weighting the results by the volume fraction, they arrive at a model for polydisperse droplets.By contrast, we use a statistically rigorous approach to account for truly random and heterogeneous distributions of droplets.

III. EXPERIMENTAL MATERIALS AND METHODS
We conduct experiments using populations of soda-lime silica glass spheres (Spheriglass A-glass microspheres from Potters Industries Inc.) for which the distribution of particle sizes is varied from run to run, and measured using a Coulter LS 230 laser refraction particle size analysis tool (resultant distributions given in Fig. 1).These spheres have a glass transition onset temperature of 824 K on heating at 10 K min −1 , a stable mass over the full experimental temperature range, and a nominally temperature-independent interfacial tension with air of = 0.3 N m −1 [29].The temperature dependence of viscosity has been measured and is well approximated by an empirical parametrization of the Vogel-Fulcher-Tammann equation where log 10 (μ) = A + B/(T − C) where A, B, and C are −2.63,4303.36, and 530.6, respectively [23,3] when T is given in Kelvin.These glass spheres were loaded into a die and a pushrod was used to compress a cylinder of particles of ∼3 mm height and ∼1.45 mm radius where the force of the pushrod was gauged at 1.5 N to be sufficient to reach a reproducible packing for a given polydispersivity (initial porosities given in Fig. 1).When pushed from the die, the cylinder is freestanding and so edge effects are removed.The cylinder was loaded into a Hesse Instruments EM-201 optical dilatometer which consists of a lamp, furnace, and chargecoupled device (CCD) camera in series, where the camera records the time evolution of the approximately rectangular silhouette of the cylindrical sample.The samples were heated linearly at 10 K min −1 to a high final temperature, at which a minimum volume was reached.A type-S thermocouple sits in the sample plate and records the temperature, calibrated to within 1.6 K of the absolute temperature by measuring the solidus of gold.The two-dimensional (2D) images are converted to three-dimensional (3D) volumes V by rotating the radius r (recorded in pixels of edge ∼10 µm) of the initially cylindrical sample around a central axis and integrating each resultant disk volume over the sample height H c by V = H c 0 πr 2 dy where y is the vertical coordinate in the cylinder.Volume is then converted to φ by measuring the final gas volume fraction φ f using x-ray computed tomography on the postexperimental sample and assuming that all volume changes arise from changes in φ (justified by a constant sample mass throughout the experiment).X-ray computed tomography was performed using a Phoenix Nanotom E with 80-kV beam, a 0.1-mm Cu filter, and capturing 1440 projections at 1.42-1.59-µmpixel resolution.5)-( 8) to compute the evolution of φ, and dashed lines refer to the model where a i (tuned monodisperse) is a fitting parameter in the solution of φ(t) or where a i is the iterated result of the tuning of the polydisperse model (tuned polydisperse; see Fig. 3).
Errors on the determination of φ arise from the pixel resolution in (1) the optical dilatometry method and (2) the x-ray computed tomography.In (1) we achieve pixels that have a maximum edge length of 10 −5 m and in (2) we achieve voxel edge lengths of 1.4 μm [23].These errors are reported in Figs.2-4.

IV. RESULTS
Experimental results for a range of droplet size distributions show that the decay of φ with time is nonlinear, with a nonzero final φ at long t, which corresponds to the percolation threshold φ c below which remnant gas bubbles are isolated and the system becomes impermeable (Fig. 2).Under   2) with Eq. ( 4) where a i is a fitted length scale].Solid line shows the polydisperse model using Eqs.( 2)-( 5), for which a i is calculated using a best-fit Ri , as described in the main text.The normalized gas volume fraction φ or φ = as a function of dimensionless time tb where the data and models are collapsed by using (a) the best-fit monodisperse a i , or (b) the best-fit polydisperse a i .Inset: the approximate equivalence between the best-fit monodisperse characteristic radius a i and the best-fit mean of the polydisperse characteristic radii a i for all droplet polydispersivity S.
isothermal conditions, φ c represents volume equilibrium; however, temperature changes after the time at which the percolation transition occurs can change V by expansion (increasing temperature) or contraction (decreasing temperature) of the isolated bubbles.We halt our experiments when φ = φ c to preserve the equilibrium final gas volume fraction, and we neglect minor changes in φ that occur on quenching to the glass transition temperature interval.This is justified by considering that, on cooling from peak temperatures of 1123-1173 K to the region of the glass transition at ∼814 K, bubbles will shrink by no more than 1% of their initial volume based on an ideal gas at 0.1 MPa [23].We take this as an additional source of error.At times greater than the time at which the system reaches φ c , the sample becomes impermeable; in the vented bubble model paradigm, this is equivalent to closure of the hypothetical vent through which the gas can escape from the bubbles.Therefore, we restrict the analysis that follows to the case where t O(λ b ), for which the vented bubble model is valid.

A. Testing the monodisperse vented bubble model
We compare the monodisperse model with experimental data in Fig. 2(a).We first compute the evolution of φ with time using the a i calculated from Eqs. ( 5)-( 8) (φ i is shown in Fig. 1 for each polydispersivity).In this case, the model [solid line in Fig. 2(a)] always predicts that the sintering occurs more quickly than is shown by the data.A much better fit is found between the data and the monodisperse vented bubble model when a i is treated as a fitting parameter [dashed line in Fig. 2(a)].We term this model, which uses a best-fit a i , the "tuned monodisperse vented bubble model."However, in detail this model gives a φ(t) curve that is too steep compared with the data-it overpredicts φ during the early stages of sintering, and underpredicts φ during the later stages.In Table I we report the values of a i and a i relative to the measured R i for each polydispersivity S.

B. Testing the polydisperse vented bubble model
As a second step, we compare the polydisperse solution of the vented bubble model with the experimental data [Fig.2(b)].The procedure for obtaining this solution is illustrated in Fig. 3, using the data for S ≈ 0.4 as a worked example.First, the droplet size distribution is measured and the mean of the distribution R i is found [Fig.3(a)].Then Eqs. ( 5)-( 8) are used to find F (ζ ) where ζ = a i / R i and the mean pore size is found, a i [Fig.3(b)].Equation ( 2) is then applied to find φ(t) for each value of a i [Fig.3(c)], where a i = ζ R i ; this, in turn, is used in Eq. ( 4).The polydisperse solution is then found by the integral technique described by Eq. ( 9) and is given as the solid line in Figs.2(b) and 3(c).Here we convert the solution to the same form as the experimental data using φ(t) = φ = (t)φ i .The resulting polydisperse model [solid line in Fig. 2(b)] provides a good approximation of the data at high polydispersivity (S ∼ 0.4) but is systematically less good at capturing the data as S → 1.As for the monodisperse case, the polydisperse model predicts that the sintering occurs more quickly than is shown by the experimental data.
As for the monodisperse case, we can greatly improve the fit of the model to the data by including an empirical shift.To do this, the procedure described above is iterated for different values of the mean particle size R i until the polydisperse solution matches the data with a maximum coefficient of determination.The best-fit characteristic particle radius is then R i [which is given superimposed on the original distribution of R i in Fig. 3(d)].The new pore size distribution and resultant mean pore radius are then given in Fig. 3(e) and the final result of the data and model are given in Fig. 3(f).This procedure is repeated for all of the datasets, and the results of this "tuned polydisperse vented bubble model" are shown as dashed lines in Fig. 2(b).The polydisperse solution matches the data well across all S up until some critical point at long times (see Sec. V).

V. ANALYSIS AND DISCUSSION
The experimental data and model fits can be normalized and nondimensionalized (Fig. 4).For the monodisperse case, this is achieved by using a i in λ b (T ) via Eq.( 4) and produces an excellent collapse across all datasets even when a i is used to normalize the tuned polydisperse model curves [Fig.4(a)].Nondimensionalizing time using the polydisperse radii in λ b (T ) is nontrivial as there is no longer a single λ b (T ) for each distribution of radii, but a range.However, we find that using the mean of the best-fit polydisperse radius distribution a i to compute tb also produces excellent agreement across all experimental data [Fig.4(b)].We show that the agreement between both nondimensionalization types (using a i or a i ) is a result of the approximate equivalence between the fitted monodisperse radius and the best-fit mean of the distribution of polydisperse radii (Fig. 4, inset).
The freely varying R i , used to shift F (ζ ) from that which is predicted from Eqs. ( 5)- (8), is justified when we consider the number density of individual model bubbles in the system.In a polydisperse pack of spherical droplets at the initial porosity φ i , the droplet (or particle) number density is Similarly, the number density of model bubbles that fit between the droplets is N b = 3φ i /4π a 3 i .Note that N p and N b are defined as number density per unit total volume, rather than per unit liquid volume, as was the case for the number density N introduced earlier.The theoretical solution for the bubble size distribution F (ζ ) via Eqs.( 5)-( 8) finds the maximum size of spherical bubbles that would fit between the spherical droplets.In reality, however, the pores are irregular; hence they have a larger volume than the largest spherical bubble that they can contain.In general, therefore, we expect that N b will be overestimated by this method.We note that this provides a physical explanation for the overestimate of the rate of sintering predicted by the vented bubble model when F (ζ ) is used-the model bubbles are too small, hence also too numerous.The magnitude of the overestimate of N b will depend on the geometry of the droplet packing, which is a function of its size distribution.The droplet number density also depends on this size distribution, and provides a convenient, and conceptually intuitive parameter by which to nondimensionalize the bubble number density.We adopt N b /N p as a metric by which we scale the ratios a i / a i and a i / a i in order to provide a quantitative correction to the vented bubble model for systems of sintering droplets.
The data show that both the best-fit monodisperse a i and the mean of the best-fit distribution of polydisperse radii a i are approximately equal to the predicted value of a i for highly polydisperse systems of droplets.However, there is an apparent deviation from unity of both a i / a i and a i / a i with increasing N b /N p (i.e., for more monodisperse distributions).This suggests that the number density and characteristic radius of the model vented bubbles are well predicted by the polydispersivity model presented in Sec.II B only for highly polydisperse droplet populations, but that scaling is required for more monodisperse populations.We propose an empirical scaling in which a i / a i and a i / a i are exponentially dependent on N b /N p .We fit for the two parameters in the exponential form used and find that a i / a i = simple empirical relationships can be used to correct the predictions of the vented bubble model by correcting the overestimation of the bubble number density derived from the prediction of the bubble size distribution.
While this adjustment is purely empirical, scaling on the basis of a number density ratio is conceptually robust.In summary, we propose that in the polydisperse limit, the surfaces of the pore spaces between the droplets approach smooth bubbles as smaller droplets occupy the narrow clefts between larger droplets.Thus it is not unexpected that our fully polydisperse vented bubble model accounts for the densification dynamics well in this limit, even without adjustment.By the same logic, in the monodisperse limit, it is not surprising that the interdroplet pore spaces are poorly approximated as spherical bubbles; hence our micromechanical model based on bubble geometry breaks down.
As the sintering process approaches completion (at high tb ), the porosity evolution diverges from both the adjusted monodisperse and adjusted polydisperse models.We propose that this is a manifestation of a second-order transition from the permeable to the impermeable state as φ → φ c .This transition involves the evolution of pore clusters such that, at high φ, there exists a single permeable pore cluster, but as φ → φ c , the cluster number density increases.Our experimental design did not permit assessment of the time dependence of the cluster number density.

VI. CONCLUDING REMARKS
We develop a robust microphysical model for sintering of polydisperse populations of droplets under arbitrary nonisothermal conditions.This constitutes an important general-ization of previous models since many natural and industrial applications involve nonisothermal sintering of polydisperse droplets including high temperature terrestrial environments such as volcanic interiors [2][3][4][5], sticking of atmospheric dust or volcanic ash particles in gas turbines [30][31][32], fabrication of ceramics [8,9], and coalescence of droplets in planetesimal formation [1].We validate our model against experiments on populations with a wide range of polydispersivities, and find that an empirical tuning is required to describe the sintering of distributions of droplets that approach monodisperse, but that the unadjusted model is a very good description of polydisperse droplet sintering dynamics.
We find that the model breaks down in the final stages of sintering, as φ → φ c , indicating that further work is required to explain the microphysical behavior of sintering systems near the percolation threshold.Investigation of the progressive isolation of clusters of connected bubbles as φ → φ c , and exploration of the congruence of dynamic percolation transitions with those found for geometrical simulations [33,34], are identified as frontier topics.
We provide raw and processed experimental data [35].

FIG. 2 .
FIG. 2.Comparing the vented bubble model with the experimental data for (a) the monodisperse model and (b) the polydisperse model.In both cases the solid lines refer to the model using a i calculated from Eqs. (5)-(8) to compute the evolution of φ, and dashed lines refer to the model where a i (tuned monodisperse) is a fitting parameter in the solution of φ(t) or where a i is the iterated result of the tuning of the polydisperse model (tuned polydisperse; see Fig.3).

FIG. 3 .
FIG. 3. A worked example for a highly polydisperse droplet size distribution (S ≈ 0.4).(a) The droplet size distribution with the mean droplet size R i labeled.(b) The resultant pore size distribution as calculated with the measured mean of the pore sizes labeled as a normalized value a i / R i .(c) The prediction of the sintering model using the distribution of initial pore sizes.The color scale refers to the monodisperse solutions for each contributing pore radius in the polydisperse distribution.The black curve is the polydisperse solution and the data points are as measured.(d) The droplet size distribution as in (a) but for which the mean is shifted to the best-fit value R i .(e) The resultant distribution of pore sizes using the shifted R i and the new mean a i labeled.(f) As in (c) but using the new distribution of pore sizes in (e).

FIG. 4 .
FIG. 4. Comparison between the mono-and polydisperse tuned vented bubble model and experimental data.Dashed line shows results of the monodisperse model [Eq.(2) with Eq. (4) where a i is a fitted length scale].Solid line shows the polydisperse model using Eqs.(2)-(5), for which a i is calculated using a best-fit Ri , as described in the main text.The normalized gas volume fraction φ or φ

1 .
FIG.5.The scaling factor between the best-fit characteristic radii controlling the coalescence of many viscous droplets and the ratio of the number density of model bubbles to the number density of droplets (see text).

TABLE I .
Scales used to model experimental data.