Undrained behavior of binary granular mixtures with different fines contents

Accepted Manuscript Undrained behavior of binary granular mixtures with different fines contents

Wei Zhou, Wei Wu, Gang Ma, Tang-tat Ng, Xiaolin Chang PII: DOI: Reference:

S0032-5910(18)30749-6 doi:10.1016/j.powtec.2018.09.022 PTEC 13696

To appear in:

Powder Technology

Received date: Revised date: Accepted date:

22 January 2018 4 August 2018 10 September 2018

Please cite this article as: Wei Zhou, Wei Wu, Gang Ma, Tang-tat Ng, Xiaolin Chang , Undrained behavior of binary granular mixtures with different fines contents. Ptec (2018), doi:10.1016/j.powtec.2018.09.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Undrained behavior of binary granular mixtures with different fines contents Wei Zhou, Wei Wu, Gang Ma*, Tang-tat Ng, Xiaolin Chang 1. Wei Zhou. Professor, State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, 430072, Wuhan, China. E-mail: [email protected]

PT

2. Wei Wu. Ph.D. candidate, State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, 430072,Wuhan, China. E-mail: [email protected]

RI

3. Gang Ma. Associate Professor, State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, 430072,Wuhan, China. (corresponding author)E-mail: [email protected]

SC

4. Tang-Tat Ng. Professor, Department of Civil Engineering, University of New Mexico, Albuquerque, NM 87131. E-mail: [email protected]

MA

NU

5. Chang Xiao-lin. Professor, State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, 430072, Wuhan, China. E-mail: [email protected]

AC C

EP T

ED

*Corresponding author. Tel.:+86 27 6877 3778; E-mail address: [email protected] (Gang Ma)

ACCEPTED MANUSCRIPT Abstract Binary granular mixtures are comprehensively explored from both macro- and microscopic perspectives. In this study, the undrained behavior of binary granular mixtures constituted by fine and coarse particles is investigated using the discrete element method (DEM) under different fines

PT

contents (FCs). The peak deviatoric stress ratio has a linear relation with the effective void ratio. The

RI

increase in the FC to a certain value (approximately 10% in this work) can change the material

SC

response from strain hardening to limited “liquefaction”. The enhanced strain hardening behaviors were observed when the FC exceeds 15%. The location of the critical state line (CSL) depends on

NU

the fines content. Voronoi tessellation was adopted to explore the mesoscale structures of binary

MA

mixtures. The local isotropy, local anisotropy, and sphericity of the Vorono i cells are clearly related to the FC. The radial distribution functions (RDFs) also exhibit clear correlations with the FC. The

ED

abundance and mechanical contributions of different types of contacts, e.g., coarse-coarse,

EP T

coarse- fine, and fine- fine contacts, are analyzed. For these binary systems, the proportions of the sliding contacts with respect to the FCs are contact-type dependent. The coaxiality between the loading direction and fabric orientation of the different contact types depends on the relative

AC C

proportions of the contact types. The evolutions of the global anisotropies are also explored, and the contact-based anisotropy parameters extracted at peak states under various FCs are compared to the global anisotropy parameters. The stress- force-fabric relationship is evaluated in various subnetworks of the mixtures and discussed in terms of the strong-weak partition and nonsliding-sliding partition networks. The dominant role of geometrical anisotropy derived from the strong subnetwork is emphasized, and the linear relationship between the global stress ratio and the geometrical anisotropy within strong- nonsliding subnetwork is evidently observed under various

ACCEPTED MANUSCRIPT conditions of fines content.

AC C

EP T

ED

MA

NU

SC

RI

PT

Keywords: binary mixtures; fines content; DEM; anisotropy; coaxiality; Voronoi

ACCEPTED MANUSCRIPT

1. Introduction In geotechnical engineering, the particles in natural sand-silt mixtures have been traditionally classified into two categories: coarse and fine particles [1]. Numerous laboratory tests have been

PT

conducted to explore the behavior of binary mixtures, and attention has been paid to the mechanical

RI

effects of the coarse particles and fine particles on the mixture instability, shear strength, particle

SC

breakage as well as the steady state and liquefaction conditions [2-8]. The fines content (FC) in granular mixtures is the most direct and simplest factor with which to investigate the responses of the

NU

mixtures under various loading conditions. The basic mechanical properties of binary granular

MA

materials are significantly influenced by their fines content [9-13]. In addition to laboratory tests, several discrete element method (DEM) simulations were performed to explore the characteristics of

ED

granular mixtures [14, 15]. Most of the previous numerical simulations utilized the direct shear and

EP T

conventional triaxial compression loading procedures. Xiao et al. [16] conducted a series of drained triaxial tests to analyze the effects of fines content on the stress-dilatancy relationship and the peak and critical-state friction angles. A constitutive

AC C

model was adopted by Xiao et al. [17] to characterize the fines-dependent behaviors of granular mixtures, and the critical state lines (CSL) was found to get elevated with decreasing fines content. Lade and Yamamuro [18] discovered that the fines content contributes to a unique particle structure that exhibits tremendous effects on the instability phenomenon with different extents of liquefaction resistance. Meanwhile, different studies showed contradictory conclusions on the effects of fines on liquefaction potential. Some work [19, 20] discovered positive influences, but other work [21, 22] identified negative influences. Liquefaction of uniform- grade soils is usually achieved with lower

ACCEPTED MANUSCRIPT densities at low or high confining pressures [23]. For gap-graded coarse-fine mixtures, samples can even liquefy at high densities due to the distinct mechanical contributions of the fines and coarse material. Hence, it is necessary to adopt a three-dimensional undrained loading condition to investigate the softening and hardening behaviors of binary mixtures with different fines content.

PT

Although some previous work has concentrated on the particle-scale characteristics of these

RI

mixtures [24-26], few studies have undertaken undrained DEM research on the unique geometrical

SC

and mechanical characteristics of fine particles. Lopez et al. [24] and Minh et al. [25] utilized the force transmission within mixtures to quantify the roles of the coarse and fine fractions. They found

NU

that three contact types (coarse-coarse, coarse- fine, and fine- fine contacts) occupied different

MA

contributions to the deviatoric stress. Gong and Liu [26] analyzed the effects of the differences in the friction among various contact-types during drained biaxial tests and provided the relationship

ED

between the peak-state anisotropy and the percentage by weight of the coarse material. Ng et al. [27]

EP T

conducted simulations with ellipsoidal particles and stressed the important roles of the particle shape and fines content on the peak friction angle and critical state void ratio. Rahman and Lo [28] adopted several experiments to analyze the undrained behavior of fines- in-sand mixtures with limited

emphasized.

AC C

conditions of fines content, and the effects of fines on the compressibility and critical state are

However, few studies have concentrated on the meso-scale properties of binary mixtures, and the scope of fines content are not large enough in some previous work. To further understand the effect of external stresses on mixtures with different fines contents (FCs), we employed several mesoscale and microscale statistical quantities to investigate the local pattern arrangement and contact force transmission. At the mesoscale, the morphological method, in terms of the Voronoi

ACCEPTED MANUSCRIPT polyhedron cell, is employed. The morphological characteristics of the Voronoi cells and the local isotropy and anisotropy parameters are both considered regarding the influence of the change in the fines content, providing a comprehensive understanding of the local structural properties of the particle assemblies. The radial distribution functions (RDFs) are also obtained to characterize the

PT

pair correlations with different FCs. Unlike the conventional anisotropies of the global contact

RI

network, we adopt a contact-type-based geometrical anisotropy and mechanical anisotropy to explore

SC

the unique structures and contributions of the three contact types. Two partitioned methods, strong-weak networks and nonsliding-sliding networks, are also adopted to investigate the

NU

characteristics of anisotropies and the effects of contact types within different fines contents. The

MA

force transmission and sliding fraction of different contact types are also explained. The tendency toward coaxiality [29], representing the divergence of the fabric and loading direction, is explored

ED

under different contact types and various FC conditions during shearing. On the basis of the

EP T

four-partitioned method (strong- nonsliding, strong-sliding, weak- nonsliding, and weak-sliding subnetworks), the relationship between the geometrical anisotropy that derives from the

conditions.

AC C

strong-nonsliding subnetwork and the global stress ratio is investigated under different FC

2. Simulation preparation and theoretical background

2.1 Sample generation and properties Allowing for the fact that increasing the number of fine particles requires a larger sample size which can increase the computational costs [26, 30], the open source code LIGGGHTS [31], standing for LAMMPS improved for general granular and granular heat transfer simulations, is used

ACCEPTED MANUSCRIPT to perform DEM simulations in this study. We use a periodic boundary to guarantee a RVE condition, and the rolling resistance is also incorporated into the simulation to numerically consider the mechanical effects of particle shape. Instead of using a pure binary mixture, we use a gap- graded particle size distribution to model natural mixture systems. The particle size distributions are shown

PT

in Fig. 1(a). The gap ratio, defined as the ratio of the minimum diameter of coarse particles to the

RI

maximum diameter of fine particles, is set to 4, which is much larger than 2.4 used in the study by

