Topological inversions in coalescing granular media control fluid-flow regimes

Fabian B. Wadsworth,1,* Jérémie Vasseur,1 Edward W. Llewellin,2 Katherine J. Dobson,1,2 Mathieu Colombier,1 Felix W. von Aulock,3 Julie L. Fife,4 Sebastian Wiesmaier,1 Kai-Uwe Hess,1 Bettina Scheu,1 Yan Lavallée,3 and Donald B. Dingwell1 1Earth and Environmental Sciences, Ludwig-Maximilians-Universität, Theresienstrasse 41, 80333 Munich, Germany 2Department of Earth Sciences, Durham University, Science Labs, Durham, DH1 3LE, United Kingdom 3School of Environmental Sciences, University of Liverpool, Jane Herdman Building, Liverpool, L69 3GP, United Kingdom 4Swiss Light Source, Paul Scherrer Institut, 5232 Villigen PSI, Switzerland (Received 30 April 2017; published 25 September 2017)


I. INTRODUCTION
Sintering of viscous droplets is an important process in both natural and industrial scenarios.These include planetesimal formation [1], welding in shallow volcanic interiors [2][3][4], formation of deposits of large meteorite impacts and volcanic supereruptions [5], adhesion of volcanic ash in commercial jet engines [6,7], and industrial fabrication of ceramics and glasses [8][9][10][11].In most sintering scenarios of practical interest, the viscosity of the droplets is high and the droplets are small, such that the Ohnesorge and Eötvös numbers are high and low, respectively.In such a regime, sintering behavior is dominated by the balance between the capillary stresses that drive flow, and the viscous stresses that arise as a result.When the geometry of the sintering system is simple, as in the case of coalescing pairs [12,13] or small clusters [14] of droplets, it is possible to calculate explicitly the evolution of the surface area, and consequently of the driving capillary stress.In such systems the sintering time scale is found to scale with the initial droplet radius.However, for more complex heterogeneous and random systems of many droplets [15], the change of the surface area is not easily calculated, leading to the development of approximate descriptions of sintering, based on bulk geometric parameters such as pore number density and bulk porosity [2,8,16,17].In those systems the sintering time scale is found to scale with the droplet size at early times, consistent with the models of simple droplet coalescence, but with the size of the interdroplet pore spaces at longer times, when the pore space is approximated as spherical bubbles in liquid shells [16,17].We infer, therefore, that this change in scaling represents the topological inversion-from spherical droplets surrounded by gas, to spherical bubbles surrounded by liquid.* Fabian.wadsworth@min.uni-muenchen.deAttempts to unify these descriptions have been made [8], but without a quantitative view of the changing internal geometry of such sintering materials, it is impossible to reconcile this hypothesized inversion with quantitatively constrained dynamic models.Here, we exploit recent technical developments in high-temperature, high-speed x-ray tomographic imaging to demonstrate the inversion of topology implicit in all models of viscous sintering, and show that this manifests itself as a measurable transition in fundamental bulk system properties.

II. MATERIALS AND METHODS
We pack spherical glass beads of mean radius R i = 41 μm to form freestanding cylinders of ∼2.98 mm height and ∼3 mm diameter.These are mounted on the synchrotronsource tomography beam line (TOMCAT) [18] and heated by lasers [19] at ∼5 K min −1 to isothermal temperatures T 0 above the glass transition interval T g , yielding highviscosity droplets; experiments at different T 0 yield different viscosities.The glass beads do not lose mass, and do not show liquid-liquid immiscibility or crystallize on the experimental time scales or at the temperatures tested, which has been confirmed previously by differential scanning calorimetry and thermogravimetry [2].The glass transition onset, determined on samples of the same size and at 5 K min −1 , is T g = 810 K (where we take the onset of the transition as a threshold temperature; see Fig. 1).The temperature dependence of the viscosity (used in the nonisothermal t; see below) and the particle size distribution used here are given in Fig. 1 for reference.
Each three-dimensional (3D) dataset (1440 projections, 1.9 μm voxel resolution) was acquired in 10.7 s.The error on t dominantly represents the 10.7 s taken to collect a 3D tomographic dataset.The spatial resolution of the x-ray computed 3D tomography is 1.9 μm, which permits accurate quantification of pores 6 μm radius [20], smaller than FIG. 1.The physical properties of the glass beads used in this study.(a) The relative heat flow rate (arbitrary units) from which we find the glass transition temperature determined at 5 K min −1 heating rate in differential scanning calorimetry.(b) The temperature dependence of liquid viscosity of a dense liquid of the same composition and the parametrization by the Vogel-Fulcher-Tammann equation form and a comparison with published viscosity data for this material [2].(c) The size distribution of glass beads, which become viscous droplets at high temperature, and the predicted pore size a i .We mark the mean droplet radius R i .the smallest pore in the final isolated spherical pore size distribution.Image visualization and analysis is performed using AVIZO.Segmentation is performed using a standard gradient-based algorithm using the moments of the intensity distribution.All pores with volumes <125 voxels are discarded from the analysis as below this value the error on absolute volumes exceeds 5% [20] and these objects only comprise ∼0.01-0.02vol % of the sample.Pore volumes were calculated for remaining segmented features.
During the time-resolved high-temperature scanning, droplets interact and sinter, and manifest as bulk densification of the compact at the expense of the pore gas volume.The pore and liquid phases (Fig. 2) are segmented from the central region of each sample to avoid edge effects.

