Role of inter-particle friction in granular materials under three dimensional conditions

Role of Inter-particle Friction in Three Dimensional Mechanics of Granular MaterialsS. Yang, W. Zhou, J. Liu, T. Qi and G. Ma VI International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2019 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, & Gang. Ma Role of inter-particle friction in granular materials under three dimensional conditions –– PARTICLES 2019 Shuhan. Yang¹, Wei. Zhou¹, Jiaying. Liu¹, Tianqi. Qi ¹, Gang. Ma¹ ¹ State Key Laboratory of Water Resources and Hydropower Engineering Science Wuhan University 8 Donghu South Road, 430072 Wuhan, China E-mail: [email protected](Shuhan Yang); [email protected] (Wei Zhou) Key words: DEM, Inter-particle Friction, True Triaxial Tests, Strong and Weak Contact System, Coordination Number, Anisotropy. Abstract. The inter-particle friction is known to be an important contributor to the strength and deformation characteristics in granular materials. The mechanism of inter-particle friction to the macroscopic responses can be explained by microscopic investigations. Based on the discrete element method (DEM), a series of true triaxial tests for the cubic granular assembly are carried out and the effects of inter-particle friction coefficient (μ) on the evolutions of macro- and micromechanical parameters of granular materials are studied. The macroscopic stress, the distribution of coordination numbers and contact force with regard to strong and weak contact networks are concerned, as well as the corresponding fabric tensor and anisotropies. Findings indicate that increasing inter-particle friction sharpens the peak value of deviatoric stress and enhances the degree of dilatancy of the granular assembly at the macroscopic level. From the microscopic perspective, the distribution of the coordination number of the weak contact system varies dramatically, while the number of particles with smaller coordination number in the strong contact system changes little with different μ. Besides, the difference between strong and weak contact networks is enlarged, and anisotropy indicators are significantly enhanced, which strengthen the bearing ability of anisotropic stresses in granular materials. 1 INTRODUCTION Granular materials are closely linked to our daily life. Many important infrastructures concerning the quality and safety of our everyday life, such as rockfill dams, dikes and foundations, are built with geotechnical friction-dissipative granular materials. The performance of granular materials under external loads directly affects the design, construction and operation of these structures [1]. Under shear, granular materials exhibit extremely complex mechanical behaviors at the macroscopic level while the typical multiscale features contribute a lot to the overall complexity [2]. The granular materials are composed of discrete particles whose macroscopic responses are vitally interrelated to their interaction. Inter-particle friction as one of the key factors affecting the interplay between particles, may significantly contribute to the macro- and micro-responses of granular materials. DEM (Discrete Element Method) has demonstrated its ability in reproducing the 388 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 2 macroscopic responses and exploring the microscopic mechanism of granular materials [3]. Based on DEM simulations of granular assemblies, contact orientation [4, 5], force transmission [6, 7], contact networks [14] are regarded as the important characterization of the microstructure. Besides, the fabric tensor [9, 10, 11, 12] has been defined to measure the intensity and orientation of the contact texture anisotropy. The interplay between particles is the intrinsic factor that causes the macroscopic responses of granular materials. To figure out the effects of inter-particle friction on granular materials, Rothenburg et al. [13] observed the relationship between void ratio and coordination number along different inter-particle friction by biaxial tests. Huang et al. [14] evaluated the sensitivity of critical state behaviors to inter-particle friction. Maya et al. [15] identified that increasing inter-particle friction promotes the formation of straighter chains and a greater degree of branching in the force chain network by two-dimensional simulations. Antoy et al. [16] investigated how particle-scale friction affect the mechanism of mobilization of macroscopic shear strength by conventional triaxial tests. However, true triaxial tests are rarely used in these simulations, which can reflect the real situation in nature or engineering. The main purpose of this paper is to investigate the effects of inter-particle friction on the characteristic behaviors of granular materials under true triaxial tests. The stress–strain relationships and peak state of deviatoric stress are presented. From the view of strong-weak contact network, the micro-responses are explored, including the coordination number, the contact force, the fabric tensors and the anisotropies. The contribution of micro-investigations to macro-mechanical properties is further discussed. 2 DEM SIMULATIONS OF TRUE TRIAXIAL TESTS 2.1 Stress-strain invariants The stress invariants describing true triaxial stress state are generalized shear stress q, average hydrostatic pressure p and Lode angle θσ: 123 3 p    （1） 2 22 12 23 31 1( )( )( ) 2 q        （2） 3 tan 2 b b    , 23 13 b     （3） The intermediate principal stress ratio b responses the relative magnitude of the intermediate principal stress σ2 between the minor principal stress σ3 and the major principal stress σ1. The range of b is 0≤b≤1. The Lode angle θσ defines the angle between the major principal stress and the deviatoric part of stress in the principal stress space. The true triaxial test with b = 0.0 represents the triaxial compression test (σ2 = σ3, θσ=0) and b = 1.0 indicates the traxial extension test (σ1 = σ2, θσ = π/3). 2.2 Sample preparation and loading paths The numerical sample is composed of 31 253 non-contacting spherical particles with Gaussian distribution in a 400mm×400mm×400mm cube. The minimum and maximum 389 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 3 particle sizes are 2 and 8.5 mm, respectively. Fig.1 shows the particle size distribution (PSD). In order to generate isotropic samples in the initial fabrics, we use the displacement control to compress the samples with a constant-speed in all directions until the target size. After that, isotropic consolidation is carried out and three-dimensional isobaric stress is applied to the specimen until the desired confining pressure 0.5 MPa is reached. Finally, the initial void ratio of the specimen is 0.582. The Hertz–Mindlin contact model is adopted. The inter-particle friction coefficient μ during loading process is set as 0.05, 0.1, 0.2, 0.3, 0.4 and 0.5, respectively, and other numerical parameters are shown in Table 1. Fig.1 Radius distribution Fig.2 Stress path under true triaxial tests Table1: Microscopic parameters of numerical simulation Parameters Value Density ρ/(kg/m3) 2600 Particle radius d/mm 2-8.5 Young’s modulus E/GPa 65.0 Poisson’s ratio ν 0.4 coefficient of restitution e 0.95 sliding friction coefficient μ 0.05,0.1,0.2,0.3,0.4,0.5 In this paper, a series of true triaxial tests are under the constant-p and constant-b loading condition, and the stress paths are shown in Fig.2. In order to eliminate the boundary influence and avoid the occurrence of strain localization, these tests are carried out with periodic boundary. The stress in each direction of the specimen is controlled by the migration rate of periodic boundary. The intermediate principal stress σ2 and minor principal stress σ3 can be estimated from Eq. (4) and Eq. (5) respectively: 1 21 3 (1 ) (1 ) (2 ) pb bb b       （4） 1 3 3 (1 ) (2 ) pb b     （5） 390 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 4 3 MACRO-SCALE RESPONSES For various inter-particle friction, the evolution curves of generalized shear stress q and volumetric strain εv along deviatoric strain εd are shown in Fig. 3 (the positive volumetric strain represents contraction, negative value represents dilation). Increase in inter-particle friction μ steadily raises the shear strength of a granular assembly (the peak deviatoric stress) while the peak value occurs at a similar strain level (about 6%) irrespective of various μ. The deviatoric stress changes to a downward trend after reaching the peak value, which means strain softening. The larger inter-particle friction can enhance the strain softening behavior. In addition, diverse μ values lead to different volume responses of a granular assembly. Specifically, after a small compression deformation, the specimens are all in the dilatation condition. The particle assembly with large inter-particle friction starts the dilatation state earlier and shows a more evident dilatation degree at later stage of loading. Fig.3 Macroscopic response of granular assemblies with different inter-particle friction: evolution of deviatoric stress (a) and volumetric strain (b) versus deviatoric strain εd Fig.4 Peak state of deviatoric stress along different friction coefficients 391 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 5 Corresponding to different inter-particle friction, the three-dimensional stress surfcae in π plane is plotted in Fig. 4. For each given μ, the curve presents similar shape like triangular cone but differs from one another in size. The stress evolution in π plane reflects the peak deviatoric stress of the granular assmebly at different b and μ values. It can be seen that, the peak deviatoric stress decreases with the increase of b but improves with the increase of μ. Moreover, the sensitivity of the specimen to inter-particle friction debases along with μ, for the smaller increaseing extent of deviator stress. The stress characteristics shown in Fig. 4 conform to the general rule of true triaxial test and are also consistent with the conclusions of other numerical tests [16, 17, 18]. 4 COORDINATION NUMBER The geometric stability of a granular assembly under mechanical loading is generally studied in the matter of its apparent coordination number Z (i.e., average number contacts of per particle) at a given stage of loading. As a mesoscopic scalar index of granular materials, average coordination number Zt is given by 2Nc/Np where Nc and Np are the total number of contacts and particles in the specimen respectively. It can also reflect the overall volume change of a granular assembly. In generally, higher average coordination number is corresponding to smaller porosity, closer interaction between particles, and more stable/isotropic internal structure. On the contrary, the degree of anisotropy increases. Radjai et al. [5] divided the whole contact network into two complementary sub-contact nectworks, a strong contact subnetwork and a weak contact network that carry normal contact forces larger and less than the average. The two contact networks have different mechanisms for the geometric stability of granular materials. In presenting the results, we have hereby separated the contributions of the weak contact network from the strong contact network, along with the total value. The strong contact network is denoted by Γstrong, the part from the weak contact network by Γweak, and the entire contact network by Γtotal. Fig. 5 shows difference performance among the three contact networks in terms of the relationship between the average coordination number and the inter-particle friction in the triaxial compression condition. Under the given loading path, each contact network presents a law that Zt decreases with the increase of μ, which confirms the macroscopic phenomenon that the dilatation degree of the specimen is more obvious with large μ. Besides, average coordination number in Γweak is larger than that of Γstrong, reflecting stronger internal stability and higher isotropy of Γweak. Considering the relationship between the average coordination number and the volume deformation of granular materials, it indirectly shows that the weak contact network contributes greater to the deformation of a granular assembly. In a granular media, the coordination number Zt indicates an average over particles within contact networks while the number of contacts neighbors n varies from particle to particle. Bratberg et al. [19] recognized that the connectivity of a granular assembly should be given by the fraction Pn, defined by particles having n contact neighbors. From the study of Liu et al. [20], particles with fewer coordination numbers are in lower degree of restraint condition, corresponding to higher anisotropy. Conversely, particles are more stable and show isotropy feature. From the distribution of coordination number at stress peak state in Fig. 6, the coordination number Z in accordance with the abscissa denotes Pn, and the corresponding ordinate refers to 392 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 6 the number of n adjacent particles in contact with a single particle. In the total contact network (Fig.6a), the maximum particle content changes from P5 to P4, meanwhile, all distribution curves seem to move left entirely with the increase of μ value. Evidently, some particles with n≥5 contacts lose neighbor contacts, causing the number of particles with n≤4 contacts increase. The distribution of coordination number within weak contact network (Fig.6b) is similar to that of the total contact network, while the maximum particle content moves from P4 to P3. From that, the anisotropy feature of Γweak is stronger than that of Γtotal and both of them diminish owing to the increase of μ. It is interesting to take a quantitative look at the strong contact network (Fig.6c) whose distribution of coordination number is quite different from that of the total contact network. Thanks to the increase of μ value, the number of particles (Pn≥P3) decreases. However, μ value appears to have little effect on the number of particles whose coordination number is 1 and 2, whereas P1 and P2 constitute the end and main body of the force chain in the contact network. It shows that the number of particles forming the force chain accounts for 1/3 of the total number and does not fluctuate with the change of μ value, as if satisfying the most advantageous condition of particulate distribution. This demonstrates that the shear resistance of granular assembly can be mobilized by selecting relatively few benificial orientated contacts to transmit the greater than average contact force. Compared with strong, weak and total contact network (Fig.6d), the three curves can be roughly considered as the left translation of the total contact network curve. For the total contact network, the weak contact network makes more contribution to the coordination number distribution than the strong contact network. The influence of inter-particle friction on coordination number distribution of Γweak is greater than that of Γstrong. The weak contact network tends to be isotropic, balancing hydrostatic pressure of the granular media, dominating overall volume deformation; while the strong contact network has obvious anisotropic feature, and bear most of the deviatoric stress of the granular assembly (shown in next chapter). Fig.5 The relationship between the average coordination number and the inter-particle friction 393 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 7 Fig.6 The distribution of coordination number at stress peak state: Number distribution of coordination number Z in the total contact network (a) the weak contact network (b) the strong contact network (c) and the comparison (d) while μ=0.2 5 FABRIC TENSOR AND ANISTROPY 5.1 Contact force Radjai et al. [6] pointed out that the strong and weak contact network in a granular media have different mechanical mechanisms. Strong contacts have a decisive influence on the contribution of macro-mechanics. Fig. 7 shows the evolution of the average normal contact force and tangential contact force in the strong, weak and total contact network during true triaxial loading. The average normal contact force is much larger than the tangential contact force, indicating it’s leading position. There is great disparity between the average normal contact force of strong and weak contact network, and the numerical relationship between them is about 4 times. For the total contact network, the evolution contribution of the average normal contact force is dominated by the strong contact network. With the increase of μ value, the normal contact force and tangential contact force in each contact network increase 394 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 8 obviously, while the average coordination number of them decreases. This reasonably explains the increase in the mechanic anisotropy feature of the granular assmebly (see Fig. 8), and highlights the supporting effect of the strong contact network on balancing external load. Fig.7 The relationship between contact force and inter-particle friction (a) Normal contact force (b) Tangential contact force 5.2 Definition of Contact Fabric Tensor Based on the work of Yimsiri et al. [21] in quantitatively analyzing the anisotropy of a particle assembly, two classes of anisotropy with different mechanisms are differentiated: geometric anisotropy and mechanical anisotropy. The geometrical anisotropy is defined by the local orientation of contact plane which can be expressed by the distribution of contact normals and branch vectors. The mechanical anisotropy, however, is mainly created by the external forces and can be splited into normal force anisotropy (caused by normal contact forces) and tangential force anisotropy (induced by tangential contact forces). This paper adopts the following expression of fabric tensor proposed by Satake [22] and Oda [23]: () ij i j E nnd      （6） where n is the unit vector along the normal direction of the contact surface. Θ characterizes the orientation of n relative to the global coordination system. E(Θ) is the distribution probability density function, which can be expressed in terms of the second-order Fourier expansion: 1 ( ) (1 a ) 4 c ij i j E nn    （7） where the second-order anisotropy tensor a c ij is deviatoric and symmetric, which characterizes the geometrical anisotropy induced by normal contact. ' a 15/ 2 c ij ij   （8） where ' ij  is the deviatoric part of the fabric tensor ij  . 395 Shuhan. Yang, Wei. Zhou, Jiaying. Liu, Tianqi. Qi, and Gang. Ma 9 Only when the particle system is composed of non-spherical particles, the branch vector contributes greater to geometric anisotropy. The distribution of branch vectors can be expressed in the way similar to formula (7) and formula (8): 1() 4 ij i j d d nnd      （9a） 0 ( ) (1 a ) l ij i j d d nn  （9b） where 0 ' a =15/2 / l ij ij dd（ ） is the contribution value of branch vectors to the geometric anisotropic tensor. 0 = ii dd is the average length of branch vectors in the domain Θ and ' ij d is the deviatoric part of ij d . The distribution function of normal contact force and tangential contact force are expressed respectively as (10b) and (11b). Correspondingly, the normal contact force tensor and the tangential contact force tensor are defined by equation (10a) and (11a). 1() 4 ij n n ij f nnd      （10a） 0 ( ) (1 a ) nn ij i j f f nn  （10b） and 1() 4 ij t t ij f tnd      （11a） 0 () ( ) ttt ik k kl k l i i f f an annn     （11b） wherein ' 0 15 2 ij n ij n a f   ， ' 0 15 2 ij t ij t a f   0n ii f   is the average normal force in the domain Θ. The second-order anisotropy tensors ,,, ij ij ij ij cl nt aaaa are deviatoric and symmetric, then, the degree of anisotropy can be quantified by using the second invariants of the three anisotropy tensors defined above: ** * 3 2 ij ij a aa （12） where the sub/super-script * stands for contact normal, c, branch vector, l, normal contact force, n, and tangential contact force, t, respectively. 5.3 Characteristics and Evolution of Anisotropy The characteristics of anisotropy and its evolution during the loading process are depicted in Fig.8. Owing to the use of spherical particles with a relatively narrow size distribution, the branch vector is almost in the same direction as contact normal vector and contributes so little to the total strength that can be ignored. The anisotropic coefficient of the initial specimen increases rapidly with the increase of the deviation strain, and the increase rate becomes more sharply with larger inter-particle friction. It can be seen evidently that the peak value of the normal contact force anisotropy an (about 3% of the deviatoric strain) comes earlier than that of the contact normal anisotropy ac (about 6% of the deviatoric strain), and the peak value is higher than the latter. With the increase of μ value, the difference is more obvious. Interestingly, the contact tangential force anisotropy at reaches its peak within 1% of the 396