SC

Kumara et al. [32]. The diameters of the fine particles range from 1.3 to 1.5 mm. The diameters of the coarse particles vary from 6 to 9.8 mm. The particle size ratio (PSR), defined as the ratio of the

NU

mean diameter of coarse particles to the mean diameter of fine particles, is 5.65 in all of the mixtures.

MA

From the perspective of crystallography, the particle diameter required to occupy the tetrahedral voids within the coarse particles of radius d should be smaller than 0.225d [33]. Therefore, PSR is

ED

large enough for the fines to fit into the voids created by coarse particles.

EP T

The coarse particles play a primary role within the mixtures, especially when the fines content is relatively low [30]. However, many earlier studies changed the number of the coarse particles when increasing the fines content [24-26], thus mixing the effects of the assembly size and fines content

AC C

during loading. The number of coarse particles remains constant when the FC is less than 50% and decreases when the FC exceeds 50%. The assembly sizes of the mixtures with FCs less than 50% are approximately the same. When the FC is sufficiently high, i.e., greater than 50%, the coarse particles are completely surrounded by fines. Therefore, a relatively small representative volume element can be extracted to simulate the material behavior. One typical granular mixture of FC = 25% is shown in Fig. 1(b). Detailed information of the granular mixtures, including the number of each particle type, assembly sizes, and initial effective coordination numbers, are listed in Table 1.

ACCEPTED MANUSCRIPT The numerical samples are prepared by randomly generating particles with reduced sizes in a prescribed domain. Then, the particles are expanded slowly to the target particle sizes [34]. The friction coefficient for the sample preparation is set to 0.35 to obtain moderately dense samples. When this process terminates, the samples are relaxed and then isotropically compressed to a

PT

pressure of 200 kPa. Because the gravity may cause segregation in the granular mixtures, the

RI

gravitational constant is set to zero, which can also help to avoid the erosion phenomenon [27]. The

SC

Hertz-Mindlin contact model and periodic boundary are adopted. The rolling friction is also used to consider the interlocking effect [35]. The input parameters, such as Young’s modulus Y, Poisson’s

NU

ratio  , and the coefficient of restitution re , are summarized in Table 2 [36]. To guarantee a

I  z d

 p  10

-3

[37], where

d

MA

quasi-static loading condition, the axial strain rate is selected to ensure the resulting inertia index is the mean particle diameter, 

is the particle density

EP T

2.2 Void ratios of the mixtures

ED

given in Table 2, and p represents the isotropic compression pressure.

The void-based classification of binary mixtures was first discussed by Vallejo [38], and followed by many other researchers [12, 15, 39]. The minimum void ratio of binary mixture

m ix

e m in

AC C

and the corresponding fines content FCmin can be approximately predicted by the void ratios of the pure fine and pure coarse particles: pc

e m in 

e e

m ix

e

F C m in 

pc

 1  e e

e

pf

pc

pf

 1

(1)

pc

 e

pf

1

where e p c and e p f represent the void ratios of the pure coarse and fine particle systems, respectively. In this study, e p c and e p f

are 0.654 and 0.603, respectively. Therefore, the

theoretical minimum void ratio of the granular mixtures

m ix

e m in

is approximately 0.18, and the

ACCEPTED MANUSCRIPT corresponding FCmin is about 30%. The inter-coarse and inter-fine void ratios e c and e f can be taken as the first-order indices of the effective coarse and fine contacts [1]. They can be calculated by the void ratio of the granular mixture e m ix and FC: c

e

f



e

m ix

 FC

PT

e 

1  FC e

m ix

RI

FC

(2)

SC

The variations in the global, inter-coarse, and inter-fine void ratios with FC are presented in Fig.

NU

2, in which a typical “V” shape is observed in the e m ix versus FC curve. Instead of a sharp transition in the evolution of e m ix , a smooth transition zone is observed in most experimental tests

MA

and numerical simulations [40]. As shown in Fig. 2, the prediction of Eq. (1) is less than the experimental data, especially around FCmin . Both e c and e f research their maximum values at

ED

FCmin . A critical FC threshold commonly employed in previous literature is used to quantitatively

EP T

characterize the FC-based zones [1, 3, 41]. When FC is lower than FCth , the coarse particles play a primary role in controlling the shear strength; when FC becomes greater than FCth , the shear strength

AC C

is simultaneously controlled by both the fines and coarse particles. The threshold FCth remains relatively constant in different mixture systems, i.e., FCth≈20~30% [25, 42-44]. The FCth can be determined from the intersection of the FCth line ( F C

 e

m ix

pf

e m ax

, where

pf

e m ax

represents the

maximum void ratio of the pure fine assembly) and the e m ix versus FC curve [24]. The FCth can be directly calculated from the e m ix versus FC curve using the assumption that

pf

e m ax

equals 1.0 [45].

In this study, FCth is about 31%, which is close to FCmin . Hence, FCth can also be used as an effective indicator of the transitional zone. Fig. 3 illustrates the typical local structures of mixtures with different FC values. When the FC

ACCEPTED MANUSCRIPT is low, the fine particles disperse among the voids between the coarse particles. It can be observed that the floating fine particles can effectively fill the voids generated by the coarse-particle skeleton in Fig. 3(a). Hence, the mixture can be regarded as a coarse-supported structure (Fig. 3(a)). When the FC increases, both the coarse and fine particles affect the sustaining structure, corresponding to the

PT

transitional zones (Figs. 3(b) and (c)). In Figs. 3(b) and (c), the fines can evidently separate the

RI

coarse particles, thus elevating the ability of fines to sustain the whole structures. In the

SC

coarse-supported transitional zone (Fig. 3(b)), the fines effectively fill the voids. However, in the fine-supported transitional zone, the fines can overfill the voids and begin to form a weaker particle

NU

skeleton (Fig. 3(c)). As shown in Fig. 3(d), the packing structure is dominated by the fine particles

MA

when the FC is sufficiently high.

Although the coarse particles play a dominant role in the particle structure when FC is lower

ED

than FCth , the effects of the fine particles cannot be neglected [1]. Intuitively, an assembly with fines c

EP T

is stronger than that without fines, even if the assemblies have the same inter-coarse void ratio e

[30]. The global, inter-coarse and inter- fine void ratios are not sufficient to characterize the mechanical behavior of a mixture. Binary mixtures can be regarded as inter-granular structures with e eq ,

AC C

an equivalent void ratio

which has different expressions for different FC ranges [46]. The shear

responses of mixtures when FCFCth at void ratio

e eq

e eq .

e eq .

e eq

can be similar to that of the pure coarse assembly with

Likewise, the mechanical behavior of binary granular mixtures when

can be expected to be analogous to that of the pure fine assembly with the same

The equivalent void ratio of a binary mixture can be calculated as follows [30, 46]:

ACCEPTED MANUSCRIPT e   FC  1  FC  m ix eeq =  e  1  FC  FC  m Rd  m ix

F C  F C th

(3)

F C  F C th

PT

where  and m are the FC-based material coefficients. It can be understood that the fine particles can fill the voids that are generated by the coarse skeleton or separate the coarse contacts when FC is

RI

increased. The increment of voids in granular mixtures,  V v , can be defined as  V f where

V

f

SC

represents the volume of fine particles [46]. The incremental change of void volume can characterize

NU

the effects of fines on the whole volume behavior. The term    1 means that the fine particles can totally fill the voids within the coarse skeleton;   0 .0 indicates that all the fines particles

MA

participate in supporting the coarse particle structures. R d is the ratio of the mean coarse diameter Dc50 to the mean fine diameter Df50 . Eq. (3) indicates that when the FC is 0.0% and 100%, e e q is

ED

equal to the void ratio of the pure coarse and pure fine assemblies, respectively. The coefficients 

EP T

and m can be obtained from the minimum void ratios: e m in  e m in  1  F C m ix

AC C

 

pc

e m in

and

pf

e m in



,

FC

F C  F C th

e   ln  mp fin  F C   e m in  m ix

ln  1  F C

m 

where

pc



ln R d

(4) ,

F C  F C th

represent the minimum void ratios of the pure coarse and fine particles. To

obtain the minimum void ratios for the various FCs, we set the friction coefficient to zero and then consolidate the samples to designed pressure (200kPa). The equivalent void ratios for various FCs are listed in Table 1. According to Eq. (4) when FC
e eq

is equivalent to the inter-coarse void ratio when 

to -1. When FC keeps increasing towards 100%, the

e eq

is close

gradually approaches the inter- fine void

ACCEPTED MANUSCRIPT ratio, which corresponds to the configuration in Fig. 3(d). 2.3 Mechanical definitions for various subnetworks The average stress tensor within the granular assembly can be expressed as follows [47-49]:

fi

c

V



c

fi d

c

(5)

j

c N

and d cj represent the contact force and branch vector connecting the centers of two

contacting particles, respectively.

N

is the total number of contacts, and V is the volume of the

RI

where

1

PT

 ij 

SC

particle assembly. If the global contact network is partitioned into several subsets, the product of the

where

N

k