III. RESULTS AND ANALYSIS
The tomographic data (Fig. 2) show that the droplets are initially spherical and the interdroplet pore phase is a continuous complex network.As time at high temperature progresses, the pore and liquid phase geometries invert, so that the pore phase is finally a suspension of spherical gas bubbles in a liquid continuum.The in situ nature of the experiment permits us to distinguish the components of the pore phase that are percolating across the system length from those that are isolated and nonpercolating (respectively, gray and green volumes in Fig. 2).This reveals an additional sintering phenomenon that is not accounted for in theoretical models, namely: progressive partitioning of the pore phase, from continuous and fully connected, to discontinuous and fully isolated.This evolution is quantified in Fig. 2(e), which tracks the pore fraction that is connected across the domain φ p , normalized by the total pore fraction φ, as a function of dimensionless time t = t/λ.Here, λ = lμ/Γ is the capillary time, where l is a characteristic length scale, μ is the droplet viscosity, and Γ is the interfacial tension (in detail we use a nonisothermal extension of λ which accounts for the variation in viscosity during heating [16] and the polydispersivity of the initial droplet size distribution [21]).This nondimensionalization, taking l as the mean initial size of the pores between the droplets [16,17], collapses the data for experiments at different viscosities such that a deviation from fully connected (φ p /φ = 1) occurs at a consistent critical dimensionless time t ≈ 0.25.We term this transition 1 .A second transition, at which the pore phase becomes fully isolated (φ p /φ = 0), occurs at a consistent critical dimensionless time t ≈ 1.4.We term this transition 2 .We hypothesize that these transitions, which are microstructural in origin, influence bulk sintering dynamics and associated permeability evolution.

