Grain-size Effects on Mechanical Behavior and Failure of Dense Cohesive Granular Materials †

The grain sizes can significantly influence the granular mechano-morphology, and consequently, the macro-scale mechanical response. From a purely geometric viewpoint, changing grain size will affect the volumetric number density of grain-pair interactions as well as the neighborhood geometry. In addition, changing grain size can influence initial stiffness and damage behavior of grain-pair interactions. The granular micromechanics approach (GMA), which provides a paradigm for bridging the grain-scale to continuum models, has the capability of describing the grain size influence in terms of both geometric effects and grain-pair deformation/dissipation effects. Here the GMA based Cauchy-type continuum model is enhanced using simple power laws to simulate the effect of grain size upon the volumetric number density of grain-pair interactions, and the parameters governing grain-pair deformation and dissipation mechanisms. The enhanced model is applied to predict the macroscopic response of cohesive granular solids under conventional triaxial tests. The results show that decreasing grain-sizes can trigger brittle-to-ductile transition in failure. Grain size is found to affect the compression/dilatation behavior as well as the post-peak softening/hardening of granular materials. The macro-scale failure/yield stress is also found to have an inverse relationship with grain-sizes in consonance with what has been reported in the literature.


Introduction
The macro-scale behaviors of dense cohesive granular materials are directly linked to the grain-scale through their granular mechano-morphology, which is represented by their collective arrangements (microstructure) and their interaction mechanisms (micromechnics). The mathematical description of macro-scale response of granular materials must begin from the conception of grain-interactions. From this point of departure, either discrete or continuum descriptions can be elaborated. In discrete models, the macro-scale behavior is inferred by tracking the motion of each grain in response to the external action. Continuum models, in contrast, aim to describe behavior at a spatial scale in which the grains and their motions remain concealed within a material point of a continuous system. While discrete analysis schemes can provide results with detailed description of grain motions, inter-granular deformations, they are limited in their applicability to large-scale problems due to large computational cost, and the challenges of grain-scale verification and validation. On the other hand, classical continuum models, while computationally affordable, suffer from losing characteristic features due to neglecting grain scale phenomena (in terms of both the microstructure and its evolution, as well as the inter-granular mechanisms). In this regards, it is noteworthy that higher-gradient continuum models represent a way to capture the effect of microstructure and its evolution during loading (Placidi L. and Barchiesi E., 2018;Placidi L. et al., 2018a;2018b), an approach that can be traced to the early development in continuum mechanics (dell'Isola F. et al., 2014;Eugster S.R. and dell'Isola F., 2017).
Granular micromechanics approach (GMA) provides a paradigm that bridges the grain-scale models to appropriate continuum model Poorsolhjouy P. and Misra A., 2019). In this approach, the constitutive relation for a continuum material point is derived by utilizing the directional average of variously oriented grain-pairs. In this way, GMA incorporates data from the microstructure and micromechanical phenomena at grain scale into the continuum model . It has been shown that by incorporating the grainpair interactions in different directions, GMA is able to capture various behavior, that are difficult to attain with classical continuum models, using computational demand which is significantly lower than that associated with discrete models . In a form of GMA that leads to Cauchy-type continuum model with coupled damage and plasticity, Helmholtz free energy density and dissipation potential functions for the continuum model are obtained as the volume averages of corresponding functions defined for the grain-pair. Grain scale constitutive equations are calibrated to ensure that the model replicates the observed behavior. Since the grainpair behavior is nonlinear, the grain-pair oriented in different directions undergo different unique loading paths as a continuum material point experiences a macro-scale loading sequence. The consequence is that the response of the granular system is highly loading-path dependent. Indeed, when a macro-scale body is subjected to a loading-path, the variously oriented grain-pair can be in different loading stages including loading, unloading, or reloading. The material response, therefore, at any stage of loading is a function of the entire loading history rather than the current loading state.
Since grain-pair interactions are key in GMA-based continuum models, the grain size effect can be incorporated by considering its influence on the grain-pair behavior. Here, we enhance the non-linear damage-plasticity grainpair interaction previously proposed for cementitious materials by incorporating the effects of grain sizes. We introduce simple scaling relationships, based on power laws, to describe the effect of grain size on grain-pair stiffness, grain-pair damage dissipation, grain-pair branch length, and the volumetric number density of grain-pair interactions in the granular material system. Using these relationships, size-effect relationships for macroscopic response of the material (macroscopic stiffness tensor and the Cauchy stress) are predicted. These results are calculated for parameters that are calibrated using experimental data available in the literature. The predicted results show that grain-sizes can influence brittle to ductile transition and softening-hardening behavior. The results also show that macro-scale failure stress under axisymmetric loading scales with grain-sizes as has been observed from experiments on cohesive granular materials.