1  c1   fi d V  c N 1

c1 j



 c N

c2

fi d

c2 j



2



 c N

ck

fi d k

ck j

  

(6)

MA

 ij =

NU

force and branch vectors can be performed within several subnetworks:

represents the contact number in the kth subnetwork.

ED

The definitions of four anisotropies have been given in many previous studies [50-52]. The c ,k d ,k anisotropies of a subnetwork, e.g., the fabric anisotropy tensor a ij , branch anisotropy tensor a ij ,

EP T

n ,k t ,k and normal and tangential contact force anisotropy tensors a ij and a ij , are defined in a similar

way to the global contact systems [53-56]. The symbol with superscript k indicates a value derived

AC C

from the kth contact subnetwork, and the term without the superscript k is derived from the global contact system. Once the geometrical and mechanical anisotropy tensors are obtained, we can calculate the anisotropic measures of each subnetwork. a

* ,k



3 2

*,k

*,k

a ij a ij

(7)

where * strands for c, d, n, and t, corresponding to the four types of anisotropy. Radjai et al. [57] classified the contact network into a two-partition system, i.e., the strong-weak contact system. Specifically, the contacts with a normal contact force greater than the mean force are

ACCEPTED MANUSCRIPT strong contacts, and the rest of the contacts are weak contacts. Aloso-Marroquin et al. [58] noted that the sliding contacts are very important to the plastic behavior of granular material. Partitioning the contact network into a nonsliding-sliding contact system is another way to understand the binary mixture behavior. Therefore, the kth network mentioned above can be either of the two networks. We

PT

depict the two partitioned methods in Fig. 4. It shows that the sliding contacts are scattered across the

SC

sustain greater contact forces within the granular assembly.

RI

assembly. Generally, compared to the weak and sliding contacts, the strong and non-sliding contacts

NU

3. Macroscopic behavior

MA

The particle assemblies with different FCs are sheared under constant volume axisymmetric triaxial loading conditions. This loading path corresponds to the undrained conventional triaxial test

ED

[59]. The deviatoric stress q and mean effective pressure p are calculated from Eq. (5), and the stress

EP T

paths are given in Fig. 5(a). When the FC is below a certain value (10% in this case), the deviatoric stresses at both the peak state and phase transformation (PT) points decrease as the FC increases, coinciding with the experimental observations of Lade and Yamamuro [18] and Yin et al. [30]. As

AC C

shown in Fig. 5, the limited “liquefaction” phenomenon is observed in the sample with FC=10%. According to Fig. 5(b), it indicates that a strain hardening behavior is evident when the FC is greater than 10%.

The relationship between the peak deviatoric stress ratio q/p and the FC resembles an overturned “V” shape, as shown in Fig. 6(a). The decrease of shear strength when FC increases from 30% to 100% demonstrate the weakening effect of FC in shear strength mobilization. The phenomenon that the global void ratio, which characterizes the packing density, increases at high

ACCEPTED MANUSCRIPT FCs, shown in Fig. 2, can help to explain this strength-decreasing trend at high FCs. Meanwhile, it is interesting to note that the q/p is linearly related to the equivalent void ratio at various FCs, which is illustrated in Fig. 6(b). The critical state theory has been interpreted as the basic framework to establish the constitutive

PT

model [29]. The fundamental concept in the critical state theory is the critical state line (CSL) that

RI

can describe the relationship between the void ratio and mean effective pressure. Hence it is

SC

meaningful to explore the effects of fines content on the critical sta te characteristics. Considering the fact that the steady-state line within undrained tests is equivalent to the critical state line in this paper

NU

[60], a series of drained triaxial compression tests with constant  3 have been conducted under

MA

different conditions of fines contents and confining pressure (0.5 MPa, 1.0Mpa, 2.0 MPa, 4.0 MPa, 7.0MPa, 10 MPa, and 16 MPa). The CSL observed in the triaxial tests can be expressed as [61]:

e0

 e0   cs

 p     pa 



corresponds to the limiting value at

EP T

where the coefficient

ED

e

m ix

(8) p  0 .0

.

e0

also represents the

location of the CSL, and  c s is donated as the slope of the CSLs. As shown in Fig. 7, the CSLs are e0

and  c s are depicted in the

AC C

exhibited under different conditions of fines contents. The terms

inserted figure of Figure 7. Although the CSLs are apparently parallel to each other, the term  c s firstly decreases with increasing FC, and then experiences an ascending trend. It demonstrates that when FC increases the CSLs firstly exhibit an anticlockwise rotation, and then present a clock-wise rotation. The location of the CSL, characterized by

e0

, firstly decreases with increasing FC from 0.0%

to 30%, and then gets increased when FC exceeds 30%. This observation is consistent with the results of previous work [17, 28]. In general, the larger particles sustain the greater contact forces. At low FCs, the coarse-coarse

ACCEPTED MANUSCRIPT contacts dominate. When the FC changes, the dominant contact type and the local structures also alter. Therefore, an in-depth analysis of the meso- and microscale structural characteristics of the coarse and fine particles is required to deepen the understanding of the effects of fines content.

PT

4. Voronoi structure and radial distribution function (RDF)

RI

4.1 Voronoi tessellation of the granular mixtures

SC

The Voronoi tessellation algorithm and radical plane polyhedron are employed to characterize the local structure of binary mixture and elucidate the topological characteristics of the particle- void

NU

system [62-65]. A representative Voronoi cell is shown in Fig. 8(a), in which n is the normal vector

MA

of the outer face. We can define several measurements to characterize the Voronoi cells. To quantitatively analyze the shape characteristics of the Voronoi cells, we adopt the sphericity derived

ED

from the ratio of the surface area related to a sphere of the equivalent volume of the Voronoi cell and

EP T

the surface of the Voronoi cell [66]:

3

R 

2

3 6 π V c e ll

(9)

S c e ll

AC C

where S c e ll and V c e ll are the cell surface area and cell volume. To better investigate the solid parts of the Voronoi cells, the local volume fraction is taken into account and is defined as

  V p a r tic le V c e ll

. The local volume fraction can be used to describe the

isotropic properties of the local pore structures [67, 68]. The Minkowski functionals W 1 and tensor 0 ,2

W1

were used to explore the morphology of the Voronoi-based structures [69, 70]: W1  0 ,2

W1

1 3 



K

1 3

dA

(10)



K

n  ndA

ACCEPTED MANUSCRIPT The tensor with unity trace,

0 ,2

W1

W1

, is adopted to characterize the local anisotropic properties

[71]:

where



is the local anisotropy index and

  1 .0

 m ax

(11)

 m in

indicates an isotropic cell.  m a x and  m i n

are the maximum and minimum eigenvalues of the tensor

0 ,2

W1

W1

.

PT

 

RI

To explore how does the local structure vary with FC, we estimate the probability density

SC

functions (PDF) of different measurements of the Voronoi cells and their statistics. Figure 8(b) shows

NU

such a PDF related to the cell sphericity for increasing FC. The peak points of PDFs gradually rise up when FC keeps increasing, and then approach a steady level when FC are larger than 30%. This

MA

result demonstrates that the extent of the sphericity related to the Voronoi cells becomes more prominent as the fines content increases. When FC is located in the fine-dominant zone, the

ED

magnitude of sphericity nearly keeps constant when FC continuously increases, which can be

EP T

observed from the change of the value of sphericity when the peak points are achieved. Meanwhile, it exhibits second smaller peak points of various PDFs marked by the dashed rectangle in Fig. 8(b).

AC C

The behaviors of the smaller peaks are amplified in the inserted figure. It shows that these smaller peak points gradually decrease and nearly disappear when FC exceeds 30%. The emergence of the second peak can be attributed to the binary components whose sizes are quite distinctive. The PDFs of the local anisotropy index



are depicted in Fig. 8(c). It is clearly shown that the

distribution narrows with increasing FC. The local anisotropy value corresponding to the peak of PDF gradually decreases with increasing FC and approaches a constant value. The PDFs nearly coincide with each other when FC is larger than 30%. It can be explained that with increasing FC the void ratio reaches maximum at FC=30%, thus making the arrangement of particles well-distributed.

ACCEPTED MANUSCRIPT For FCs below 30%, a smaller peak can also be detected when



is located between 1.05 and 1.15,

as shown in the inset of figure 8(c). The smaller peak values also decrease when FC increases. The PDFs related to local volume fraction are plotted in Fig. 8(d). The PDF shifts to right with increasing FC. Similar to the PDFs of sphericity, the  values where the peak points are obtained

PT

gradually increase towards a constant value. The peak points of PDFs of  decrease first, followed

RI

by a slight increase until a stable distribution is reached. The second peak values are evidently

SC

observed when FC is below 30% and decrease with increasing FC, as shown in the inset of Figure 8(d). Therefore, the increase of FC can render the distributions of local volume fraction similar to

NU

each other. It can also be concluded that the local structures can keep stable when FC approaches the

MA

fine-dominant zone (FC>30%).

4.2 Structural characteristics of the binary mixtures

ED