A. Sintering kinetics
We first consider bulk sintering, by comparing data for φ(t) with two end-member models: an extension of the Frenkel model [10] (FM), which considers the dynamics of neck formation at the contact point between coalescing droplets; and the vented bubble model [17] (VBM), which approximates the pore space as an array of spherical bubbles in a liquid shell, which can freely vent their gas.The FM has been extended for bulk sintering [8,22] to give [16] where φ i is the initial pore fraction, φ = φ/φ i , and t takes l = R i in the definition of λ, where R i is the initial droplet radius.Figure 3(a) shows that the FM [Eq.( 1)] captures the early part of bulk sintering dynamics, consistent with its physical basis in neck formation at droplet contact points [8,16].We observe that the model breaks down close to the transition 1 .
The VBM is [16,17] d φ where, in this case, t takes l = a i in the definition of λ, where a i is the initial bubble radius between the droplets (predicted using a model for pore-space geometries between spheres [15]).Figure 3(a) shows that the VBM [Eq.( 2)] is almost indistinguishable from the FM at an early time, before 1 , but captures better the evolution of φ at longer times.There remains, however, a discrepancy between the VBM and the experimental data after the 1 transition.Moreover, the discrepancy grows toward the second transition 2 , at which the pore phase becomes fully isolated, and φ reaches a fixed, finite value that is not captured by either model.Our time-resolved experimental method permits us to measure the connected pore fraction φ p as a function of time [Fig.2(e)].The connected pores are, by definition, the parts of the pore space that are able to continue to vent their gas out of the system beyond 1 .It is reasonable, therefore, to replace φ with φp in Eq. ( 2), where φp = φ p /φ i .Figure 3(b) demonstrates that this modified VBM is in excellent agreement with experimental data throughout the entire sintering process.
Based on this analysis, we draw the following conclusions regarding the role of the transitions in sintering dynamics.First, the failure of the droplet-based FM beyond 1 , taken together with the success of the bubble-based VBM, leads us to associate this transition with the topological inversion from spherical droplets surrounded by gas, to spherical bubbles FIG. 3. Models of sintering dynamics for different microstructural topologies.(a) The evolution of φ with t showing that the droplet geometry model breaks down at 1 and that neither model captures the data after 2 .Here the Frenkel model [12] uses R i in λ to find t and is termed "droplet geometry," while the vented bubble model [17] uses a i in λ to find t and is termed "bubble geometry."Inset: a schematic of the controlling length scales.(b) Demonstration that, when only the interconnected clusters of pores φp are considered, the model cast for bubble geometry describes the data very well throughout both 1 and 2 .Inset: the data for φp normalized by the model result φ * p , for the droplet geometry (top) or the bubble geometry (bottom).surrounded by liquid.Second, the requirement that we replace φ with φp in Eq. ( 2) beyond 1 , combined with the direct observations encapsulated in Fig. 2(e), lead us to associate the onset of pore-space isolation with the same topological inversion.Finally, the agreement of the modified VBM with the experimental data around 2 is consistent with this being the percolation transition at which the connected porosity drops to zero and sintering halts.In our experiments this transition occurs at a critical pore fraction φ c = 0.025 ± 0.009, consistent with ex situ estimates [16] and numerical simulations [23].

B. Permeability of the pore space between packs of sintering droplets
We now explore the impact of the transitions on permeability evolution during sintering.The fluid permeability is simulated through the pore network internal to the densifying liquid using the LBFLOW code [24].This approach applies a simulated pressure gradient across the fluid within the pore network.Average flow velocity is determined when steady flow is reached, defined as when a convergence criterion is met such that the velocity of the fluid averaged over the entire fluid lattice does not vary by >10 −5 (simulation units) over 50 time steps.Low Reynolds number, low Mach number, Stokes flow, and scale independence are all obeyed (Fig. 4; see Ref. [24] for more information).An example of the evolution of average velocity for one experiment is given in Fig. 5 as porosity evolves during sintering for a constant pressure gradient.The subvolume of the sample selected for simulation is 475-μm-square edge length (Fig. 2).
At high φ, before the topological inversion identified with 1 is reached, the LBM data are well described by the Stokes permeability [15] for flow around spheres and FIG. 5.The absolute evolution of the average fluid velocity in each direction in a sample during sintering.The example shown is for T 0 = 908 K and the colors correspond to the evolution of φ or, equivalently, of t.Each simulation is therefore for a single 3D dataset recorded from the x-ray tomography.(a-c) The velocity in each of three orthogonal directions.FIG. 6.The topological inversion is manifest as a shift in the bulk permeability of the system from the droplet geometry scaling to the bubble geometry scaling.(a) Validation of the scaling laws for permeability k normalized by the Stokes permeability k s to give k as a function of the gas volume fraction φ which in the low-φ limit is affected by the emergence of a percolation transition φ c at 2 .Shown are the dilute expansions [25] and simulations using dilute cubic packs (SC-simple cubic; BCC-body-centered cubic; FCC-face-centered cubic) from close to the single-sphere Stokes solution at φ → 1 down to moderate φ.Simulations using polydisperse droplet size distributions identical to the droplet size distribution used in the time-resolved experiments agree with the cubic packs in the dilute regime.The simulation data are normalized using k = k/k s .The scaling in the concentrated regime is simply k = (φ − φ c ) ē where ē takes the theoretical value 4.4 (Ref.[27]).expansions thereof [25] for cubic lattices [Fig.6(a)], consistent with the approximately spherical geometry of the droplets.In the high-φ limit, the Stokes permeability approaches the Stokes limit where R n is the nth moment of the sphere radius distribution, allowing application to polydisperse sphere packs [15].This provides a convenient normalization for calculated values of permeability, such that k = k/k s and is used in Fig. 6(a) for the simulation data far above 1 .At lower φ, below the 1 topological inversion, these expansions for flow around spheres do not describe the data well.Instead the data collapse toward an apparent asymptote where k → 0 at φ = φ c , which is reached at the percolation transition 2 .In the regime between the two transitions 1 and 2 , the data are well described by a percolation permeability k = k r (φ − φ c ) ē, which is consistent with the bubble geometry of the pore space.Here k r is a reference permeability, ē = 4.4 is a theoretical percolation exponent [26,27], and the reference permeability takes the form [26,28] We use s as the specific surface area of the internal liquid-gas interface measured from the 3D datasets.All experimental data in Fig. 6(a) are normalized so that k = k/k r .
This analysis allows the topological transition 1 to be additionally defined as the point where k(φ) diverges from the expansions of Stokes theory, toward the power law percolation theory (Fig. 6).This has intuitive physical meaning when compared with the established change from droplet to bubble geometry (Figs. 2 and 3).The densification halts as φ → φ c /φ i , below which there is no pore cluster that percolates to the outside of the sample [see Fig. 2(e)].The success of permeability scaling (Fig. 6) means that we can combine k/ ki (φ p ) with φ p ( t) modified from Eq. ( 2) to predict the time evolution of permeability during sintering, even when the interdroplet pore space undergoes this topological transition.

