Consolidation Behavior of Inhomogeneous Granular Beds of Ductile Particles using a Mixed Discrete-Continuum Approach

In this paper we study the consolidation behavior of inhomogeneous granular beds of elasto-plastic particles by means of a discrete-continuum formulation, which systematically bridges the micro and meso scales. The methodology is particularly suitable for describing the post-rearrangement regime where consolidation proceeds mostly by elastic and inelastic deformation. This formulation is able to provide the quantitative estimates of the evolution of macroscopic variables, such as pressure and density, while following microlevel processes, such as local coordination number and loading paths. This methodology is used to simulate binary mixtures composed by particles with different nonlinear properties. The predictions are in general agreement with the experimental data during both loading and unloading cycle.


Introduction
Consolidation of granular materials is a complex process operated by several mechanisms including particle rearrangement, elastic and plastic deformation, fracture and pulverization [1]. The particle rearrangement is a void-filled-in process which is characterized by minute driving or applied forces and by large non-affine particle motion. In contrast, the post-rearrangement stage or compaction stage is identified by the inability of the particles to reduce the volume by relative motion, requiring larger driving forces to deform or fracture the particles for further consolidation [2,3]. While for brittle particles fracture is the dominant mechanism of consolidation Ҁ after the initial and localized elastic deformation stage Ҁ, for ductile particles inelastic or plastic deformation provides a source for reducing the interparticle porosity, and consequently the total volume. For systems dominated by frictional forces, rearrangement and compaction tend to operate simultaneously at different spatial locations of the specimen. In those cases, molecular dynamics methods [4] or discrete element methods (DEM) [5] provide a tool to simulate these conditions. However, due to the intrinsic nature of these formulations, which trace each individual particle, only relatively small samples can be effectively simulated. Often the permissible sizes are too small to characterize the macroscopic properties of the material. In the present work, we restrict our attention to thoroughly lubricated systems where interparticle friction is weak. For this case, the compaction may be uncoupled from the particle rearrangement as a follow-up process. In this scenario, the rearrangement is simulated first, followed by the subsequent compaction. A description of the particle rearrangement methodology is briefly described in Section 2. This method is a modified DEM version (see for example [6,7] for traditional DEM formulations) which can simulate large samples with keeping an excellent agreement with the experimental data. The compaction regime is simulated by the application of a mixed discrete-continuum or granular quasi-continuum (GQC) formulation. The main aspects of the model are presented in Section 3 and the details are available elsewhere [8].
The present GQC, which can simulate granular beds of nonlinear spherical (3D) or cylindrical (2D) particles with an arbitrary size distribution, is utilized to investigate the influence of the composition of granular structure and material properties on the macroscopic compaction behavior. To that effect, two groups of particles with a given size distribution are generated and later mixed with different ratios to conform the samples, as described in Section 4. A para-metric study of the effect of particle size, Young's modulus, yield stress and initial microstructure is then presented in Section 4 along with a detailed discussion.

Die Filling and Particle rearrangement
The initial particle bed structure is generated by a ballistic deposition technique, which simulates the die filling process. The algorithm proceeds by dropping individual particles from the container top until they reach a stable position. The condition of stability is dictated by the cohesivity of the powder. Cohesive powders result in open structures, as shown in Fig. 1(a) for a container of 50҂60 mm 2 filled with 0.4 mm-radius monosized particles. Due to the open structure of this powder bed, particle rearrangement operates in the early stages of densification. This process is simulated by forcing the particles residing at the top surface to move down at a specified rate. At each step, the punch is displaced 0.15 mm and the particles are allowed to accommodate until a new equilibrium configuration is attained. Due to the fact that the rearrangement is an evolution process of position change, an iterative scheme is required to establish the updated equilibrium configuration. Fig. 1(b) shows an intermediate stage during the rearrangement process, where the punch has been displaced by ∆҃5.1 mm. Three clear regions are observed: an unperturbed region near the bottom, a packed region near the punch, and a narrow gap in between or rearrangement front. Massive non-affine motion is only observed on rearrangement front in close agreement with the experiment record [9].

Compaction
Due to the inability of the particles to modified their current structure, further punch motion induces particle deformation as a mechanism of volume reduction. During this stage both relative roll and slip are drastically reduced and play a negligible role on the macroscopic response of the granular bed [2,3]. Under these conditions a discrete/continuum approach can be effectively exploited to account explicitly for the interparticle interactions within a constrained kinematic field, as described succinctly in the next section and in [8].

Mixed Discrete/Continuum Approach
The central idea of this approach hinges on the ability of constraining particle motion by a coarse grid of dimensions commensurable with the sample size. In addition, the constrained displacement field can be realized by utilizing standard Finite Element techniques. Within this framework, a mesh is generated, in either 2D or 3D, and the motion of each particle is described by the element that contains the particle's center. Interaction with particles in the same and neighboring elements are allowed via an interparticle or local constitutive relation. It should be noted that this approach is equally applicable to both dynamic and quasi-static conditions, offering similar advantages. In the present case of analysis, we are interested in the evolution of the powder bed under relatively slow loading conditions Ҁ compared with characteristic wave speed in these materials Ҁ, and thus, we concentrate on a quasi-static formulation. Details of the formulation are presented in [8].