The radial distribution function (RDF), also called the pair correlation function, is adopted to

EP T

quantify the spatial structural characteristics of granular mixtures with different FCs. In statistical mechanics, the RDF in a system of particles describes how the density varies as a function of the distance from a reference particle. Hence it is significant to understand the local interdependency of

AC C

particle assembly by means of RDF. In three dimensions, the RDF is determined by calculating the distance between all pairs of particles and binning the results into a histogram, which can be expressed as [72]:

g r 

n r 2

where n  r  is the number of particles between = N V

and

V

(12)

 4πr dr r

is the average density of the particles, where

and N

r  dr

from a reference particle.

represents the total number of particles

is the assembly volume. The reference particles form an ergodic system among all the

ACCEPTED MANUSCRIPT particles. As shown in Fig. 9, the RDF is clearly related to the FC. The distance is normalized by dmin , where dmin is the minimum diameter of the fine particles. Three clear peaks, indicating high packing orders, are observed in assemblies with FC greater than 15%. The first nearest-neighbor peak at the

PT

distance of 1.1dmin is clearly greater and sharper than the second and third peaks, indicating that the

RI

packing order is maximized at the first peak. This result indicates that the probability of finding a

SC

particle at a distance that is approximately equal to the minimum diameter of the fines is the highest. However, when the FC is smaller than 15%, the distance corresponding to the first peak is slightly

NU

greater than 1.1, decreasing from 1.2 when FC=2.5% to 1.15 when FC=10% and approaching 1.1 at

MA

FC=15%.

The peak values, marked by solid points in Fig. 9, increase when the FC increases to 30% and

ED

then gradually decrease as the FC continues to increase. When the fines content is extremely low, the

EP T

packing density is lower, and the fine particles may float among the voids among the coarse skeleton. Therefore, in comparison to high-FC samples, less particles are located at the same distance from the reference particles for low-FC samples. This structural characteristic naturally corresponds to the

AC C

softening behavior exhibited at low FCs, as shown in Fig. 5. When sustaining external shear, the particle systems jam [73, 74] to behave like a solid. Because the g(r) values at low-FC conditions are relatively low, the probability of forming jamming structures among the particles in the system is not as high as that for higher-FC samples. According to the void ratio distribution shown in Fig. 2, the densest conditions occurs around FC=30%, and the fine particles completely jam within the voids of the coarse-supported structures. Therefore, the greatest peak is found at FC=30% in Fig. 9. With the continuous increase in the FC and void ratio values from FC=30%, the jamming structures are

ACCEPTED MANUSCRIPT slightly impaired, causing the decrease in the peak g(r) values. In general, the RDFs of granular materials are multimodal [72]. For FC=2.5% and 10%, the second peak in the RDF occurs at a considerable distance from the first peak and is not as sharp as the first peak. The lack of secondary and third peaks for samples with low FCs also illustrate the less

PT

intensive pair correlations. When the FC exceeds 15%, the multimodal characteristic of the RDF

RI

becomes increasingly apparent. The second and third peaks in the RDFs of the FCs greater than 20%

SC

are enlarged in Fig. 9 and emerge at approximately the distances of 2.1dmin and 2.8dmin , respectively. Meanwhile, the secondary and third peak values with respect to the FCs also experience an initial

MA

NU

increase, which is followed by a decrease from the peak points achieved at FC=30%.

5. Contact type characteristics

ED

5.1 Percentages and mechanical contributions within different contact types

EP T

Three kinds of contacts are included in the mixture systems, i.e., the coarse-coarse (c-c), coarse- fine (c- f), and fine- fine (f- f) contacts. The coarse- and fine-dominant zones are mainly occupied by the c-c and f- f contacts, respectively. Therefore, it can be expected that the c-c

AC C

percentages decrease in comparison to the increasing trend of f- f percentages when the FC increases. The evolutions of the percentiles of the c-f contacts during the loading process with varying FCs are depicted in Fig. 10(a). Unlike the monotonically decreasing and increasing trends of the c-c and f- f percentages with the increase in the FC, respectively, the percentage of the c- f contacts fluctuates with the FC; the color bar in Fig. 10(a) indicates this percentage value. For c- f type contacts under large shearing strains, the percentage first increases to a peak value, approximately 50%, when the FC approaches 25% and then decreases when the FC continues to increase. During the shearing

ACCEPTED MANUSCRIPT process, the proportion of c- f contacts exhibits different phenomena at different FCs. When the FC is below 25%, the c-f contacts increase as the shearing evolves. An initial decrease is followed by a continuous increase in the c- f proportions when FC=25%, and a monotonic decrease toward a uniform value is observed for higher FCs. This result indicates that the c-f contacts play a significant

PT

role in the force transmission within the transitional zone.

RI

The contact force of various contact types is another important measurement. As shown in Fig.

SC

10(b), the mean contact force within the c-c contact networks is significantly larger than that within other two contact systems. The mean force of the c-c contact type experiences a sharp increase with

NU

the increase in the FCs from 10%. With an increase in the fines content, the fine particles begin to

MA

separate the coarse particles, and the c-f type contacts, which bridge the fines and coarse particles, gradually participate in sustaining the whole structure. Therefore, the mean force of the c-f contact

ED

type is increased when the FC increases from 10% to 20%. When the FC increases from 20% to 30%,

EP T

the mean force of the c-f contact type slightly decreases, which is attributed to the fact that the increase in the number of fine particles between the coarse particles can reduce the contribution of the coarse particles [14] (as shown in Figs. 3(b) and 3(c)). With a further increase in the fines content,

AC C

the granular assembly is gradually subject to more fine particles, thus making the secondary reinforcement effect of the coarse particles more evident [1]. Hence, the c-f mean force further increases with FCs greater than 30%, as shown in Fig. 10(b). It is interesting to observe that the mean forces of the c- f and f- f contacts are nearly constant, irrespective of the fines content, when the FC is relatively high. Figure 11 depicts the force chain distributions of the c-c and c-f contact networks in the assembly with FC=25%, and the color of the columns represent the force magnitude. Regardless of the greater proportion of the c-f contacts, the c-c contact force transmissions are more

ACCEPTED MANUSCRIPT prominent. The dominant contact type, in terms of relative proportion, does not always undergo the greatest contact force at the microscale. It is also necessary to investigate the sliding fractions of different contact types of binary mixtures. We define the sliding contact fraction as the number of sliding contacts to the total number

PT

of contacts for each contact type. The variations in the sliding contact fraction of different contact

RI

types versus the FCs at a steady state are illustrated in Fig. 12. For the c-c systems, the sliding

SC

fraction increases with the increase in the FC and nearly reaches 100% when FC=30%. In contrast, the f- f sliding fraction exhibits a monotonic decrease as the FC increases. For the fine-dominant

NU

regime, the sliding contact fraction of the f- f system gradually approaches a steady value with the

MA

increase in the fines content. The sliding fraction of the c- f network reaches a minimum value (approximately 50%) at approximately FC=20%. This result demonstrates that for the c- f network,

ED

the sliding contacts dominate for most FCs; the sliding fraction nearly accounts for all the contacts

EP T

when the FC is either extremely high or low.

Based on Eq. (6), the stress tensor deduced from the contact-type partition can also be written as

AC C

follows:

 ij 

1  c -c   fi d V  c  c -c

c -c j



 c  c -f

fi

c -f

d

c -f j



 c  f-f

fi

f-f

d

f-f j

  

(13)

where the forces are additive within different contact types. Their relative contributions to the global mean effective stress are calculated and shown in Fig. 13. This mechanical contributions of the different contact types are similar to their proportion characteristics, as discussed above and shown in Fig. 10(a). 5.2 Coaxiality of the stress tensors and fabrics for each contact-type network The coaxiality of the fabric tensor  ij , characterizing the distribution of contact normal

ACCEPTED MANUSCRIPT direction, and the stress tensor 

ij

has been adopted as one of the critical state indexes to examine

whether the critical state is reached [29]. An anisotropy variable A [75], defined as the joint invariant of the deviatoric fabric tensor and the deviatoric stress tensor

s ij

, is utilized to characterize the

coaxiality and has been incorporated into constitutive work [76]. We employ A in the following

tr  s ik  k j



'

 pq pq '

smn smn s ij  

ij

 p  ij

(14)

and  ij'   ij   kk  ij 3 . The stress and fabric tensors used in Eq. (14) are

SC

where

'

RI

A 

PT

equation [77]:

defined as follows [50]:  ij 

1



MA

k

NU

defined for each contact type. Likewise, the contact-type partition of the fabric tensor  ijk can be

N

k

c N

k

ni n

k j

(15)

k

ED

The three contact-based anisotropy variables are denoted as Ac-c, Ac-f, and Af-f. This partition can help

EP T

us explore the influences of different contact types on the constitutive contribution. Three typical FCs, representing the three zones mentioned in Section 2, are adopted to analyze the contact-based loading direction with respect to the fabric orientation. As shown in Fig. 14, the

AC C

