Mesostructural investigation of micron-sized glass particles during shear deformation – An experimental approach vs . DEM simulation

The interdependency of structure and mechanical features of a cohesive powder packing is on current scientific focus and far from being well understood. Although the Discrete Element Method provides a well applicable and widely used tool to model powder behavior, non-trivial contact mechanics of micron-sized particles demand a sophisticated contact model. Here, a direct comparison between experiment and simulation on a particle level offers a proper approach for model validation. However, the simulation of a full scale sheartester experiment with micron-sized particles, and hence, validating this simulation remains a challenge. We address this task by down scaling the experimental setup: A fully functional micro shear-tester was developed and implemented into an X-ray tomography device in order to visualize the sample on a bulk and particle level within small bulk volumes of the order of a few micro liter under well-defined consolidation. Using spherical micron-sized particles (30 μm), shear tests with a particle number accessible for simulations can be performed. Moreover, particle level analysis allows for a direct comparison of experimental and numerical results, e.g., regarding structural evolution. In this talk, we focus on density inhomogeneity and shear induced heterogeneity during compaction and shear deformation.


Introduction
Shear flow of granular media is ubiquitous in nature and of industrial importance when it comes to handling and processing of bulk solids.The interesting shear rate dependence of granular flow is accompanied by a rate-independent regime at slow, quasi-static deformation.Strain localization within the bulk, often referred to as failure zone or shear band, represents one unique feature of this regime [1].A crucial task is the prediction of the critical shear stress a resting bulk solid can withstand.This shear stress depends on the preparation-procedure and is also known as yield locus.Utilizing non-associate flow rules, plasticity theory provides a framework to model quasi-static flow.However, in most cases only the volume fraction is taken into account to characterize the structure.
Modern technology offers a variety of powerful tools to further investigate the interplay between structure and mechanical behavior of granular media.On one hand, the discrete element method (DEM) provides detailed information on the behavior of small samples presuming a certain particle interaction.On the other hand, sophisticated experimental setups are capable of determining the properties of micron-sized particles [2] and of picturing the inner structure of a bulk solid non-destructively via computer tomography [3].Using this methodology, we critically ex-plore the potential of DEM simulations and micro-sized experiments by performing a direct comparison considering the bulk structure as well.For this task, a fully functional micro-sized shear tester was developed [3], which can be fitted into an X-ray tomography device.As will be shown, this setup allows for an in-situ shear band detection during an experiment at constant normal load with only a few μl of sample volume.

Experiments
For this study we used micron-sized particles made of borosilicate glass microspheres (CoSpheric LLC, USA).The glass spheres are characterized by a narrow monodisperse particle size distribution with median value x 50,3 ≈ 30 μm.Coequally an almost spherical shape of particles is observed.Both features benefit the subsequent computational image-based analysis of tomographic data and ease the comparison to the DEM simulation using ideal spheres.Besides, a non-negligible surface roughness and a slightly cohesive behavior is identified.

X-ray microtomography and micro shear-tester
In order to measure the spatial powder structure and behavior, a nondestructive image-based analysis is realized by using a high-resolution X-ray microtomography device (MicroXCT-400, Zeiss (Xradia), Germany).The operating parameter of the X-ray source amounts to an acceleration voltage of 50 kV with a constant current intensity of 200 μA.This choice of parameters results in high contrast images, and hence, supports the accurate computational image-based analysis afterwards.Considering a shear cell diameter of 2 mm, a ten-fold optical magnification was adequate to image the entire cell with a resolution of 2.2 μm (2x binning).In addition, an intensity gray-scale projection is formed by detecting X-rays passing through the sample.A record of 2000 projections with varying sample orientation provide the basis for a precise three-dimensional reconstruction and smooth the way for a following detailed computational sample analysis.
The experimental setup contains a torsional micro shear-tester with a cylindrical shaped sample chamber [4].In comparison to common shear devices it allows for investigating very small volumes on the order of 6 μl and features both powder densification up to 20 kPa and shear deformation of small angle increments from millidegree up to theoretically unlimited shear displacement.Vertically, the sample chamber is confined by a pair of structured pistons whereas a filigree glass capillary with wall thickness of 50 μm acts as a sidewall.Another characteristic feature of this modern shear device is the decoupled and most precise determination of the acting forces and torque behavior whereas the whole system is theoretically frictionless using advanced technology of magnetic spring and air bearing mechanisms.
In order to prevent agglomeration during the preparation process, powder is sieved into the capillary.After the filling process, the sample is compressed to the desired normal load of 0.5 kPa.During shear deformation the sidewall and the upper piston rotates stepwise in angle increments of 0.5 • up to 9.5 • shear deformation.Afterwards, the angle increment is raised up to 5 • until the final shear deformation of 39.5 • is achieved.The shear deformation is performed in a quasi-static regime using an angle velocity of ω = 0.1 • /s.Between the shear steps tomographic data is recorded.