Granular Micromechanics Approach (GMA)
GMA can be categorized as a meso-scale approach that bridges the gap between continuum models and discrete models as is seen in Fig. 1. At the smallest scales, atomic models can be utilized provided the composition and initial structures are known. Indeed, provided a faithful realization of the atomic structure, these models can yield results with high degree of exactness (i.e. fidelity). For examples of their application to atomically complex materials the reader is referred to (Misra A. and Ching W.Y., 2013;Qomi M.J.A. et al., 2014). These models however are computationally costly even with the present computation and memory capacities. Moreover, they are inapplicable to large systems with complex compositions, or disordered and defective structures given the difficulties associated with the characterization and specification of such materials. In larger scales, discrete models are utilized to investigate the material behavior by simulating the interaction of macro-particles that comprise granular system (Gonzalez M. and Cuitino A.M., 2012;Hentz S. et al., 2004;Misra A. et al., 2018;Tavarez F.A. and Plesha M.E., 2007;Turco E. et al., 2019). Discrete models consist of identifying the grains in the material, their interacting neighbors, and the behavior of the interaction. For these models one needs to be able to define grains in a given volume of interest, their shapes and sizes, their contacts with their neighbor grains, their surface properties, etc. Moreover, this identification needs to be done at any loading increment, considering the fact that during loading all properties of grains and their contacts will change due to the displacement, deformation, and rearrangement of grains. If the mechano-morphology of a given granular materials can be faithfully reproduced throughout a loading process, discrete models can be used to extract high-fidelity results as depicted in Fig. 1 about the micro-level response as well as the micro-macro correlation for the material (Gonzalez M. and Cuitino A.M., 2012;Misra A. and Poorsolhjouy P., 2015b).
GMA incorporates data from the material microstructure and micromechanics into the continuum model in a computationally affordable statistical manner. The antecedents of the granular micromechanics approach may be traced to second quarter of the 20th century in works related to the estimation of wave velocities in granular beds in which the Hertz solution for normal contact between elastic spheres was used to find equivalent elastic moduli of regularly arranged grain packing (see (Duffy J. and Mindlin R.D., 1957;Hara G., 1935) and for a brief review ). This analysis served as the inspiration for the development of stress-strain relations for random packing of grains including certain statistical representation of the packing geometry (Chang C.S. and Misra A., 1990;Digby P.J., 1981;Misra A. and Chang C.S., 1993;Walton K., 1987). In GMA, the exact evolution of grain trajectories within the material point are not investigated. Instead, the macroscopic behavior of the granular system is derived by averaging the behavior of grain-pair interactions in different directions. In our past works, we have presented a derivation of the GMA based model applicable to cementitious granular materials in the framework of thermo-mechanics (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017). Here, this thermomechanical derivation of GMA is enhanced to account for the effect of grainsizes. The paper is structured as follows: In section 2.1, a brief review of the general derivation of GMA within a thermo-mechanical viewpoint is presented for completeness. Section 2.2 focuses on the specifics of coupled damage-plasticity laws at micro-scale for deriving the correct macroscopic behavior of the material (for a more complete description of the micro-scale laws, the reader is referred to (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017)). In Section 2.3, the relationship between grain sizes on one hand and model variables (including branch vector, density of grain-pair contacts, as well as stiffness and damage variables) on the other, have been discussed and size-effect laws are implemented within the model. Section 3 is devoted to results of simulations for samples of various grains sizes and under different confinements. Particularly, effects of changing grain sizes on failure and its type (brittle or ductile) has been discussed. Finally, a summary and conclusion, along with an outline of continued research is presented in Section 4.