contact-based anisotropy variables, Ac-c at FC=15%, Ac-f at FC=25%, and Af-f at FC=60%, quickly approach 1.0 during the early loading stage and are nearly constant during shearing. The evolving rates of the anisotropy variables for the three contact types are different. In the coarse-dominant zone (FC=15%), the coaxiality is expressed earlier in the c-c contact network and later in the c- f contact network. In comparison, Af-f fluctuates more, and the extent of the coaxiality is remarkably weaker in the f-f network than in the c-c and c-f networks. For the transitional- zone mixture (FC=25%), Ac-f quickly approaches 1.0 at a faster rate than that of Ac-c. The coaxiality within the f- f network is also

ACCEPTED MANUSCRIPT the slowest to be expressed and generally exhibits the lowest coaxial magnitude when FC=25%. For the fine-dominant zone (FC=60%), Af-f is always slightly higher than Ac-f. Therefore, we can conclude that the contact-based coaxial magnitude is subject to the relative proportions of the three contact types under various FCs. The rate of increase in the approach toward the steady value of the

PT

anisotropy variable also depends on the proportions. To further construct the constitutive model of

RI

the mixtures, the different contact-based roles should be comprehensively considered.

SC

5.3 Anisotropy within different contact-type systems

Since the contact normal anisotropy and normal contact force anisotrop y are generally the main

NU

source of strength, the evolutions of a n and a c during the loading process of simulations with

MA

different FCs are presented in the contour depiction of Fig. 15. We find that a c distributes into a valley shape, where the red region represents the highlands. As discussed in Section 2.2, the granular

ED

assembly is the densest and the shear strength (in terms of shear stress ratio) is also the greatest at

EP T

FC=30% because the fine particles surround the coarse particles to the greatest extent at this FC. Hence, this granular assembly is characterized by an insignificant geometrical anisotropy and a very low a c . The a n reaches a peak point at approximately 3% of the a xial strain, and the

AC C

corresponding FC is 30%.

The contact-type classification of the anisotropies is extracted to explore the micromechanical mechanisms related to the different contact types. The contact-based and global anisotropy parameters, extracted at the peak deviatoric stress ratios, are found to closely relate to the FC. Figs. 16(a) and (b) show that the global contact normal and contact normal force anisotropies experience a minimum and a maximum at FC=30%, respectively (marked by a dashed line in Fig. 16). The different anisotropies calculated for the same contact type exhibit nearly the same relationship with

ACCEPTED MANUSCRIPT the FC. Namely, the c-c anisotropies experience a sharp increase in comparison to a gradual increase in the f- f anisotropies. For the c-f anisotropies, a valley distribution among various FCs is clearly observed for both a c and a n . We also observe that the global anisotropies nearly coincide with the anisotropy parameters deduced from the dominant contact type in the coarse-dominated zone,

PT

transitional zone, and fine-dominant zone, as classified in Fig. 13. In other words, when the FC is

RI

lower than 20%, the c-c contact network contributes most to the global anisotropy. As the FC ranges

SC

from 20% to 40%, the global anisotropy parameters are similar to the anisotropy parameters deduced from the c- f contact network. For the fine-dominant zone, the contact normal anisotropy of the whole

NU

assembly is in good agreement with that calculated from the f- f contact network. A similar

MA

relationship is not apparent between the normal contact force anisotropies of the f- f contact network and that of the whole assembly in the fine-dominant zone, but they become closer at higher FCs. It

ED

can be concluded that the anisotropy of granular mixtures depends on different contact types when

EP T

the FC varies. 5.4 Partitioned subnetworks

To connect the shear strength with the anisotropy within various subnetworks, the stress tensor

AC C

can be expressed by the anisotropy tensors [78]: 

k

k

ij



k

Nc f d 3V

k

 2  c ,k 3 t ,k  n ,k   ij   a ij  a ij  a ij  5  2   

2 35

where

 ij

(16)

 ( a p q  a p q ) a p q  ij  ( 3 a iq  4 a iq ) a p j  n ,k

t ,k

c ,k

t,k

n ,k

c ,k

   

is the Kronecker delta.

The anisotropy-based deviatoric stress ratio can be directly obtained from Eq. (16). The simplified stress-force-fabric expression is also given in the form of the anisotropy parameters [79]:

ACCEPTED MANUSCRIPT  q  2 3    ac  an  ad  at    5 2   p  a n is o tro p y

(17)

Sufian et al. [80] further deduced the stress- force- fabric relationship within various partitioned networks. The global stress ratio q/p can be expressed as the summation of (q/p)k that is calculated

q

=

p

,

k

an

,

k

ad

, and

k

k

k

k  at  2 

3

k



k

(18)

are used to measure the degree of anisotropy within the k th subnetwork.

at

The weighting parameter



 k  ac  an  ad 

RI

k

ac

 5

 k   k k



k

 k k

SC

where

2

PT

within all the sub networks in terms of the anisotropy tensors [80]:

 can characterize the contribution of the contact

NU

number and mean contact normal force in the kth subnetwork. The term

MA

the proportion of contacts of the k th sub network, in which the k th subnetwork and global network respectively. fn

k

 k  fn

k

k

and fn

N

k

 N

k

N

is donated as

are the contact number in

represents the ratio of the mean

and that in global system

fn

. Naturally, the k th

ED

contact normal force within k th subnetwork

N



EP T

network can represent the strong-weak networks, the nonsliding-sliding networks, or the contact-type-based subnetworks.

The global anisotropy-based stress ratio derived from Eq. (16) or Eq. (17) agrees well with the

AC C

magnitude of the global deviatoric stress ratio q/p based on the partitioned contact types in Eq. (18), as depicted by the black and red lines in Fig 17(a) and (b). The q/p of the c-c network calculated by Eq. (17) is found to take the largest contribution to global stress ratio when FC=25%. Although the c-f contacts account for the largest proportion when FC=25%, the q/p of the c- f network calculated by Eq. (17) is considerably less than that calculated from the c-c network. The q/p deduced from the f- f contact type is the lowest. The same phenomena can also be observed for the sample within the fine-dominant zone (FC=60%) in Fig. 17(b). The q/p of the f- f network calculated by Eq. (17) are

ACCEPTED MANUSCRIPT nearly the same as that of global network calculated by Eq. (17). Despite the less proportion of c- f contacts when FC=60%, the q/p calculated from c- f system by Eq. (17) is the largest in Fig. 17(b). Therefore, it can be concluded that the coarse-related contacts (c-c for samples within coarse-dominant zone and c-f contacts for samples within fine-dominant zone) play the dominant

PT

role of contributing to the measured global stress ratio. The stress ratio calculated from the contact

RI

type with the largest proportion is identical with the global stress ratio calculated by the

SC

stress-force-fabric relationship.

To further investigate the stress- force- fabric characteristics of granular mixtures, we take the

NU

two-partition networks, the strong-weak and nonsliding-sliding partition methods among the

MA

proportion-dominant subnetworks, into account. As illustrated in Fig. 18(a) (when FC=25%), the shear stress ratios that originate from both the strong-weak and nonsliding-sliding partition networks

ED

in Eq. (18) are similar to the global q/p. When FC=60%, the anisotropy-based q/p of global network

EP T

perfectly coincides with the stress ratios derived from the strong-weak partition network and the nonsliding-sliding network, shown in Fig. 18(b). Hence the stress ratio, derived from both the strong-weak and nonsliding-sliding partition methods within the proportion-dominant contact type,

AC C

can effectively represent the global stress ratio under various FC conditions. As discussed by several previous researches [47, 79], the mechanical anisotropy occupies the largest contribution to global deviatoric stress ratio. For nonsliding-sliding partitioned systems (when FC=25%), the

n o n s li d i n g

an

is still larger than

n o n s li d i n g

ac

as usual, shown in Fig. 19(a). However, when

the strong-weak partitioned method is adopted, the geometrical anisotropy within the strong subnetworks are significantly larger than mechanical anisotropy, which is also verified in the work of Sufian et al. [80]. As shown in Fig. 19(a), the contact normal anisotropy

s tr o n g

ac

within the strong

ACCEPTED MANUSCRIPT subnetwork is evidently larger than the contact normal force anisotropy

s tr o n g

an

. When the contact

system is divided into four partitioned subnetworks (i.e. strong-nonsliding, strong-sliding, weak-nonsliding, and weak-sliding networks), the contribution of the geometrical anisotropy s tr o n g -n o n s lid in g

ac

obtained from the strong- nonsliding network to the mobilized internal friction (i.e. the s tr o n g -n o n s lid in g

ac

and

PT

stress ratio q/p) is proved to be dominant [80]. The linear relationship between the

RI

global stress ratio is verified for various FC conditions in Fig. 19(b). When the fines content keep

SC

increasing, the slope of the fitting line gets increased first, and then experiences a descending trend. To some extent, the prediction of the mobilized internal friction can only be obtained from the s tr o n g -n o n s lid in g

ac

and global q/p, which is dependent on the fines content.

NU

relationship between

MA

To further understand the strong-weak partitions, the distribution of the contact normal direction and contact normal forces derived from the strong-weak networks are depicted in Fig. 20. It is