Local Constitutive Equations
The response between spherical particles can be divided into five regimes: i ) tension, ii ) local elastic contact, iii ) local elastic-plastic contact, iv) fully plastic, and v) finite particle deformation. Although the present approach can be extended to describe the whole range, the present study is restricted to the first four. The tension regime (cohesion) is described by a potential law with Lennard Jones' functional dependence and a cut-off. The cohesivity is then modulated by the value of the only adjustable coefficient. While the elastic response is governed by a Hertizan law, where F is the contact force, α is the relatively approaching displacement of the particle centers after contact, R and E are the equivalent radius and the equivalent modulus of the contact which are defined by where E 1 , E 2 and ν 1 , ν 2 are Young's moduli and Poisson's ratios respectively of two particles in-contact with radius R 1 and R 2 . The elastic regime only dominates the initial stages of compaction where particle deformation is minimal, and the plastic deformation quickly eclipses the elastic part. If the elastic/ plastic response of the solid materials of particles is described by assuming a linear elastic regime followed by a power-law plastic regime characterized by hardening exponent m (which is set equal to 3 in the present study) where the label i҃1, 2 refers to the particles in contact, the interaction law for these particles is then given by the similarity solution [3] ҃B(m) 2c 2 (m) and no interaction of the contacts need be considered due to the localized nature [3]. The reference yield and the equivalent reference stress Γ is given in terms of the yield stress s i of each particle i. The values of c, k and B for a given exponent are numerically estimated by [3] resulting in the following expressions The two limiting regimes of the local constitutive equation, elastic (m҃1) and elastic/perfectly plastic (m→ȍ), agree well with experimental data. For m҃ 1 the Hertzian regime is recovered for which exists ample experimental validation. For m→ȍ, a linear dependency is predicted by (4) which is in a good agreement with the recent experiments of [10] for nearly ideal plastic materials.
To account for the irreversibility of the plastic deformation, a Hertizan unloading branch is introduced, which allow us to trace elastic recovery during ejection as well as heterogeneous fields where loading and unloading are concurrent. The unloading regime is described by where α res is the residual indentation depth that left in the particle and might be obtained from both (4) and (1). Another indication of the irreversibility of the process is observed on the evolution of α res , which satisfies α · res ͧ0.

Sample Generation
In order to investigate the influence of particle size and material on the compaction behavior, different particle beds are generated by mixing two sets of particles, each one with a different size distribution and material properties. The Group I with a relatively small mean diameter and the Group II with a larger one. The size distribution follows a normal distribution characterized by an average radiusr and a vari-  Lower and upper limits cut-off are introduced to simulate particle populations within two given thresholds. The details of the selection of the size distribution and material properties are provided in the following sections.
For the purpose of this analysis, two cohesive granular beds are generated by mixing Group I with Group II. The mixture A has 30% in volume of particles from Group I and 70% from Group II, conversely, mixture B has 70% of group I and 30% of group II. With these simulated mixtures we fill a frictionless container of 50҂60 mm 2 utilizing the ballistic deposition method described in Sec. 2. The resulting structures are shown in Fig. 2, which share common features, such as voids, with the monosized distribution shown in Fig. 1(a).
The samples are then subjected to consolidation by rearrangement using the simulation approach described in Sec. 2, which provides the initial samples for consolidation by compaction. The post-rearrangement structures are shown in Fig. 3, where the overimposed mesh shows the resolution of the coarse grain simulation using our discrete-continuum approach. It is interesting to notice that the structures formed after the rearrangement of both inhomogeneous and monosized beds show some qualitative and quantitative differences, such as the appearance of grain-like structure, which is only observed in the monosized case. The inhomogeneous cases consistently exhibit a lower post-rearrangement average density, as indicated in Table 2 (compare to the value of 0.867 for the monosized case), indicating a better packing for the monosized case. During the process of rearrangement both the average density, Ҁ ρ, and the average coordinate number, Ҁ N c , increase as shown in Table 2. Also, there is a tendency to reduce the inhomogeneity on the initial or pre-rearrangement density distribution resulting from the die filling process. This effect is clearly indicated in Fig. 3. An interesting observation is that the post-rearrangement density of mixture B is not a direct function of the fraction of fines, in contrast, the average coordination number of the sample shows a direct correlation. This observation is in good agreement with the proposition that only when the size of the fine fraction is small enough to reside at the interparticle spaces of the coarse fraction (no perturbation of the packing),  the densification is improved [16]. Although Ҁ ρ for the inhomogeneous samples are lower than the monosized ones, Ҁ N c is higher. This is consistent with the obser vation that in monosized beds the low number of contacts along the disorder regions (grain boundaries) significantly reduces the average coordination number of the entire region (maximum coordination number is 6 for a perfect hexagonal arrangement). In contrast, in inhomogeneous beds large particles may be surrounded by more than 6 smaller particles, increasing the average coordination number.