Thermo-mechanical derivation
The GMA hypothesizes that the deformation energy of a volume of granular materials can be effectively expressed in terms of the deformation energies of all its grain-pair interactions. As a consequence, the macro-scale Helmholtz free energy density, W, and dissipation potential, ψ, can be obtained as volume averages of their grain-pair counterparts, which are expressed in terms of grain-scale variables as (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017)   c 1 1 , , where the summations go over all grain-pair interactions, superscript α refers to the α th grain-pair, total number of grain-pair interactions is N c , and the superimposed dots imply time derivative. The grain scale kinematic variables consist of the relative displacement and the plastic relative displacement between neighbor grains, j . Here, we note that in this model, the relative displacements between grains represents the relative displacement between the centroids of each two neighboring grains. This relative displacement can have various sources including deformation of individual grains, or relative movement between grains centroids, or deformation concentrated at the contacts, e.g. in Hertzian contact. These phenomena are not addressed separately in this model and their collective effect is lumped into the relative displacement between grain centroids. The same is true for plastic relative displacement between grain centroids. In addition, the unitless grainscale damage parameter, j  D [53] , is a vector variable which represents the loss of grain-pair stiffness coefficients. In the entire document, the summation convention over repeated indices is implied, unless explicitly noted otherwise.
In this analysis, we only account for the so-called macro-displacements and do not introduce any additional kinematic descriptors or higher-gradient terms and grain-spins/rotations that are necessary for complete description of the complex granular systems (Misra A. and Poorsolhjouy P., 2016a;2016b;Nejadsadeghi N. and Misra A., 2020;Poorsolhjouy P. and Misra A., 2019). In this case, , is related only to the first gradient of macro-displacement field. With this in mind, the macro-scale Helmholtz free energy density, W, and dissipation potential, ψ, can also be defined as functions of macroscopic variables as follows where the superscript p denote plastic, ε ij is the independent kinematic variable (strain tensor), and D ij is the damage tensor which for a given direction, n j , yields the damage vector, D i , i.e. D i = D ij n j . As standard in classical continuum mechanics, the power density of internal forces is expressed as ( , ) , where σ ij is the Cauchy stress tensor. Using the Clausius-Duhem inequality and the above definitions, the following is obtained Since the equality in Eq. 3 should be satisfied for any arbitrary set of kinematic variables, all three terms must simultaneously vanish. For the second and third terms, since the coefficients are functions of the same kinematic variables, the weaker orthogonality condition can be sufficient (Ziegler H., 1983). As a result, the following equations are obtained Further substituting Eq. 1a into Eq. 4a and applying chain rule the following expression for the Cauchy stress can be obtained In Eq. 5 the derivative of grain-pair interaction energy with respect to the displacement is defined as the conjugate force measure in the grain-scale Assuming that the macroscopic variables including strain, ε ij , plastic strain, , and damage tensor, D ij , are distributed linearly in the material point, the corresponding local variables defining the relative kinematics of two neighbor grains, i.e. the relative displacement, δ i , plastic relative displacement, p i  [57] , and the damage vector specifying the loss of stiffness between the two grains, D i , can be derived (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017). Approximating the displacement and plastic displacement of a grain using the Taylor series expansion of the displacement and plastic displacement of its neighbor grain, allows us to define the relative displacement and relative plastic displacement between two neighbor grains (δ i and p ) as the product of their macroscopic counterparts (ε ij and p ij  [58] [59][60]-入力あり ) with the branch vector. In other words, the displacements are written as the projections of their macroscopic counterparts, in the direction of the grain-pair interaction, and then scaled by the magnitude of the branch vector. In addition, based on the definition of the macroscopic damage tensor, D ij , the grain-scale damage vector, D i , is given, by definition, as the product of the damage tensor and the direction of grainpair interaction (i.e. as the projection of tensor, D ij , in the direction of the grain-pair). ; ; where l j is the branch vector joining the centroids of two neighbor grains which is given as the product of the length of the branch vector, l, and the unit vector in the direction of the contact, n j as l j = ln j . Substituting Eq. 6 and 7 into Eq. 5 the expression for Cauchy stress can be rewritten as It is remarked here that the widely used Eq. 8 is valid under specific assumptions as outlined in Eq. 7. Further, we observe that based on the definition of the damage vector and damage tensor (see Eq. 7) the partial derivative of the damage vector with respect to the damage tensor is given as ∂D k /∂D ij = n l ∂D kl /∂D ij = δ ik n j . In the same way, the partial derivative of the rate of the damage vector with respect to the rate of the damage tensor can be derived as ∂Ḋ k /∂Ḋ ij = n l ∂Ḋ kl /∂Ḋ ij = δ ik n j . The same relationships can be formulated for the partial derivative of plastic displacement with respect to plastic strain as, , and for their time rates, . It is therefore derived that . Substituting this assumption, as well as Eqs. 1a and 1b into Eqs. 4b and 4c, results in Now it should be noted that the summations in Eq. 9 should hold for any number of grain-pair interactions and thus, the equality should be satisfied in a term-by-term manner.
Combining Eq. 6 and 10 results in which leads to the following micro-scale Clausius-Duhem type inequality for the α th inter-granular interaction: where the grain-scale dissipation, d α , is defined as and the time derivative of micro-scale free energy, W α , is Thus, employing appropriate formulation for the Helmholtz energy, W α , and dissipation potential, ψ α , we can completely describe the grain-pair force-displacement relationships (using Eq. 12) and in turn, the classical stressstrain behavior of a granular system (using Eq. 8).