ED

evident from the weak distributions of both the contact normal direction and normal contact forces

EP T

that they are nearly isotropic. The horizontal contact directional distribution density within the global network (shown in Fig. 20(a)) is greater than that within the strong network (shown in Fig. 20(b)). For contact normal forces, the distribution along the vertical loading direction within the strong

AC C

network and global network are similar. However, the horizontal forces among the strong network are greater than that in the global network. Hence, the contact force distribution within the global network is more anisotropic, as shown in Fig. 20 (d). Hence, we conclude that the participation of the weak contacts can positively contribute to the contact normal force anisotropy but negatively contribute to the contact normal anisotropy.

6. Discussion The coaxiality between the loading direction and fabric orientation for different contact types

ACCEPTED MANUSCRIPT have different characteristics, which should be considered when incorporating the contact-based anisotropy variable, Ac-c for instance, into constitutive work. Within the proportion-dominant contact type, the anisotropy-based stress ratios derived from both the strong-weak and nonsliding-sliding partitioned

networks almost coincide with the global stress ratio calculated from the

PT

stress-force-fabric relationship. Therefore, it is necessary to identify the representative contact

RI

systems when evaluating the strength of a binary mixture. The representative roles of the

SC

strong-weak and nonsliding-sliding partitioned subnetworks will be verified with various loading paths and particle size ratios in our next investigation. The geometrical anisotropy within strong

NU

contact systems accounts for the largest contribution to global mobilized internal friction. Many

Considering

the evident

linear

MA

previous researches have attempted to incorporate the influence of fabric into constitutive model. relationship

between

the

geometrical anisotropy

from

ED

strong-nonsliding network and the global stress ratio, the incorporation of the fabric tensor when

EP T

establishing the constitutive model may be more efficient by means of the geometrical anisotropy based on the strong- nonsliding subnetwork. The strong-weak anisotropy distributions within different contact types may provide an effective way to understand the origin of global anisotropy,

AC C

which will be investigated in our next work. The positive and negative effects of the weak contacts on the mechanical and geometrical anisotropies will also be thoroughly explored in our future work under complex loading conditions.

7. Conclusion We adopt the three-dimensional DEM method and gap-graded numerical samples to explore a range of binary granular mixtures under undrained loading conditions. The influences of the fines

ACCEPTED MANUSCRIPT content on the macro-, meso- and microscale mechanical behaviors are systematically investigated. Both the mechanical and physical characteristics are comprehensively discussed. The main findings are summarized as follows. 1. The local structures of binary mixture are remarkably influenced by the variation in the FC.

PT

Both the local volume fraction, sphericity, and the local anisotropy parameter follow a multimodal

RI

distribution. The PDFs for sphericity and local volume fraction shift to right, respectively. For the

SC

anisotropic Voronoi cells within various mixtures, the PDFs of local anisotropy index tend to coincide with each other when FC increases. The mesoscale structures tend to be stable when FC

NU

increase towards the fine-dominant domain. For a binary mixture, three evident peaks of the RDF

MA

exist at consistent distances, independent of the FC. The intensities of the pair correlations at various FCs exhibit a gradual increase to a maximum at FC=30%, followed by a monotonic decrease as the

ED

FC continues to increase.

EP T

2. The proportions of the c-c, c- f, and f- f contacts exhibit different relationships with the FC. In terms of the force magnitude, the average c-c contact forces are greater than the forces within the c-f and f-f contacts. The relationship between the sliding fraction and FC values within the different

AC C

contact-type networks exhibit different characteristics. The sliding contacts within the c- f systems are dominant in the samples among nearly all the FCs. The mechanical contributions of the three contact types exhibit the same characteristics of their proportions. 3. The coaxiality of the loading direction with respect to the fabric orientation for three contact types are found to depend on their respective proportions. The coaxiality evolution is the fastest for the dominant contact type, and the highest degree of coaxiality is generally reached by the dominant contact type. The contact-based peak state anisotropy parameters calculated from certain contact

ACCEPTED MANUSCRIPT types with the maximum mechanical contributions are similar to the global anisotropy parameters within various FC conditions. The relationship between the FC and the anisotropy parameters within the same contact type are nearly the same but behave differently among different contact types. 4. The global stress ratio q/p perfectly coincides with the q/p calculated from the

PT

contact-type-based partitioned networks. The q/p derived from the proportion-dominant contact type

RI

is also consistent with the global q/p. When further dividing the contact system into strong-weak or

SC

nonsliding-sliding networks within the proportion-dominant contact type, the addition of the anisotropy-based stress ratios multiplied by the weighting parameter match the global q/p results

NU

well. The larger magnitude of the geometrical anisotropy within strong contact network is evidently

MA

observed. Based on the four-partition method, the linear relationship between the strong-nonsliding geometrical anisotropy and the global stress ratio is satisfied under the influence of different fines

ED

contents. The weak contacts are proved to play different roles in geometrical and mechanical

AC C

EP T

anisotropies.

ACCEPTED MANUSCRIPT Acknowledgements This work was financially supported by the National Key R&D Program of China (Grant No. 2017YFC0404801), National Natural Science Foundation of China (Grant No. 51509190 and No.

PT

51579193), and the Fundamental Research Funds for the Central Universities.

S. Thevanayagam, T. Shenthan, S. Mohan, J. Liang, Undrained Fragility of Clean Sands, Silty

SC

1.

RI

Reference

Sands, and Sandy Silts, J. Geotech. Geoenvironmental Eng. 128 (2002) 849–859. Coli, N., Boldini, D., Bandini, A., Lopes, D. S., Modeling of complex geological rock

NU

2.

MA

mixtures under triaxial testing conditions. International Society for Rock Mechanics. In ISRM International Symposium-EUROCK (2012).

W.-J. Xu, Q. Xu, R.-L. Hu, Study on the shear strength of soil-rock mixture by large scale

ED

3.

4.

EP T

direct shear test, Int. J. Rock Mech. Min. Sci. 48 (2011) 1235–1247. R. Salgado, P. Bandini, A. Karim, Shear Strength and Stiffness of Silty Sand, J. Geotech. Geoenvironmental Eng. 126 (2000) 451–462. S. Thevanayagam, S. Mohan, Intergranular state variables and stress–strain behaviour of silty

AC C

5.

sands, Géotechnique. 50 (2000) 1–23. 6.

A. Papadopoulou, T. Tika, The effect of fines on critical state and liquefaction resistance characteristics of non-plastic silty sands, Soils Found. 48 (2008) 713–725.

7.

S.L. Yang, R. Sandven, L. Grande, Instability of sand-silt mixtures, Soil Dyn. Earthq. Eng. 26 (2006) 183–190.

8.

Y. Xiao, H. Liu, Q. Chen, L. Long, J. Xiang, Evolution of particle breakage and volumetric

ACCEPTED MANUSCRIPT deformation of binary granular soils under impact load, Granul. Matter. 19 (2017). 9.

X. Jiang, P. Cui, Y. Ge, Effects of fines on the strength characteristics of mixtures, Eng. Geol. 198 (2015) 78–86.

10.

A.F. Cabalar, The Effects of Fines on the Behaviour of a Sand Mixture, Geotech. Geol. Eng.

P. V. Lade, C.D. Liggio, J.A. Yamamuro, Effects of Non-Plastic Fines on Minimum and

RI

11.

PT

29 (2011) 91–100.

12.

SC

Maximum Void Ratios of Sand, Geotech. Test. J. 21 (1998) 336–347. C.S. Chang, Y. Deng, Z. Yang, Modeling of Minimum Void Ratio for Granular Soil with

Chiu C F, Fu X J., Interpreting undrained instability of mixed soils by equivalent intergranular

MA

13.

NU

Effect of Particle Size Distribution, J. Eng. Mech. 143 (2017) 4017060.

state parameter, Geotechnique, 58 (2008) 751-755. T. Ueda, T. Matsushima, Y. Yamada, Effect of particle size ratio and volume fraction on shear

ED

14.

15.

EP T

strength of binary granular mixture, Granul. Matter. 13 (2011) 731–742. W. Zhou, K. Xu, G. Ma, L. Yang, X. Chang, Effects of particle size ratio on the macro- and microscopic behaviors of binary mixtures at the maximum packing efficiency state, Granul.

16.

AC C

Matter. 18 (2016).

Xiao, Y., Xiang, J., Liu, H., & Ma, Q., Strength–dilatancy relation of sand containing non-plastic fines, Géotechnique Letters, 7 (2017) 204-210.

17.

Y. Xiao, Y.F. Sun, H.L. Liu, J. Xiang, Q.F. Ma, L.H. Long, Model predictions for behaviors of sand- nonplastic-fines mixturesusing equivalent-skeleton void-ratio state index, Sci. China Technol. Sci. 60 (2017) 878–892.

18.

P. V Lade, J. a Yamamuro, Effects of nonplastic fines on static liquefaction of sands, Can.

ACCEPTED MANUSCRIPT Geotech. J. 34 (1997) 918–928. 19.

