Evolution of the Mesoscopic Parameters and Mechanical Properties of Granular Materials upon Loading

1Department of Geotechnical and Underground Engineering, College of Water Resources and Hydropower, State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China 2China Merchants Chongqing Communications Research & Design Institute Co., Ltd., Chongqing 400067, China 3National Engineering Laboratory for Highway Tunnel Construction Technology, Chongqing 400067, China 4State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China


Introduction
Granular materials are common and widely encountered in many fields, which have a close relationship with geotechnical engineering, such as the construction materials of rockfill dams, ground of the highways, and the cohesionless filling materials of pile body for draining purposes.Furthermore, some natural disasters, such as debris flows and landslides, can also behave as the flow of granular materials to some extent.All of these challenging problems force us to further explore the mechanical features of granular materials.Although the meso-and macromechanical properties of granular materials have been studied recently and great achievements in the relationship of stress-force-fabric have been achieved, the linking relationships between meso-and macroscales are still unclear and need further study.The mesoscopic parameters from the concept of complex network, including the coordination number, degree, clustering coefficient, and average shortest path length, can evolve and reflect the geometrical features during the process of loading and are rarely investigated on granular materials, which may establish a linking between the two scales of meso and macro ones and thus help us formulate the stress-strain relationship of granular materials.
It is widely recognized that the internal structure of granular materials [1][2][3][4][5] has a significant influence on their macromechanical properties.Adopting reasonable parameters [6][7][8] and bridging the micro-and macrofields [9][10][11] have attracted a great deal of research in recent years [12][13][14].Although many experiments are performed on the fabric of granular materials [7,15,16], which are complicated and expensive, the analysis is in dispute under some conditions.Initiated by Cundall and Strack [17], the discrete element method (DEM) for modeling the mechanical behavior of granular materials has been used to address many challenges in geotechnical engineering [5,[18][19][20].The dispersibilities of the particle assemblies in granular materials give rise to some distinctive characteristics.Casagrande and Carillo [8] first presented the anisotropy of soil, including the initial and induced anisotropy, and the latter refers to the resulting anisotropy of soil samples during the shear process.Kruyt [21] performed a series of DEM experiments for planar granular materials.By varying the friction coefficients, ranging from almost no friction to heavy roughness between the particles, we obtained the critical state of granular materials and the evolution of the coordination number of specimens during the process of a quasi-static deformation.Oda and Kazama [2][3][4] studied the fabric anisotropy of sands without disturbing the arrangements of the assembled particles and then cut the fabric from the samples and found differences in the initial fabric of soils with a microscope.It was found that the original structure had a profound influence on the macromechanical behaviors of granular assembles.By studying the statistical data of both the fabric and distribution of the contact force, Rothenburg and Kruyt [13] investigated the evolution of the induced anisotropy of granular materials in two-dimensional DEM.Kruyt [21] used DEM to simulate constant pressure tests of anisotropic samples and discussed the variations of the fabric tensor during quasi-static deformation.He concluded that the invariants of fabric tensor are the coordination number and fabric anisotropy.Kuhn [1] emphasized two potential sources of induced anisotropy in granular materials, both the anisotropy of the direction of the contacts and the contact forces.In addition, he focused on analyzing the characteristics of fabric and the material strength in a critical state.Li and Dafalias [22] proposed a new theory regarding the critical state (ACST) of the granular materials, which provided the function of the anisotropy tensor in the critical state.All of the above studies focused on characterizing the fabric tensor in assemblies of particles.In addition, Tordesillas and Muthuswamy [23] investigated the evolution of both the fabric and contact force in dense granular materials from the perspective of complex networks.
It is clear that the mesoscopic structure has a significant influence on the mechanical properties of granular materials, but its evolving rules and the relationships between mesoscopic parameters and macromechanical features are not clear and need further study, which will be carried out here in the perspective of complex network.In this paper, we mainly concentrate on the influences of mesoscopic parameters on the macroscopic characteristics of granular material according to DEM simulation.First, we simulated the macroscopic mechanical properties of granular material under different lateral pressures.Then, the features of granular materials of the mesoscopic structure were analyzed based on the concept of a complex network, which consequently relate the mesostructure to the macromechanical properties of the granular materials during the loading process, especially at critical state.

Micromechanics
In the paper, we only discuss micromechanics in 2D, which can easily be extended to 3D.Granular materials consist of dispersed particles that interact at contacts [21].For the two particles  and , as shown in Figure 1, the definition of the branch vector    is the vector connecting the centers of two particles that are in contact.The corresponding unit normal vector at the contact is    .Many researchers have defined the stress tensor in granular materials [13,24,25].For example, Drescher and de Josselin de Jong [26] first gave one definition based on the continuum-mechanical approach.They supposed the spherical assembly of volume  consisting of grains that have an arbitrary shape.The assembly is subjected to external forces,   1 ,   2 , . . .,    , on its boundary points,   1 ,   2 , . . .,    .As a result, the average stress of the collection under the loads is defined as follows: The distribution of contact orientations in the assembly of particles can be characterized by the symmetric fabric tensor   [27], which is defined as where   is the number of particles in the assembly and  is one contact in the set of contacts .According to the definition of the coordination number , the average number of contacts per particle,  can be expressed as Because       = 1, from (2) in two-dimensional conditions and (3) we can have Therefore, the first invariant of the fabric tensor is  =  1 + 3 .
In biaxial tests, the major and minor principal stresses of the stress tensor are  1 and  3 , respectively.The major and minor principal strains of the strain tensor are denoted by  1 and  3 , respectively, while the major and minor principal values of the fabric tensor are  1 and  3 , respectively.For the principal strains  1 and  3 , the invariants of the strain tensor are the volumetric strain  V and the deviator of the strain tensor   , which are defined as ( Similarly, the first invariant of the fabric tensor is the coordination number , and the deviator of the fabric tensor is the contact anisotropy .Based on the principal values  1 and  3 , the invariants of the fabric tensor in two-dimensional conditions can be defined by

Description of Numerical Tests.
For the samples simulated in biaxial tests, their size is 100 mm in height and 40 mm in width.The assemblies composed of discrete particles and discs have a uniform distribution from  min = 0.25 mm to  max = 0.50 mm, as shown in Figure 2, for a sample with porosity of 0.12 and 7967 particles.The stiffness of the lateral wall is set as one-tenth of the grain.The interparticle friction coefficient is set as 0.3 in the loading process.The sample is generated by the method proposed by Cundall and Strack [17], and many researchers employed this method when performing DEM simulation [10,12,20], during which the initial anisotropy within the sample perhaps exists and can be neglected compared with the subsequent loading process.Following the loading path of triaxial compression test, after the sample with initial porosity is generated, the isotropic stresses along vertical (axial) and horizontal (lateral) directions are applied; then the axial load is applied with constant lateral pressure until the sample fails.In addition, the lateral pressure ( 3 ), deviatoric stress ( =  1 −  3 ), axial strain, and volume change in the whole loading process are recorded to gain the macromechanical behaviors of the numerical samples.Sixteen samples of granular materials generated were tested.According to the parameters in [1,9,17], the parameters and arrangement of the tests used here are shown in Tables 1 and 2, in which both normal stiffness and tangential stiffness are 50 MPa/m, particle density is 2630 kg/m 3 , and frictional coefficient is 0.3.

Macromechanical Properties.
Figure 3 shows the deviatoric stress-axial strain curves of the samples at the same  lateral pressure with different initial porosities, in which the normalized curves are also provided.At the same lateral pressure, the samples with larger initial porosities behave strain hardening and those with smaller initial porosities behave strain softening, and all of them reach a similar value at failure or critical state.Figure 4 presents the volumetric strain-axial strain curves of the samples at the same lateral pressure with different initial porosities.It reveals that the samples with a larger initial porosity contract all the time, while the samples with a smaller initial porosity contract within small axial strain, followed by dilatancy, reaching a stable value at critical state for the different lateral pressure.In addition, the denser the samples, the greater the dilatancy.The simulated results of these numerical tests are consistent with the laboratory experiments in a reference on soil mechanics [28], which shows that the numerical model used here is reasonable for the given parameters and loading conditions.Figure 5 presents the evolution of the void ratio with the axial strain under different lateral pressures in the process of shear loading, which indicates that the void ratios of the samples with larger initial porosities decrease with the increase of the axial strain, while the void ratios of samples with smaller initial porosities reduce first and then increase to a stable value during the shear process.For the same lateral pressure, the void ratio of different samples eventually tends n = 0.12 n = 0.17   to be similar, indicating that all samples will achieve the same critical state with different initial porosities.As shown in Figure 6, the critical void ratio of the samples decreases with increasing mean stress of ( 1 +  3 )/2, indicating that the greater the mean stress, the denser the samples at critical state.

Coordination Number.
Figure 7 presents the relationship curves of coordination number  versus axial strain for 100 kPa, 200 kPa, 300 kPa, and 400 kPa lateral pressures, respectively.At different lateral pressures, the coordinate number  increases gradually and reaches a stable value at critical state for initial loose samples with larger initial porosity, but it has bigger value within small axial strain and decreases rapidly to a stable value at critical state for initial dense samples with smaller initial porosity.The reason for this phenomenon is that, for the initial loose samples, the particles will move close and have the tendency to contact more neighboring particles in the process of loading; and for initial dense samples, the particles are close to each other within small axial strain accompanying volumetric compaction, followed by dilatancy which makes the particles move apart and tend to contact less neighboring particles in the subsequent loading process.With the continuous increase of axial strain, the coordination number with different initial porosities tends to be identical in the critical state for the same lateral pressure, which indicates that the samples with different initial porosities will reach the critical state [11,21,29,30].
To explore the evolution of the coordination number with lateral pressure, two representative samples with initial porosities of 0.17 and 0.30 in Figure 8 are shown.For dense samples with small initial porosity, the coordination number under all lateral pressures sharply increases in the initial state, then gradually decreases, and eventually reaches a stable value.For loose samples with large initial porosity, however, the coordination number under all lateral pressures increases gradually with the axial strain until it reaches a steady state.For both initial dense and loose samples at critical state, the greater the lateral pressure, the larger the value of the coordination number.

The Evolving Characteristics of Mesoscopic Parameters
4.1.The Degree of the Particles.The definition of the degree of a particle for the samples simulated is given here, which is a concept from graph theory [30].According to the definition in graph theory, if the vertex V  is one endpoint of some n = 0.12 n = 0.17 edge   , then V  and   are said to be incident with each other.It is assumed that the vertex V is one point of the graph .The degree (V) of V is the number of edges of  incident with V, counting each self-loop twice.The basic unit in granular materials consists of separate particles, which contact each other.Therefore, each particle is considered to be one vertex in the numerical samples, and its surrounding contacts can also be regarded as the edge incident with this vertex.As shown in Figure 9, the degrees of the central particles are 4, 5, and 6 for (a), (b), and (c), respectively.
From the above analysis, we can confirm that the critical state is very significant to study the mesostructure of granular materials.Hence, the distribution of the degree of particles during the critical state will be investigated here.Figure 10 presents the distribution of the number of particles in the n = 0.12 n = 0.17   samples versus their average degrees for samples with different initial porosities in the critical state at a lateral pressure of 400 kPa.For all of the samples shown in Figure 10, the distribution of the number of particles with varying degrees is similar and has the largest degree of five.The number of particles relative to the degree of five is roughly distributed symmetrically, and with the reduction in the initial porosity, the maximal values of the number of particles gradually increase.This denotes that, at critical state, the particles with degrees of five, four, and six are the main patterns to bear the external loads, which can also be obtained from Figure 11.[31] is a typical property of complex network.In terms of a generic graph , this is defined as the ratio between the triangles and connected triples of the vertices (triads) in , as shown in Figure 12, that is, the fraction of connected triples nodes, which also form triangles, which is expressed as From the definition of the clustering coefficient, we can determine that it means the aggregation degree of the node in the graph.This concept is introduced here to describe the overall contact network properties in the sample.It is observable that the contact in the material is divided into two types, particle-particle and particle-boundary types.To more reasonably discuss the internal contact network, we remove the contacts between the particles and boundaries in the overall samples.In Figure 13, we show the evolution of the clustering coefficient of different compacted granular materials under four lateral pressures.The evolving trend of the clustering coefficient is almost the same as that of the coordination number  in Figure 7.When the lateral pressure is lower, the variation in the curves has no obvious rules.With increasing lateral pressure, especially for 400 kPa, we can clearly see that the smaller the porosity, the greater the clustering coefficient until the samples reach the critical state.In addition, the samples with different initial porosities tend to have identical values of clustering coefficients when the axial strain reaches a certain value and the samples are at critical state.In Figure 14, we present the evolution of the clustering coefficient with the axial strain under different lateral pressures.When the granular materials are dense or have smaller initial porosities, the clustering coefficient suddenly increases at a small axial strain; then, it decreases with increasing axial strain and gradually reaches a stable value.When the granular materials are loose or have larger initial porosities, the clustering coefficient increases gradually with axial strain and then remains stable.The evolving laws of these curves are similar to the stress-strain curves according to macroscopic mechanical behaviors of the samples.

The Clustering Coefficient. The clustering coefficient
In Figure 15, the clustering coefficient of samples with various initial porosities increases at critical state with increasing lateral pressures, indicating that the larger lateral pressures could make the internal particle contact more compact and consequently the samples have a higher strength.However, the clustering coefficient of the samples with different initial porosities is almost identical at critical state under the same lateral pressure, indicating that the internal structures of samples with various initial porosities tend to be similar at critical state.

The Average Shortest Path
Length.The shortest path [31] plays an important role in a complex network.If we want to confirm that some particles are located in a certain force chain, the shortest path could provide a simple method.As a result, it is effective in taking advantage of this concept to describe the internal structure of the samples.In the graph,   is the length of the geodesic (the minimum number of edges) from node  to .The maximum value of   is called the diameter of the graph, which is indicated as  in Figure 16 for samples at a critical state.It is clear that the diameter of the contact network has almost no change with varying lateral pressure, and the denser samples have larger diameters than looser samples in the critical state at the same lateral pressure.The average shortest paths, also known as the characteristic path length, are defined as the mean value of geodesic lengths over all node couples: Based on Figure 17, the average shortest paths decrease with lateral pressure and tend to be stable in the critical state because the lateral pressure could make the granular particles contact closely so that average paths between particles become small.The average shortest paths increase with the decreasing initial porosities of the samples, which indicates that the internal contact network of the samples with smaller porosities will be denser in the critical state.Therefore, we can conclude that lateral pressure has an important influence on the average shortest paths and little effect on the diameter of the internal network in granular materials.

Conclusions
Some of the concepts introduced from the complex network were employed here to investigate the internal structure of granular materials and the mechanical properties at the mesoscale under biaxial compressional conditions.The following conclusions can be drawn: (1) The evolution of the coordination number is obviously different between dense and loose samples.The coordination number of dense assemblies has a sudden increase at small axial strain and then decreases rapidly with the increase of axial strain until it reaches a stable value, while that of loose samples increases gradually with the increase of axial strain and then reaches a stable value.However, the coordination numbers tend to be identical under the same lateral pressure at critical state for the samples with different initial porosities.
(2) The distribution of the degree of particles is very similar for the samples with different initial porosities.The number of particles with a degree of five has the largest distribution, and the distribution of the remaining particles is roughly symmetrical.In addition, the distribution of the degrees can satisfactorily reflect the macromechanical behavior of the samples during the loading process.
(3) The clustering coefficient of simulated samples is closely related to the coordination number and reflects the mechanical behavior of granular materials.In the critical state, the clustering coefficient of all samples is almost identical under the same lateral pressure, indicating that the internal structures of the different assemblies tend to be similar.
(4) The diameter of the contact network of the granular material can be changed with lateral pressure in the critical state.The average shortest paths of the samples with the same initial porosity decrease with increasing lateral pressure and tend to be stable in the critical state.As a result, we conclude

Figure 1 :
Figure 1: Two particles interact at contact.The branch vector    is the vector connecting the centers of two particles that are in contact.The corresponding unit normal vector at the contact is    .

Figure 3 :Figure 4 :
Figure 3: The stress-strain curve of samples under different lateral pressures: (a) lateral pressure of 100 kPa; (b) lateral pressure of 200 kPa; (c) lateral pressure of 300 kPa; (d) lateral pressure of 400 kPa; (e) the normalized curves at the lateral pressure of 100 kPa; (f) the normalized curves at the lateral pressure of 200 kPa; (g) the normalized curves at the lateral pressure of 300 kPa; and (h) the normalized curves at the lateral pressure of 400 kPa.

Figure 5 :Figure 6 :
Figure 5: The relationship between the void ratio of samples and axial strain under different lateral pressures: (a) lateral pressure of 100 kPa; (b) lateral pressure of 200 kPa; (c) lateral pressure of 300 kPa; and (d) lateral pressure of 400 kPa.

Figure 8 :
Figure 8:  The relationship between the coordination number and axial strain of samples under different pressures: (a) dense state ( = 0.17) and (b) loose state ( = 0.30).

Figure 9 :Figure 10 :
Figure 9: The concept of degree: (a) degree of four; (b) degree of five; and (c) degree of six.

Figure 11 :
Figure 11:  The relationships between the number of particles with different degrees and initial porosity in the critical state at 400 kPa of lateral pressure.

30 Figure 15 :Figure 16 :
Figure 15:  The relationship between the clustering coefficient of granular material and lateral pressure with variable porosities.

Figure 17 :
Figure 17: The average shortest paths of the contact network of granular material in the critical state.

Table 2 :
The arrangement of the tests.