Grain-scale coupled damage-plasticity relationships
At this juncture, it is useful to formulate the grain-pair energy and dissipation functions in terms of grain-pair displacement decomposed into two components, one along the vector connecting the two grain-centroids (branch vector), termed as normal component, and the other in the orthogonal direction, termed as the tangential component (as depicted in the inset in Fig. 2). Accordingly, a local Cartesian coordinate system, nst, is defined composed of three unit vectors: unit normal vector n i in the direction of branch vector, l i , and two orthogonal vectors lying in the tangential plane denoted as s i and t i . The direction cosines of these local coordinate axes are given as cos ,sin cos ,sin sin using spherical coordinates (See Fig. 2). The normal and tangential components of a grain-pair are then given as It should be noted that the subscripts n and w, refer to normal and tangential components and they do not represent vector quantities. Summation convention, therefore, does not apply to them. The grain free energy and dissipation potential are written as the summations of their normal and tangential components, . To account for damage and plasticity, the normal and tangential components of inter-granular free energy take the following form (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017)   are the grain-pair plastic displacements in normal and tangential directions, respectively. Normal and tangential components of inter-granular dissipation potential are represented by the sum of the dissipations due to damage and plasticity where dissipation coefficients Y and Z denote generalized forces that are conjugates of the inter-granular damage and plastic dissipation, respectively. The forms assumed in Eqs. 17 and 18 result in rate-independent form of damage and plasticity as opposed to rate-dependent form that are also relevant to certain granular systems (Giorgio I. and Scerrato D., 2017;Misra A. and Singh V., 2014a;2014b). Substituting Eqs. 17 and 18 into Eqs. 10a and 10b rewritten for normal and tangential directions one can find the dissipation coefficients, Y and Z, duals of inter-granular damage and plastic displacements respectively.
Further, normal and tangential components of grain-pair force vector can be found by substituting Eq. 17 into Eq. 7   n n n n n n n The details of particular forms of the damage functions D n and D w and related parameters that are applicable to cohesive granular materials have been described in details in our previous work and is therefore not repeated here . Fig. 2 contains the force-displacement relationships, as well as schematic curves of the functional forms of the force-laws. Moreover, a selection of the major parameters used in the relationships and their role in the force-displacement curves are depicted in the Fig. 2. The tension-compression asymmetry of grain-pair behavior and the post peak softening due to damage under grain-pair tension and shear is clearly illustrated. The grain-pair tension-compression asymmetry is signif-icant in induced anisotropy evolution in granular systems (Jia H. et al., 2017). Table 1 lists all the model parameters and their numerical values used in this study. For more details and implementation of loading-unloading-reloading criteria, ratio between failure force in normal compression to that in normal tension (captured in parameter R), effect of grain-scale relative normal displacement on the grainscale damage parameter in shear direction (captured in parameters α 1 , α 2 , and α 3 ), the effect of confinement on the compressive behavior of grain-pair interactions (captured in the parameter α 4 ), as well as the incorporation of permanent deformation in grain-pair interactions (captured in parameters β n and β w ), the reader is referred to (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017).
Substituting Eqs. 7 and 20 into the equation for Cauchy stress (Eq. 8) and taking the strain tensor out of the summation, the equation for the fourth rank stiffness tensor can be derived as where N c is the total number of grain-pair interactions within the material point. Total number of grain-pair interactions can be calculated as half of the product of the total number of grains, N, and the mean coordination number, z. Using two new parameters, namely the volumetric number density of grain-pair interactions ρ c = N c /V and the directional density distribution of grain-pair contacts, ξ(θ,ϕ), this summation over all contacts can be transformed into an integration over all different orientations as  where l is the average value of grain-pair branch length in the whole RVE. Using the same approach, the Cauchy stress will also be derived in the integral form as For more detailed description of the procedure used to transform the summations into integration see (Misra A. and Poorsolhjouy P., 2016c).