J.A. Sladen, R.D. D’Hollander, J. Krahn, D.E. Mitchell, Back analysis of the Nerlerk berm liquefaction slides, Can. Geotech. J. 22 (1985) 579–588.

20.

Troncoso, J. H., Verdugo, R., Silt content and dynamic behaviour of tailing sands. In Proc., XI

H.B. Seed, I.M. Idriss, I. Arango, Evaluation of Liquefaction Potential Using Field

RI

21.

PT

Int. Conf. on Soil Mechanics and Foundation Engineering (1985) 1311-1314.

22.

K. Tokimatsu, Y. Yoshimi, Empirical correlation of soil liquefaction based on spt N-value and

NU

fines content., Soils Found. 23 (1983) 56–74.

J. a Yamamuro, P. V Lade, Static liquefaction of very loose sands: Reply, Can. Geotech. J. 34

MA

23.

(1997) 905–917.

R. de Frias Lopez, J. Silfwerbrand, D. Jelagin, B. Birgisson, Force transmission and soil fabric

ED

24.

SC

Performance Data, J. Geotech. Eng. 109 (1983) 458–482.

25.

EP T

of binary granular mixtures, Géotechnique. (2016) 1–6. N.H. Minh, Y.P. Cheng, C. Thornton, Strong force networks in granular mixtures, Granul. Matter. 16 (2014) 69–78.

J. Gong, J. Liu, Mechanical transitional behavior of binary mixtures via DEM: Effect of

AC C

26.

differences in contact-type friction coefficients, Comput. Geotech. 85 (2017) 1–14. 27.

T.-T.. Ng, W.. Zhou, X.-L.. Chang, Effect of particle shape and fine content on the behavior of binary mixture, J. Eng. Mech. 143 (2017).

28.

M.M. Rahman, S.R. Lo, Undrained Behavior of Sand-Fines Mixtures and Their State Parameter, J. Geotech. Geoenvironmental Eng. 140 (2014) 4014036.

29.

X.S. Li, Y.F. Dafalias, Anisotropic critical state theory: role of fabric, J. Eng. Mech. 138

ACCEPTED MANUSCRIPT (2012) 263–275. 30.

Z.Y. Yin, J. Zhao, P.Y. Hicher, A micromechanics-based model for sand-silt mixtures, Int. J. Solids Struct. 51 (2014) 1350–1363.

31.

Kloss C, Goniva C., LIGGGHTS: a new open source discrete element simulation software. In:

PT

Proceedings of the Fifth International Conference on Discrete Ele ment Methods, London, UK.

J. Kumara, K. Hayano, Y. Shigekuni, K. Sasaki, Physical and mechanical properties of

SC

32.

RI

(2010)

sand- gravel mixtures evaluated from DEM simulation and laboratory triaxial test, Int. J.

T. Ng, FABRIC EVOLUTION OF ELLIPSOIDAL ARRAYS WITH DIFFERENT

MA

33.

NU

GEOMATE. 4 (2013) 546–551.

PARTICLE SHAPES, J. Eng. Mech. 127 (2001) 994–999. Krishna P, Pandey D, Taylor C A., Close-packed structures. International Union of

35.

EP T

Crystallography (1981)

ED

34.

J. Ai, J.F. Chen, J.M. Rotter, J.Y. Ooi, Assessment of rolling resistance models in discrete element simulations, Powder Technol. 206 (2011) 269–282. LIGGGHTS user manual. DCS Computing Gmbh. Austria (2015)

37.

P. Jop, Y. Forterre, O. Pouliquen, A constitutive law for dense granular flows, Nature. 441

AC C

36.

(2006) 727–730. 38.

L.E. Vallejo, Interpretation of the limits in shear strength in binary granular mixtures, Can. Geotech. J. 38 (2001) 1097–1104.

39.

T.T. Ng, W. Zhou, G. Ma, X.L. Chang, Macroscopic and microscopic behaviors of binary mixtures of different particle shapes and particle sizes, Int. J. Solids Struct. 0 (2017) 1–11.

ACCEPTED MANUSCRIPT 40.

C.C. Furnas, Grading Aggregates: I—Mathematical Relations for Beds of Broken Solids of Maximum Density, Ind. Eng. Chem. 23 (1931) 1052–1058.

41.

S. Thevanayagam, Liquefaction potential and undrained fragility of silty soils, Proc. 12th World Conf. Earthq. Eng. (2000) 1–8. R.J. Fragaszy, J. Su, F.H. Siddiqi, C.L. Ho, Modeling Strength of Sandy Gravel, J. Geotech.

PT

42.

C.S. Chang, J.Y. Wang, L. Ge, Modeling of minimum void ratio for sand-silt mixtures, Eng.

SC

43.

RI

Eng. 118 (1992) 920–935.

Geol. 196 (2015) 293–304

C.S. Chang, M. Meidani, Dominant grains network and behavior of sand-silt mixtures:

NU

44.

MA

Stress-strain modeling, Int. J. Numer. Anal. Methods Geomech. 37 (2013) 2563–2589. T.W. Lambe, R. V Whitman, Soil Mechanics, Soil Mech. Found. Eng. 365 (1969) 576.

46.

C.S. Chang, Z.Y. Yin, Micromechanical modeling for behavior of silty sand with influence of

ED

45.

47.

EP T

fine content, Int. J. Solids Struct. 48 (2011) 2655–2667. L. Rothenburg, R.J. Bathurst, Analytical study of induced anisotropy granular materials in idealized, Geotechnique. 39 (1989) 601–614. J. Christoffersen, M.M. Mehrabadi, S. Nemat-Nasser, Micromechanical Description of

AC C

48.

Granular Material Behavior., J. Appl. Mech. Trans. ASME. 48 (1981) 339–344. 49.

W. Zhou, L. Yang, G. Ma, X. Chang, Y. Cheng, D. Li, Macro–micro responses of crushable granular materials in simulated true triaxial tests, Granul. Matter. 17 (2015) 497–509.

50.