Computational analysis -Image segmentation and particle tracking
The advanced computational analysis of tomographic data bridges the gap between experiments and DEM simulations regarding particle positions, and hence particle movement over time.First, particle positions (and radii) are extracted from the image data by using a marker-based watershed transformation [5].In this context, a Gaussian filter and the IsoData algorithm from ImageJ [6] are applied to smooth and binarize the image data, respectively.Besides, due to few hollow particles, small disconnected pores are filled with solid artificially.On the basis of particle centroids calculated from these segmented binary images, a particle tracking is realized.Several time steps with a small angle increment of 0.5 • are used to apply a tracking algorithm that approximately minimizes the sum of squared displacements for each time step [7].To speed up calculations, particle displacements which exceed a predefined threshold are not considered to be possible during optimization.In this study the tracking threshold is set to s = 19.8μm.This simplifying assumption is valid for small shear angle increments of 0.5 • , where the largest particle displacements between two time steps are expected to be 9 μm (at the capillary sidewall).However, it does not hold for later time steps of the experiment, when shearing was performed in steps of 5 • .Here tracking is still possible, but it requires prior information.For this purpose, we estimate a local angle of rotation directly from the image data as a function of the radius and height in the sample.The basic idea for this approach, which has been proposed in [4], is to rotate each slice of an image stack by a number of different angles and to determine for which angle it agrees best with the corresponding slice of the subsequent image stack.We refined this approach by subdividing each image slice into disjoint rings of equal area which are rotated independently, and captures variations in radial direction.The quality criterion for comparing two image slices is their correlation.Using this image-based local shear deformation and an adaption of [7], we can track a large number of particles even in steps of 5°.We validated the method based on the time step from 5 • to 10 • of shearing, where information in small steps is available, and found that 98.9 % of the tracks identified by our method agreed with the tracking based on 0.5 • steps.

Contact model
A soft particle DEM approach is utilized to simulate the system of interest, i.e. trajectories are calculated by numerically solving Newton's equations of motion.Particles interact via pairwise forces ( F i j ) and torques, which, in general, depend on their relative position ( r i j = r i − r j ), velocity, angular velocity, diameter (d i , d j ) and the contact history to account for friction.In consideration of the external load, we use a damped Hertz contact model [8] and add adhesion as well as rolling and torsional friction.The normal part of the contact model in dependence of the particle overlap ξ is given by with F adh being a constant, ) the effective elastic modulus, γ n a viscous constant and r eff = 1 2 d i d j /(d i + d j ) the reduced radius of the particles i and j.
Although a linear dependence of adhesion force on r eff is employed, we expect minor influence due to the narrow particle size distribution.
In order to implement Coulomb friction, together with rolling and torsional friction, we use the approach presented in [9], which suppresses each contact mobilization mode (sliding, rolling and twisting) individually by using a linear spring-dashpot, until a threshold is exceeded.This threshold, as well as the resisting force (or torque) in case of mobilization, is proportional to |F n + F adh |.For simplicity, the same friction coefficients are used in the static and dynamic case.

Parameters and calibration
In principle, all microscopic parameters could be determined by sophisticated experiments [2].However, the cumbersome work of measuring every parameter experimentally can be circumvented by utilizing microscopically undetermined parameters to calibrate the contact model.This enables a correct description of the bulk's macroscale behavior [10], e.g.stress response to applied strain.Is it possible to predict the structural behavior, e.g.Poisson's ratio ν as well?We pursue this task by intentionally keeping the contact model and the parameter choice as simple as possible, exploring its predictive power as well as its limitations.
Particle stiffness and adhesion force are set according to the experimentally determined values.As we study the bulk's response to a quasi-static deformation, viscous parameters should have a negligible influence.For computational convenience, we use the high damping coefficient γ n = √ 4/3 E eff m eff , where m eff is the reduced mass of the colliding particles.A typical choice for the tangential stiffness is k t = 2/7k n , where k n (ξ) = ∂F (el)  n /∂ξ with F (el)  n being the elastic part of equation ( 1).Here we use k n = 4/3E eff x 50,3 /2, approximating the contact radius with the median particle radius and therefore overestimating the stiffness.All springs utilized for friction are damped, using a viscoelastic constant of 2 √ k t m eff .As twisting and sliding are closely related frictional modes, the torsion friction coefficient, μ tor /r eff , is set equal to the Coulomb friction coefficient μ.Furthermore, a load dependent rolling friction coefficient μ rol = 2r eff ξ is employed.This choice of parameters leaves the Coulomb friction coefficient, μ, as a free parameter to calibrate the contact model.
We model the experimental particle size distribution with a log-normal distribution and draw particle diameters according to the corresponding probability function.Only diameters between 25 μm and 50 μm are accepted, resulting in a slightly sharper cut-off at small particle sizes and an interquartile distance of x 75,3 − x 25,3 ≈ 5 μm, in comparison to experiment x 75,3 − x 25,3 ≈ 10 μm.
Calibration of the contact model is done by executing plane shear simulations under variation of μ at a fixed normal stress σ = 15 kPa and comparing the steady state macroscopic friction coefficient μ macro to the experimental findings.These calibration simulations are done with N = 10.000particles, confining walls perpendicular to the shear direction and periodic boundary conditions in all other directions (see [3,10] for details).Fitting an exponential saturating function for μ macro (μ) as suggested in [11] results in μ ≈ 0.58.