Scaling laws for grain-size effect
Let us consider two different granular materials whose average grain sizes are d 1 and d 2 , where without loss on generality say d 1 < d 2 , and the subscripts refer to the material type 1 and type 2. A set of size scaling relationships are hypothesized, inspired by the principles of renormalization, to compare the two materials based upon physical arguments. For example, it is reasonable to expect that the average grain-pair branch length of material 2 is related to that of the material 1, through the following simple scaling (in the subsequent equations (Eqs. 24-28) the subscripts 1 and 2 will refer to the two granular materials and are not tensor quantities) Similarly, the volumetric number density of grain-pair interactions, ρ c , of the two types of granular materials can be related to each other through the following scaling relationship Where the exponent, α, can vary widely depending not only upon the average grain sizes but also the grain size distributions as well as the packing density. The exponent, α, will also depend upon the problem dimension (that is whether the granular material is a packing of 2D grains or 3D grains) and the grain geometry. For the ideal case of 3D packings of monodisperse spherical particles in a regular closed-pack assembly, the exponent, α, will reach the lower bound of -3 (negative 3). It is evident that the value of the exponent α should be negative since the volumetric number density of grain-pair interactions is inversely related to grain size. Calibrating the value will be complicated for random packings of cohesive particles with irregular geometries and complex particle size distributions. However, the exponent will be bound within the bound -3 < α < 0.
Further, the deformation behavior of grain-pair depends not only upon the grain-shape, its surface characteristics as well as the cohesive mechanisms, but also on the size of the grain. While the grain size dependence of grain-pair stiffness is complex, the following size scaling is expected even for grains of similar composition and shapes         2 2 n n w w 2 1 2 1 1 1 where parameter β determines the relative scaling of grainpair initial (unloaded) stiffness coefficients depending on grain-size. Considering size-effect laws discussed in the literature (Cordero N.M. et al., 2010;Fredrich J.T. et al., 1990), the numerical value of parameter β is expected to be negative, that is the larger grains result in less stiff grain-pair in comparison to smaller grains. In addition, the complete force-displacement behavior of a grain pair, particularly its failure behavior, will also be scaled with the change of grain size. The following scaling relationships are introduced for those model parameters that describe the relative displacement where the grain-pair forces reach their peak Based on experimental results in the literature, it is evident that coarser grains show a more ductile behavior, which means they can store more energy. As a result, the damage parameter used in our model (in normal tension B n , in normal compression, B nc , and in tangential direction, B wo ) should be proportional to the grain size. That is, aggregates of larger sizes will fail at a larger displacement. This implies that the size-effect parameter should have a positive value, γ > 0. Fig. 3 presents an illustrative parametric study to demonstrate how the size-effect parameters influence the grain-pair behavior. The schematic curves in Fig. 3 present force-displacement relationships for normal interaction of grain-pairs, under tension, for varying grain-sizes and for varying size-effect parameters. The effect of the size-effect laws on tangential interaction as well as the normal compressive interactions is similar and not presented here for brevity. It is seen in Fig. 3a that negative β values lead to stiffer behavior for smaller grains. It is also seen that by keeping β constant and increasing γ, smaller grains show smaller energy storage capacities and smaller failure stress. Inversely, keeping β constant and increasing γ results in larger grains showing a more ductile behavior.
It is also observed (from both Figs. 3a and 3b) that for positive values of γ, if γ = -β, smaller grains show larger stiffness, however, they fail at the same force value and smaller displacement values compared to larger grains. Accordingly, larger grains have a more ductile behavior and a larger work of failure (more energy storage capacity) which is in agreement with typical empirical observations reported in the literature.
It is noted that here that the choice of parameter set for simulating the behavior of the specific material investigated in this paper, was done by considering the macroscopic behavior of the material. At the macroscopic level, combining (i) Eq. 20, which describes the mathematical form of the grain-pair stiffness coefficients as function of the stiffness and damage parameters, (ii) Eq. 24, which describes the effect of grain sizes on the magnitude of the branch vector, and (iii) Eq. 26, which describes the effect of grain sizes on the volumetric number density of grainpair interactions and substituting them into Eq. 22 which derives the macroscopic stiffness tensor, we can derive the ratio between initial macroscopic stiffness tensor of a material with grain size of d 2 to that of a material with grain size of d 1 as follows since the magnitude of the branch vector, l, is scaled linearly with size, volumetric number density of grain-pair interactions, ρ c , scales with grain size through exponent α, and the initial stiffness coefficient of interacting grains at the unloaded state is scaled by β. Accordingly, the initial macroscopic stiffness tensor of a material with grain size d 2 can be derived as a function of the stiffness tensor of the material with grain size of d 1 as follows It should be noted that this simple close-form solution is only available for the components of stiffness tensor at the initial (virgin) state. Once the material is subjected to loading, every grain-pair will experience a unique history of relative displacements according to its direction and therefore its force-displacement relationship will be scaled in a nonlinear manner dictated by Eq. 27. As a result, the macroscopic behavior, represented by the evolution of macroscopic tangent stiffness tensor, will be scaled by grain size in a nonlinear fashion which depends on the loading path. In what follows, numerical results of the effects of changing grain size on the stress-strain behavior of a specific granular material is presented to demonstrate the effectiveness of the model and the scaling laws presented.