Size Distribution
In order to focus on the effect of post-rearrangement structure, the Group I and II are assumed to have the same solid material behavior which is characterized by a Young's modulus of 10 GPa, a Poisson's ratio of 0.33 and yield stress of 100 MPa. The global pressure-deformation responses and unloading paths are plotted in Fig. 4 Table 2 Average density, coordinate number and particle bed region before and after particle rearrangement. to monosized bed is also included for comparison purposes. The pressure p is computed as the applied vertical force per unit of effective area. The effective area is estimated as the multiplication of the width of the container (50 mm) by the weighted average thickness of each particle in the bed, where the particle diameter is used as the weighting factor. The deformation is characterized by the ratio of the punch displacement over the initial height of particle bed (60 mm). An evident conclusion from the figure is that size distribution has small effects on the macroscopic response as indicated by the parallelism among the three curves. The only observed difference (horizontal shift) is due to the disparity of structures induced during the rearrangement process. This observation can be understood with the aid of (1) for the elastic regime and (4) for the plastic one. The average local or particle pressure may be equated to F/πR 2 , which is a function (nonlinear) of the average particle deformation α /R. Thus, the average pressure vs. average deformation is then independent of the particle radius, in accordance with the simulation results. The present approach provides, in addition of the macroscopic pressure-deformation response, the local distribution of the average stresses and their temporal evolution. These distributions are of key importance in tracking regions with potential defects and imperfections. In Fig. 5 shows the spatial distribution of s yy (normal stress in the direction of the compaction) for mixtures A and B, and the monosized recorded for the same applied displacement of the punch. The samples show different distributions with more uniformity for the mixture B where a large number of fine particles present. The non-uniformity in the monosized particle bed agrees with the exis-tence of the heterogeneous structure as grains and boundaries.

Elastic Behavior
The elastic properties affect not only the local stress distribution but also the global pressure-deformation response. In order to investigate the inf luence of Young's modulus on the compaction behavior, we conduct a parametric study by adopting h E ҃E I /E II ҃ [0.1,1,10], where E I and E II are the modulus of the two groups. The global pressure-deformation curves are plotted in Fig. 6. The Young's modulus of the particle group with higher volume percentage controls the average responses, even in the plastic regime. In that regime, (4) shows that the interparticle pressure is proportional to the effective stress Γ which in turns relates to the Young's modulus by (5). In the elastic regime (initial loading or unloading), (1) and (8) provide a direct connection between pressure and Young's modulus.
The local stress distributions are highly dependent of the elastic properties as shown in Fig. 7, where the contour of s yy is plotted for both mixture A and B. The stress may have different limits due to the variety of Young's modulus, but the f lood contours are drawn by setting the difference of the limits the same as it is done in Fig. 5, which is 50 MPa. At the same punch displacement as in Fig. 5, as expected, the stress distribution is more heterogeneous than the case of uniform material.

Plastic Behavior
To quantify the effect of the interparticle plastic behavior on compaction response, we vary, at constant Young's modulus, the parameter h sy ҃s yI /s yII    Fig. 8, where there is no crossing between the load and unload path on graphs of the load-displacement curves. The local stress distribution is also shown in Fig. 9 for both mixtures and for h sy ҃0.1 and 10. Compared with Fig. 5, it may be found that decrease of yield stress in either of the particle groups will lead to more uniform of stress distribution.

Conclusions
The present mixed discrete-continuum approach provides a methodology to effectively link different length and time scales for nonlinear heterogeneous particle systems. This technique resolves the kinetics at interparticle level while constraining the kinematics by a grid with suprapartcle resolution. This formulation is utilized to study the consolidation behavior of a series of samples after the pre-consolidation by means of particle rearrangement. The process of rearrangement, which consolidates the highly open structures resulting from simulation of die filling of cohesive powders, is simulated using a modified DEM approach. Model systems are studied to quantify the role of the size distribution as well as the elastic and plastic properties. In these studies, different composition of binary mixtures of two particle sets are consolidated to relatively high pressures. Each particle set is endowed with a particle size distribution (characterized by the mean and variance) and the particle's elastic module (E), yield stress (s y ) and hardening exponent (m). The results clearly indicate the significant effect of each of these characteristics not only on the overall behavior of powder bed but also on the development of local stress states. It is interesting to notice that the post-rearrangement preconsolidation structures are significantly affected by the disparity on the size distribution between the 2 particle sets. These structures prevail during the subsequent consolidation giving rise to dissimilar local states even when the overall applied force against punch displacement remains mostly unchanged. This observation shows some noteworthy features of the current methodology: 1) the ability to follow the history of the material before consolidation, which is indicated by the heterogeneity of the local fields and 2) the ability of computing consistent averages of these local states, which is indicated by the similar force-displacement response. Finally, due to the introduction of non-elastic effects the different paths during loading and unloading can be readily traced in accordance with the experimental observations.