On the relevance of the mesoscale in an upscaling approach for granular materials

This paper deals with the interest of upscaling from a meso-scale understand and model the mechanical behavior of granular materials. The meso-domains associated to the meso-scale are defined as closed loops of particles in contact. Six sets of meso-domains with different local textures called "phases" have been defined. The analysis of the behavior of these phases through 2D DEM simulations of a compression test followed by an extension test has allowed the discriminant evolution of the structure and then the role of the different phases at the meso-scale to be evidenced. Constitutive modeling of granular materials can advantageously be developed from the meso-scale where the principles of continuum mechanics apply.


Introduction
So far, a straightforward and reliable modeling of granular materials remains an open issue.The difficulty lies in the very discrete nature of these materials composed of grains in contact.In this study, we explored a new way to perform an uspcaling approach considering a meso-scale, defined at the level of clusters of particles [1,2].Following [3,4], six sets of meso-domains with different local textures called "phases" has been defined.An analysis of the behavior of these phases through 2D DEM simulations is provided.It tends to emphasize the key role of the meso-scale in an upscaling approach to model the mechanical behavior of granular materials.
The numerical biaxial compression test is achieved using the commercial software PFC2D.The sample is composed of 25, 000 circular particles, whose diameter is uniformly distributed from 4.6 mm to 9.2 mm.The normal and tangential stiffnesses at contact are equal to 10 9 N/m, the intergranular friction ratio is equal to 0.57, the walls are rigid and frictionless, the numerical damping are taken equal to 0.3, and the simulations are gravity-free.The compaction process leads to an isotropic stress state of 100 kPa with a porosity of the sample equal to 0.167.From this isotropic state, the biaxial test consists of a loading stage followed by an unloading stage (Figure 5 -gray lines).Herein, the compression direction 1 ′ is such that: where Ėd 1 is the diagonal component of the incremental deviatoric strain tensor Ėd .

Definitions for the meso-scale
We propose to define an intermediate scale at which both the local strain and stress fields can be consistently defined and the local texture can be consistently described [3,5,6].
To this end, we rely upon methods proposed in the literature for partitioning a granular medium into local structures.It must be noted that a significant number of particles in the sample do not have any contact or have only one or two contacts because of the gravity-free assumption.These particles, called rattlers, do not participate in supporting the external loading.They are thus disregarded for all the following analysis.
The sample is first paved with triangles connecting the centers of mass of three neighboring particles using a Delaunay triangulation [7].Each edge of this triangulation is then considered either as a closed branch if connecting two particles in contact, or as an open branch if connecting two neighboring particles without contact.Triangles sharing a common open branch are then associated to obtain local polygonal structures, here called the meso-domains (Figure 1).The meso-scale is then defined as the scale of these local structures.Upscaling from basic Delaunay cells cannot provide a correct upscaling because they do not hold a physical meaning as a local pore around which a mechanical transfer of forces is made possible.As a consequence, they cannot capture the local change of volume.
The structure at the meso-scale is characterized by two meso-variables: the meso-density and the meso-texture.The meso-density is usually defined through the porosity, which is the ratio of the void volume inside the mesodomain m, V v m , to its total volume, V m .For each meso domain m, the meso-texture is characterized by the loop tensor L m , from which two scalars are defined [3,8]: its orientation, θ m , and its elongation degree, β m : where r m is the valence of the meso-domain m, l k is the length of the closed branch k and n k its outward unit vector, D m is the deviatoric part of the loop tensor L m and m m its the major principal direction, and 1 ′ is the compression direction.To analyze the role of the meso-domains texture throughout a loading path, the sample is divided into sets of meso-domains sharing the same texture characteristics.This sets, here called the phases, are obtained considering two sorts of classes: • 3 orientation classes: the class 1 ′ of meso-domains nearing the compression direction 1 ′ , the class 2 ′ of meso-domains perpendicular to the compression direction, and the class 0 ′ of meso-domains oriented in an intermediate direction; • 2 elongation classes: the class W of meso-domains weakly elongated (β m < β m V , where • V is the volume-weighted operator) and the class S of mesodomains strongly elongated (β m > β m V ).The six phases are then denoted by: S 1 ′ , S 0 ′ , S 2 ′ , W1 ′ , W0 ′ , W2 ′ and the the mean orientation of a given phase p is denoted by θ p .The relevance of these phases to characterize the meso-scale and to be a framework for a classical volume average based change of scale are justified through the analyse of their spatial distribution which showed no spatial localization of any phase [4].
The definition of the strain tensor increment of a mesodomain m, δε m , is based on the definition of the incremental gradient tensor for a triangular element e proposed by [9,10] and used in [3][4][5][6]11]: where s ∈ {1, 2, 3} is the sides index of the cell e, δu s is the difference of the incremental displacements between the two vertices of a side s, ∆x s is the vector defined by  side s, and I is the identity matrix.The incremental displacement gradient tensor for a meso-domain m, δa m , is then calculated as the volume average of δa e and the strain tensor increment of a meso-domain, δε m , is obtained as the symmetric part of δa m .The definition used for the stress tensor in a mesodomain m is the one proposed by [5] and used in [4,6,10,11], based on the external contact forces acting on the particles of the given meso-domain: where K m ext is the set of external contacts for the mesodomain m, x c is the position vector of contact c, f c is the external contact force acting at contact c, and Vm is the volume delimited by the boundary on which the contact forces f c applied.Then, the following assumption is done: σ m in volume V m is equal to σm defined in volume Vm .The validity of this assumption was proved in [5].
For a given phase p, the strain tensor increment, δε p , and the stress tensor, σ p , are obtained as the volume average of δε m and σ m over all the meso-domains of the considered phase, respectively: 3 Evolution of the meso-structure of the loading path, which corresponds to an isotropic state of stress for the sample, one can note that the meso-domains are not strictly isotropically oriented due to boundary effects.During loading, the meso-domains are mainly oriented in the major principal compression direction.Then, during unloading, the major principal compression direction changes from the vertical direction to the horizontal direction, and the meso-domains are subjected to a reorientation according to this new major principal compression direction.Similar results can be observed when considering the two elongation classes W and S.
To quantify the evolution of this anisotropy, the texture tensor of a given phase defined in [4] is used: where V p is the volume of the phase p, N m the number of meso-domains in phase p, m m is defined by Equation (2).The second invariant of the texture tensor is used to quantify the intensity of anisotropy.Whatever the loading direction, the phases of the orientation class 1 ′ develop a high level of anisotropy T p d , while the phases of the orientation class 2 ′ have a weak contribution to overall anisotropy (Figure 3(a)).On the contrary, the anisotropy of the phases 0 ′ hardly evolves and remains stable throughout loading and unloading.During the loading, new mesodomains are created and others are destroyed.This leads to an evolution of the volume of the six phases during the loading-unloading test.
In terms of percentage by volume of the different phases, the initial state is homogenous.Throughout loading, the volume of the meso-domains of orientation class 1 ′ increases, reaching a maximum and then decreases slightly, while the volume of the meso-domains of orientation class 2 ′ decreases quickly before stabilizing (Figure 3(b)).The volume of the meso-domains of orientation class 0 ′ does not evolve significantly.In fact, new contacts are preferentially created in the compression direction 1 ′ , which generates many meso-domains in this direction and increases their total volume.Throughout unloading, the same observations hold if we consider the new compression direction which has changed (see on Figure 2).Whatever the loading direction, the granular material mainly bases its mechanical resistance on the meso-domains of orientation class 1 ′ .