Results
Inspired by the experimental results presented in (Fredrich J.T. et al., 1990), the enhanced GMA based model described in section 2 was applied to simulate the conventional triaxial behavior of cohesive granular materials for a range of grain sizes. The nominal grain sizes of for these simulations ranged from 0.02 m to 0.00001 m. Furthermore, the confining pressure for these conventional triaxial simulations ranged from 0.0 MPa to 400.0 MPa. In these simulations, the parameter values used in our previous work on cohesive granular materials (Misra A. and Poorsolhjouy P., 2015a;Poorsolhjouy P. and Misra A., 2017) are used (See Table 1).
Calibration of size-effect parameters was done by investigating their combined effect on the grain-scale behavior as well as the macroscopic continuum level behavior. By considering the physical meaning of each of the exponents, as it was already discussed in Section 2.3, it is clear that α < 0, β < 0, and γ > 0. However, the choice of the exact value of the parameters is done based upon their combined effect in the overall behavior of the material and by mimicking experimental results. In this paper, we calibrate model parameters based on the experimental results presented in (Fredrich J.T. et al., 1990). We first observe that, for the material under consideration, the initial stiffness coefficients are inversely related to grain-sizes. Considering this and using Eq. 29, it is clear that for this material we need to confine the size-effect parameters such that 2 + α + β < 0. However, it is also seen in (Fredrich J.T. et al., 1990) that this inverse relationship is not very strong and therefore, |2 + α + β| should be a small number.
Once these the limits stated above are set, the calibration of the exact values of the parameters, is done by comparing the model results to the stress-strain curves for conventional triaxial test under different confining pressures for materials with different grain sizes. Just like in the grain-scale (as shown in the parametric studies in Fig. 3), the effect of size-effect parameters on the macroscopic behavior of the material is also complicated. Accordingly, the choice of the parameters used in this study and for the material under consideration (presented in Table 2) is done for deriving an overall fit between the results. A detailed investigation of the influence of the various parameters remains a topic for further research and investigation.
The simulated conventional triaxial compressive experiments are typically performed in two, sequential phases. In the first phase, termed as the confinement phase, the material specimen is subjected to an incrementally increasing hydrostatic confinement to the desired confining pressure, such that (σ 11 = σ 22 = σ 33 = σ conf ). In the second phase, the material specimen is subjected to deviatoric loading by incrementally increasing the stress in the vertical direction. In this deviatoric phase, the vertical strain, ε 11 , is increased while the two horizontal stresses (σ 22 and σ 33 ) are kept constant.
It is well known for cohesive granular materials that increasing confinement results in increasing the material tangential stiffness as well as the failure stress under triaxial loading. The simulations described in this paper, show that reducing grain-size can also result in increasing stiffness and failure stress. For a clear visualization, the results of the triaxial compression simulations are categorized into two cases: case 1 for small confinements (0.0, 5.0, and 20 MPa), and case 2 for large confinements (50,100,200,and 400 MPa). This distinction has been made to highlight the significant change in the nature of stress-strain relationship as the confinement becomes larger than a value lying between 20 and 50 MPa. The different natures of the stress-strain relationship are clear from the plots given in Figs. 4 and 5 for the two cases. In Fig. 4, we give the evolution of the deviatoric stress, q = σ 11 -σ conf , as the applied vertical strain, ε 11 , is increased for granular materials of different grain sizes for the case of small confinements (< 20 MPa). Further, in Fig. 5, we show the same set of results for granular materials of different grain sizes for the case of large confinement (> 50 MPa). It is clear from Figs. 4 and 5, that both the stiffness and strength increase as the grainsizes decrease. Moreover, this stiffening and strengthening effects are greater at the higher levels of confinements as seen from Fig. 5.
More importantly, the results in Figs. 4 and 5 illustrate the transition in the nature of stress-strain behavior from "brittle failure" to "ductile failure". The observation that failure behavior transforms from brittle to ductile with increasing confinement is well-known from experiments (Bažant Z.P. et al., 1996;Mahboubi A. and Ajorloo A., 2005;Sfer D. et al., 2002) as well as from simulations using GMA that support the experimental results (Misra A. and Poorsolhjouy P., 2015a). A clear peak in the stress-strain relationship, defined by the point at which the stiffness (slope of the stress-strain curve) reaches zero, characterizes brittle failure. On the other hand, the lack of a peak but a clear change in slope of the stress-strain curve characterizes a ductile failure. Granular materials showing ductile response have a yield point (defined here as the stress and strain point at which there is a significant reduction in stiffness), followed by ideal plastic behavior or a hardening-type behavior.
As expected for a given grain size and changing confinement, the simulation results presented in Figs. 4 and 5 confirm the above observations. It is noteworthy, however, that the transition from brittle to ductile failure can also be affected by decreasing the grain size for a given confinement as explicated in both Figs. 4 and 5. Indeed, if we examine the behavior for different grain sizes at 20 MPa confinement given in Fig. 4, we can observe that for grain sizes d > 0.0001 mm, there is a clear peak in the stress-strain curves which implies a brittle failure. We note that in this context, the term "brittle failure" signifies a loss of ability of the material to support additional stress. The described "brittle failure" is distinct from the term "brittle behavior", which, as opposed to "ductile behavior", implies a low capacity to store energy (i.e. a low work of failure). In contrast, for the same confinement (20 MPa), a "ductile failure" can be observed in Fig. 5 characterized by a yield point followed by strain-hardening for the material. The transition between brittle and ductile failure happens for material with smaller grain sizes of d ≤ 0.0001 at a confining pressures between 5 MPa and 20 MPa (see Fig. 4 for d = 0.0001, 0.00002 mm and 0.00001 mm). For larger grains, this transition is observed at confinement between 20 MPa and 50 MPa as seen from comparison of results in