IV. CONCLUDING REMARKS
We explore the kinetics of sintering of droplets in light of time-resolved 3D high-temperature reconstructions of the microstructure of the process in complex packs of droplets.This provides a quantitative prediction of the evolution of pore-space permeability throughout sintering useful for industrial applications of ceramic filter production and for better understanding gas flow through welding volcanic particulate deposits.Outstanding questions remain concerning the quantitative prediction of the time dependence of porespace isolation.Solutions to this problem would make these observations of wide utility and remove the necessity for in situ techniques to find connected pore volume fractions.

FIG. 2 .
FIG.2.The internal structure of a pack of many interacting droplets as they sinter to volume equilibrium.(a-d) 3D renderings of the pore space between droplets initially packed to an interdroplet pore fraction φ i = 0.46 (where connected pore space is rendered in gray and isolated in green) at (a) t = 10 −2 , (b) t = 10 −0.3 , (c) t = 10 0 , and (d) t = 10 1 , at which time φ ∼ φ c .The inset silhouettes show two-dimensional (2D) cross sections of the bulk sample set on a ceramic plate with the green region of interest (475 μm edge) shown on which rendering and quantitative analysis were performed.(e) The proportion of the pore fraction that is connected to a sample edge φ p /φ as a function of dimensionless time t showing the two transitions observed here: 1 and 2 .

FIG. 4 .
FIG. 4. The dependence of the fluid permeability on the volume or position of the simulation domain, or equivalently, the number of fluid nodes.(a) Showing convergence when the domain exceeds a threshold size (shown here for a sample of φ = 1 and a driving pressure gradient of ∇P = 10 −2 Pa m −1 ).For samples of φ < 1, the scale dependence becomes less significant.Tested here are multiple positions of the subvolume in the larger sample volume.(b) The effect of the driving pressure gradient the resultant permeability shown here for the largest domain size used (L 3 > 10 7 voxels).

1
represents the approximate topological inversion point at φ ∼ 0.29 and 2 occurs at the percolation transition observed at φ c = 0.025.Permeability is modeled using the time-resolved experimental data in three orthogonal directions k xx , k yy , and k zz .The data are normalized using k = k/k r (see main text).(b) The combination of the scaling in (a) with the kinetics of φp in Fig. 3 permits the models for bubble geometry and droplet geometry to be cast as a function of t as k/ ki = φ ē p .(c-f) The interdroplet fluid velocity vector distribution at the same dimensionless times as in Figs.2(a)-2(d) simulated using a lattice-Boltzmann technique, where warm colors refer to high velocities under a driving pressure gradient.Here visualized are the x-direction velocity u x vectors.
These data agree with the model prediction in which k/ ki ∼ φ ē p [Fig.6(b)].