Analysis of meso-stress/strain
The analysis of the behavior of the different phases is performed by means of the evolution of the stress ratio and of the incremental strain ratio respectively defined by: where σ p i ′ is the stress component of phase p in the principal direction i ′ of the macroscopic stress and where δε p i ′ is a component of the incremental strain tensor for phase p in the principal macroscopic stress direction i ′ .A positive (respectively negative) value of δε p i ′ corresponds to a dilating (respectively contracting) phase.Figure 4(a-b) show the evolution of these ratios.For class S phases, the initial stress state in the six phases is not necessarily isotropic and is very different from one to another.The six stress ratios R p σ reach very different plateaus with values greater for phases 1 ′ and smaller for phases 2 ′ .Furthermore, these values are independent of the considered stress path (corresponding to the critical state).The behavior of the phases is closely related to their orientation with respect to the loading path as well as to the distance between the current anisotropy level and the one at the critical state for the considered loading path.In Figure 4(b), at the beginning of the loading stage, the incremental strain ratio for class S take very different values for the different phase, then these values increase and reach more or less stabilized values.This ratio remains positive (dilation) for phases 1 ′ and negative (contractancy) for phases 2 ′ .The conclusions that can be drawn for the unloading stage are similar to those obtained in the loading stage but with the adequate definition of orientation 1 ′ .These conclusions are consistent with

Upscaling process
The upscaling process from the meso-scale to the macroscale is based on the previously defined phases.The consistency of the mechanical behavior of these sets with respect to the loading direction has been here-before demonstrated.The percentage by volume of these sets evolves accordingly.So it is possible to define the evolution of the stress ratio and the volumetric strain at the sample level by considering the volumetric means of these two variables: These two volumetric means are compared to the values obtained from the boundary conditions (Figure 5).The very good agreement between these two results implies that: • the meso-scale is an appropriated local scale at which the local mechanical behavior can be defined; • the proposed meso-variables defined in the framework of continuum mechanics are relevant static and kinematic variables, even though the definition of the mesostress raises some intrinsic issues; • given the behavior of each phase, the knowledge of the percentage by volume occupied by each of them is a sufficient information to perform the change of scale.

Conclusion
A DEM simulation of a biaxial loading and unloading has been carried out in this work.Six sets of meso-domains called "phases" have been defined.The mechanical behavior of all the considered phases was found clearly dependent on their orientations with respect to the loading direction and on their current anisotropy level.Furthermore, the volumetric amount of phases oriented in the compression direction (respectively extension direction) increases (respectively decreases) during the loading.The very good agreement between the mechanical behavior defined from the sample boundary and the one defined from the volumetric average of the mechanical behavior of the set of phases proves that the complex behavior of granular material at a sample scale can be understood by an upscale from the meso-scale.

Figure 1 .
Figure 1.(a) Set of polygonal meso-domains and (b) Triangular elements, with open and closed branches, defined by the Delaunay triangulation and meso-domain (bold lines) obtained by association of the triangles sharing an open branch.

Figure 2 .
Figure 2. Polar distributions of the number of meso-domains for the reference states.Loading path: (a) initial state, (b) peak state.Unloading path: (c) initial state, and (d) critical state.

Figure 2 1 Figure 3 .
Figure2depicts the evolution of the distribution of orientations of all meso-domains for some usual reference states: the initial state, the maximum shearing resistance state (peak state), and the critical state.At the initial state

1 Figure 4 .
Figure 4. (a) Stress ratio and (b) Strain ratio for phases of elongation class S and for different orientation classes.Blue lines: 1 ′ ; Red lines: 2 ′ ; Green lines: 0 ′ ; Black lines: sample.Solid lines: loading path; Dashed lines: unloading path.