Figs. 4 and 5.
This transition from brittle to "ductile failure" is also observed experimentally as seen in (Fredrich J.T. et al., 1990). The above observation can be further elucidated by investigating the evolution of failure stresses and yield stresses with the grain-sizes. For the case of small confinements, we can identify a failure point as the peak of the stress-strain curve. Fig. 6a gives a plot of the failure deviatoric stress versus grain sizes for small confinements. On the other hand, for large confinements, it is possible to identify a yield deviatoric stress as the point of rapid slope change of the stress-strain curve. In Fig. 6b, the yield deviatoric stress is plotted against grain sizes for large confinements.
For both the failure and yield stress we see a clear trend of inverse relationship with grain sizes from Figs. 6a and 6b. Interestingly, both failure stress and yield stress increase proportionally with increasing 1/ d as a good approximation (note the values of R-squared in the graphs). It is remarkable that the same linear patterns describing the dependency of the failure and yield values on 1/ d [72] has been observed experimentally for granular geomaterials (Fredrich J.T. et al., 1990). The inverse relationship of yield stress with grain-size have been shown through finite element modeling using generalized continuum models (Cordero N.M. et al., 2012) as well as experiments (Jia D. et al., 2003).
Finally, it is interesting to investigate the effect of grain-sizes on the evolution of lateral deformations during triaxial tests. In Fig. 7, the evolution of volumetric strain (ε v = ε ii = ε 11 + ε 22 + ε 33 ) is plotted against the applied vertical strain, ε 11 , during the deviatoric phase. In these plots, Fig. 4 Stress-Strain in Small Confinement: Evolution of the deviatoric stress, q = σ 11 -σ 22 , during triaxial loading (application of ε 11 ) for initial confining pressures of 0.0, 5.0, and 20.0 MPa for all different grain sizes. the compressive volumetric stress is shown as positive as is standard in granular mechanics. These results show that increasing grain size (as well as increasing confinement) result in a transition from dilation to compaction during triaxial testing. Triaxial tests with low confinement experience compaction (compressive volume strain) over a small interval of axial strain before transitioning to dilation (tensile volume strain). For larger confinement, however, the volumetric strain evolution tends towards compressive region. In fact, for very large confinement, no dilation is observed and the volume monotonically decreases as the loading progresses.
More importantly, these results show the clear effect of grain sizes on the volumetric strain evolution. As grain sizes decrease, the volumetric strain evolves increasingly to the tensile region. Indeed for the smallest grain size (d = 0.00001 mm) used in the simulations, the volumetric strain if found to be tensile for all the simulated confinement levels. These results predict that granular materials composed of small grains tend to be more dilative than those composed of large grains. Moreover, the dilative behavior of small-grained materials appears to have limited effect of confinement (particularly for the granular system simulated in this work). We note that in the present study, the focus upon dense cohesive granular materials such as rocks, concrete, and other cementitious materials, in which the effect of initial density variation on the volumetric strain is typically limited. It is noteworthy though that in the present model the effect of initial porosity is incorporated within the parameters of grain-pair interactions. Here, we have evaluated the effect of confining stress on volumetric strain for different grain sizes as shown in