H. Ouadfel, L. Rothenburg, `Stress-force-fabric’ relationship for assemblies of ellipsoids, Mech. Mater. 33 (2001) 201–221.

51.

K. Iwashita, M. Oda, Mechanics of Granular Materials: An Introduction, in: Fundam. Mech.

ACCEPTED MANUSCRIPT Granul. Mater., 1999: p. 400. 52.

G. Ma, X.L. Chang, W. Zhou, T.T. Ng, Mechanical response of rockfills in a simulated true triaxial test: A combined FDEM study, Geomech. Eng. 7 (2014) 317–333.

53.

T.G. SITHARAM, J.S. VINOD, B.V. RAVISHANKAR, Post- liquefaction undrained

PT

monotonic behaviour of sands: experiments and DEM simulations, Géotechnique. 59 (2009)

W. Zhou, J. Liu, G. Ma, W. Yuan, X. Chang, Macroscopic and microscopic behaviors of

SC

54.

RI

739–749.

granular materials under proportional strain path: a DEM study, Int. J. Numer. Anal. Methods

W. Zhou, J. Liu, G. Ma, X. Chang, Three-dimensional DEM investigation of critical state and

MA

55.

NU

Geomech. 40 (2016) 2450–2467.

dilatancy behaviors of granular materials, Acta Geotech. 12 (2017) 527–540. W. Zhou, W. Wu, G. Ma, Y. Huang, X. Chang, Study of the effects of anisotropic

ED

56.

EP T

consolidation on granular materials under complex stress paths using the DEM, Granul. Matter. 19 (2017) 76. 57.

F. Radjai, D.E. Wolf, M. Jean, J.-J. Moreau, Bimodal Character of Stress Transmission in

58.

AC C

Granular Packings, Phys. Rev. Lett. 80 (1998) 61–64. F. Alonso-Marroquín, S. Luding, H.J. Herrmann, I. Vardoulakis, Role of anisotropy in the elastoplastic response of a polygonal packing, Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 71 (2005). 59.

G. Gong, C. Thornton, A.H.C. C han, DEM Simulations of Undrained Triaxial Behavior of Granular Material, J. Eng. Mech. 138 (2012) 560–566.

60.

D.C. Bobei, S.R. Lo, D. Wanatowski, C.T. Gnanendran, M.M. Rahman, Modified state

ACCEPTED MANUSCRIPT parameter for characterizing static liquefaction of sand with fines, Can. Geotech. J. 46 (2009) 281–295. 61.

Y. Xiao, Y. Sun, H. Liu, F. Yin, Critical state behaviors of a coarse granular soil under generalized stress conditions, Granul. Matter. 18 (2016) 1–13. K. Bagi, Stress and strain in granular assemblies, Mech. Mater. 22 (1996) 165–177.

63.

C.S. Chang, A. Misra, J.H. Xue, Incremental stress-strain relationships for regular packings

RI

PT

62.

64.

SC

made of multi-sized particles, Int. J. Solids Struct. 25 (1989) 665–681. Finney J L., Structure and properties of gra nular materials: guidelines from modelling studies

NU

of liquids and amorphous solids. Advances in the Mechanics and the Flow of Granular

65.

MA

Materials 1 (1983).

Y. Guan, E. Wang, X. Liu, S. Wang, H. Luan, The quantified characterization method of the

ED

micro- macro contacts of three-dimensional granular materials on the basis of graph theory,

66. 67.

EP T

Materials (Basel). 10 (2017).

H. Wadell, Volume, Shape, and Roundness of Quartz Particles, J. Geol. 43 (1935) 250–280. M. Roozbahani, R. Borela, J. Frost, Pore Size Distribution in Granular Material

68.

AC C

Microstructure, Materials (Basel). 10 (2017) 1237. N. Guo, J. Zhao, Local fluctuations and spatial correlations in granular flows under constant- volume quasistatic shear, Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 89 (2014). 69.

McMullen, P., Isometry covariant valuations on convex bodies. SUPPLEMENTO AI RENDICONTI CIRC MAT PALERMO, 50 (1997) 259-271.

70.

G.E. Schröder-Turk, W. Mickel, S.C. Kapfer, F.M. Schaller, B. Breidenbach, D. Hug, K.

ACCEPTED MANUSCRIPT Mecke, Minkowski tensors of anisotropic spatial structure, New J. Phys. 15 (2013). 71.

G.E. Schröder-Turk, W. Mickel, M. Schröter, G.W. Delaney, M. Saadatfar, T.J. Senden, K. Mecke, T. Aste, Disordered spherical bead packs are anisotropic, EPL (Europhysics Lett. 90 (2010) 34001. R.C. Hidalgo, I. Zuriguel, D. Maza, I. Pagonabarraga, Granular packings of elongated faceted

PT

72.

M.P. Ciamarra, M. Nicodemi, A. Coniglio, Recent results on the jamming phase diagram, Soft

SC

73.

RI

particles deposited under gravity, J. Stat. Mech. Theory Exp. 2010 (2010) P06025.

Matter. 6 (2010) 2871.

K. Wang, C. Song, P. Wang, H.A. Makse, Angoricity and compactivity describe the jamming

NU

74.

75.

MA

transition in soft particulate matter, EPL (Europhysics Lett. 91 (2010) 68001. Z. Gao, J. Zhao, X.S. Li, Y.F. Dafalias, A critical state sand plasticity model accounting for

Y.F. Dafalias, X.S. Li, A constitutive framework for anisotropic sand including

EP T

76.

ED

fabric evolution, Int. J. Numer. Anal. Methods Geomech. 38 (2014) 370–390.

non-proportional loading, Géotechnique. 54 (2004) 41–55. 77.

Z. Gao, J. Zhao, Y. Yao, A generalized anisotropic failure criterion for geomaterials, Int. J.

AC C

Solids Struct. 47 (2010) 3166–3185. 78. Chantawarangul K., Numerical simulations of three dimensional granular assemblies. Ph.D. Thesis, University of Waterloo, Waterloo, Ontario, Canada (1993) 79.

N. Guo, J. Zhao, The signature of shear- induced anisotropy in granular media, Comput. Geotech. 47 (2013) 1–15.

80.

A. Sufian, A.R. Russell, A.J. Whittle, Anisotropy of contact ne tworks in granular media and its influence on mobilised internal friction, Geotechnique. (2017) 1–14.

AC C

EP T

ED

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT Fig. 1 (a) Grading curves of assemblies with varying FC values; (b) Binary mixture assembly at FC=25%. Fig. 2 Evolutions of the global, inter-coarse, and inter-fine void ratios with FC. Fig. 3 Mesoscale descriptions of particle structures with different FCs: (a) coarse-dominant zone (underfilled); (b) coarse-supported

transitional zone

(interactive-underfilled);

(c)

fine-supported

transitional zone

(interactive-overfilled); and (d) fine-dominant zone (overfilled).

PT

Fig. 4 Partitions of the contact systems (FC=25%): strong-weak networks and nonsliding-sliding networks. Fig. 5 The macro-mechanical responses under various conditions of fines content: (a) The stress paths; (b) The

RI

evolution of the deviatoric stress.

SC

Fig. 6 The characteristics of the peak deviatoric stress ratio and the equivalent void ratio: (a) The peak deviatoric stress ratio versus FC; (b) The relationship between the peak deviatoric stress ratio and equivalent void

NU

ratio.

MA

Fig. 7 The critical state lines (CSLs) under various conditions of fines contents. Fig. 8 (a) Voronoi tessellation with radical plane polyhedron of the particles and a representative Voronoi cell with

various FCs; (c) PDFs of 

ED

normal vectors on the outer faces; (b) Probability distribution functions (PDFs) of the sphericity R for within Voronoi cells (local anisotropy index); (d) PDFs of the local volume

EP T

fraction  (isotropic quantity).

Fig. 9 Radial distribution functions for various FCs after sample consolidation.

AC C

Fig. 10 (a) The evolutions of the proportion of c-f contacts under different FC values; (b) The mean contact force within different contact-type networks under various FCs. Fig. 11 The force chain networks within different contact networks at steady state (FC=25%): (a) Within a global assembly; (b) Within a coarse-coarse (c-c) network; (c) Within a coarse-fine (c-f) network. Fig. 12 The sliding fractions within different contact-type networks and FCs. Fig. 13 The mechanical contributions of the three contact-type networks to the global effective mean stress. Fig. 14 The evolution of the coaxiality parameter A under different contact-type systems: (a) FC=15%; (b)

ACCEPTED MANUSCRIPT

FC=25%; (c) FC=60%. Fig. 15 The evolutions of the anisotropy parameters under different FCs: (a) Contact normal anisotropy parameter a c . (b) Contact normal force anisotropy parameter a n .

Fig. 16 Contact-based and global peak anisotropy parameters for various FCs: (a) Peak contact normal anisotropy

PT

parameters among various contact-type networks; (b) Peak contact normal force anisotropy parameters

RI

among various contact-type networks.

SC

Fig. 17 The stress-force-fabric relationship calculated from the global network and the different partitioned subnetworks with different contact types: (a) When FC=25%; (b) When FC=60%.

NU

Fig. 18 Comparison between the q/p calculated from the global network and the anisotropy-based stress ratios

MA

based on the strong-weak partition and nonsliding-sliding partition methods: (a) Within c-f network (FC=25%); (b) Within f-f network (FC=60%).

ED

Fig. 19 The characteristics of the geometrical and mechanical anisotropies within strong and nonsliding

EP T

subnetworks: (a) The contact normal and contact normal force anisotropies within the strong and nonsliding networks respectively (FC=25%); (b) The relationship between the global stress ratio and the geometrical anisotropy calculated from the strong-nonsliding subnetworks based on the four-partition method.

AC C

Fig. 20 The three-dimensional directional distributions of the contact normal and contact normal forces among the global network (FC=25%): Contact normal distributions : (a) within the global contact network; (b) within the strong contact network; and (c) within the weak contact network. Contact normal force distributions: (d) within the global contact network; (e) Within strong contact network; and (f) within the weak contact network.

ACCEPTED MANUSCRIPT

Table 1 Specimen properties Cube length

Total particle

Number of coarse

Number of fine

Initial effective

(FC)

(mm)

number

particles

particles

coordination number

0.0%

71.8 72.0 71.8 72.2 72.3 72.5 73.4 75.0 77.6 80.0 72.5 54.7 51.0 48.5

1003

1003

—

5177

1000

4177

9547

1000

8547

19066

1000

18066

29629

1000

28629

41719

1000

40719

55104 70095

1000 1000

54104 69095

4.74 4.69 4.61 4.71 4.67 4.62 4.54 4.89

88953

1000

87953

109954

1000

108954

100000

606

99394

5.96

50000

204

49796

6.03

50850

76

50774

6.05

50000

—

50000

5.81

20% 25% 30% 35% 40% 50% 60% 80% 100%

t

Time step

Particle density

AC C

Sliding friction

Rolling friction

EP T

Coefficient of restitution Poisson’s ratio

Y re 

Mean particle size ratio Initial confining pressure

ED

Table 2 Experimental details of the DEM simulations Parameter Symbol Unit Young’s modulus

GPa

5.32 5.74

Value 65 0.95

0.12 s

µ µr 

RI

15%

SC

10%

NU

5%

MA

2.5%

1e-7 0.35 0.02

kg/m3

2600

kPa

5.65 200

PSR

Initial void ratio

Equivalent void ratio

0.654 0.616 0.579 0.507 0.432 0.356 0.315 0.320 0.336 0.344 0.376 0.415 0.50 0.603

0.654 0.698 0.686 0.676 0.645 0.608 0.612 0.573 0.600 0.596 0.580 0.578 0.576 0.609

PT

Fines content

ACCEPTED MANUSCRIPT Highlights Liquefaction can be intensified by the change of fines content to some extent. Voronoi cells and pair correlation are sensitive to fines content. Strength contributions of different contact types were found.

 

The coaxiality within various contact types depend on their proportions. The contact-based anisotropies behave differently as usual.

AC C

EP T

ED

MA

NU

SC

RI

PT

  

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

Figure 19

Figure 20