Setup and preparation
The simulation setup mimics the experimental micro shear-tester: A cylinder with diameter of 1.78 mm is used to confine the sample in horizontal direction, while two pistons (one movable, one fixed) close the simulation cell in vertical direction.Both pistons are structured similar to the experiment, i.e., with six vanes arranged in a star shape.As there is a small gap between vanes and outer cylinder in the experiment (≈ 0.025 mm), we introduce a finite distance between the outer edge of the piston and the cylinder wall of ≈ 0.1 mm.This slight gap enlargement is introduced to avoid jamming of single particles between the piston and capillary.Particle-wall interaction is described by the visco-elastic part of the contact model Equation (1), using the same material parameters.
The preparation procedure includes the following steps: A random, porous particle configuration with volume fraction ≈ 0.2 is created.Potential overlaps are removed in a relaxation run that uses a viscous background friction.After relaxation, this omnipresent viscous term is removed and a normal load of 0.5 kPa is applied to the lower piston, compressing the material.Following compression, a constant angular velocity ω is added in order to shear the sample.To achieve a quasi-static deformation, we use the criterion ω ρ s x 2 50,3 /σ ≈ 10 −4 .

Results and discussion
We compare density profiles of the initial (compacted) and final (sheared) configuration as well as the temporal evolu- sim.
-0.8 -0.6 -0.4 -0.2 0.0 tion of shearband position and width.Assuming rotational invariance, all quantities are given as a function of height (z) relative to the vanes top edge and distance to the central axis of the chamber (r).
The number density n(r, z) of particle centers in axial and radial direction is displayed in Fig. 1.Due to the different filling height, only the range 0 mm ≤ z ≤ 0.7 mm was considered for the radial distribution.
Densities in experiment and simulation are of the same order of magnitude.The initial states differ to much higher extent due to different filling procedures and configurations for experiment and simulation.In contrast, after steering, the differences between experiment and simulation are quite small, showing that the steering can be well described by the simulation.However, the number density increases with identical gradient according to height in experiment and simulation.
This gradient reflects the outer wall's influence, which seems to be equal at first sight.However, the radial profile reveals a denser region in the chamber's center in the experiment, which is even more compacted during shear.Simulation results show the opposite: Slightly increasing density in the vicinity of the sidewall and dilation during shear.Of course, in both experiment and simulations, a dilation zone (z ≤ 0.2 mm) as well as a compaction zone (z > 0.2 mm) can be identified.Layering is observed in both experiment and simulations.Besides wall influences, the systematically higher density suggests that the contact model is not able to reproduce the structure, in terms of volume fraction, fully correct.A model refinement addressing this issue is work in progress.
The shear-induced dilation zone (z ≤ 0.2 mm), which can be identified in the axial density profile, corresponds to the shear band position.This can be shown by analyzing particle displacements.We determine the shear-band position (z sb ) and width (w sb ) by fitting [1] to the normalized angular velocity at fixed radial distances r to center of rotation, which deviates slightly from the central axis in the experiment.The temporal evolution of both quantities, averaged over r, is depicted in Fig. 2. The Inset shows the shearband position as a function r in the steady state.Deviations in shear band width at small deformation might be correlated to the experiment control, while the deviation in the radial dependence of the shearband position is caused by the larger gap used in simulations.Despite these differences, the shear band width (≈ 8x 50,3 ) and position is identical and seem to be robust with respect to details of the contact model parameters.
Note that the number density inside the shear band is also equal in experiment and simulation.
In this study a correlation between local density reduction and a shearband formation above the lower piston is verified.Moreover, we could image the temporal sequence of shearband formation on a particle level using an efficient tracking algorithm that allows for a direct comparison of micron-sized glass particles in experiment and DEM simulations.Besides comparing mechanical properties of the bulk, this structural analysis increases the validation depth for the chosen contact model.

n]Figure 1 .
Figure 1.Axial and radial density profiles of the initial and final configuration.The density gradient is represented by the dotted lines.z = 0 marks the vanes top edge