Summary and conclusions
The described work represents an effort to systematically investigate the influences that grain-sizes have upon the macro-scale or collective response of granular systems with very large number (> 10 6 ) of grains. It was shown that   (0,5,20,50,100,200,and 400 MPa). Curves for higher confinements (50 MPa to 400 MPa are too close to distinguish). Note that in these plots compression is taken as positive. the GMA based model for damage and plasticity of cohesive granular materials, which was previously derived and validated by the authors, can be enhanced to incorporate the effects of grain size. To this end, a micromechanical viewpoint is utilized to implement the grain-size effect. The set of scaling laws was incorporated into the GMA based model. These include scaling laws for: volume density of grain-pair interactions, grain-pair branch lengths, grainpair elastic stiffness in normal and tangential directions, and the grain-pair damage parameters in normal and tangential directions. In the presented modeling approach, these microscopic (grain-scale) scaling laws are found to profoundly influence the macro-scale stress-strain response of a granular system.
The derived model is applied to simulate triaxial compressive tests for cohesive granular materials for a range of grain sizes and confinements. The simulation results demonstrate that reduction in grain size results in a stiffer stress-strain behavior, a transition from brittle to ductile type of failure characterized by peak or yield points in the stress-strain curves, and an increase in the failure and yield stresses. These effects, in some ways, mirror the behavior of granular materials for a given grain-size with change in confinement. It is remarkable, therefore, that change of grain-sizes show similar variation in the material macro-scale behavior. Furthermore, the volumetric behavior of the samples are also seen to depend, in addition to the amount of initial confinement, on the grain sizes. Samples of smaller grain sizes demonstrate more tendency towards dilation (tensile volume strain) compared to samples of larger grain sizes. At very small grain sizes, no amount of compaction is sufficient to keep the volumetric strain in compressive domain. Future work will consider further details of the effect of grain-size on deformation and failure mechanism of granular system, including the formation of shear, compression-tension or mixed-model localization bands and their orientation with respect to the loading directions.

Nomenclature
B n , B w = damage parameters in the normal and tangential directions C ijkl = fourth rank stiffness tensor d 1 , d 2 = grain-sizes D i , D n , D w = grain-pair damage vector and normal and tangential components D ij = continuum damage tensor E n = normal stiffness parameter f i , f n , f w = grain-pair force vector and normal and tangential components G w = tangential stiffness parameter k ij , k n , k w = grain-pair stiffness tensor and normal and tangential components l i , l = grain-pair branch vector and branch length n i = unit vector in the direction of grain-pair branch vector N = total number of grain N c = total number of grain-pair interactions R = compression-tension ratio s i = unit vector on the plane orthogonal to grainpair branch vector forming t i = unit vector on the plane orthogonal to grainpair branch vector forming W, W n , W w = macro-scale, normal and tangential grain-pair Helmholtz free energy density Y = generalized force conjugate to grain-pair damage Z = generalized force conjugate to grain-pair plastic displacments α = size scaling parameter for grain-pair volumetric number density α 1 , α 2 , α 3 , α 4 = model parameters controlling damage β = size scaling parameter for grain-pair initial stiffness γ = size scaling parameter for grain-pair damage parameters δ i , δ n , δ w = grain-pair displacement vector and normal and tangential components ε ij = infinitesimal strain tensor ϕ = spherical coordinate σ ij = Cauchy stress tensor θ = spherical coordinate ψ, ψ n , ψ w = macro-scale, normal and tangential grain-pair dissipation potential