A tensor-based analysis of stress variability in granular media subjected to various loading conditions

Powder Technology 356 (2019) 581–593

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

A tensor-based analysis of stress variability in granular media subjected to various loading conditions Xin Huang a, Qinghua Lei b, Zixin Zhang a, Hui Qin c,⁎ a b c

Key Laboratory of Geotechnical and Underground Engineering, Ministry of Education, Tongji University, 1239 Siping Road, Shanghai 200092, China Department of Earth Sciences, ETH Zürich, Zürich, Switzerland School of Civil Engineering, Dalian University of Technology, Dalian 116024, China

a r t i c l e

i n f o

Article history: Received 30 May 2019 Received in revised form 20 August 2019 Accepted 31 August 2019 Available online 01 September 2019 Keywords: DEM Critical state Stress variability Local stress fluctuation Multi-fractal analysis

a b s t r a c t This paper uses a novel tensor-based approach to characterize the spatial variability and localization of stress fields in granular media derived from three-dimensional discrete element simulations. We use the Euclidean distance of the stress tensor of each particle to the bulk mean stress tensor to quantify local stress fluctuations and use the effective variance of the stress field to quantify the overall dispersion. We observe that the evolutions of both the local stress fluctuation and the bulk stress dispersion strongly depend on the stress state, packing density and shearing process. However, they become constant, and are uniquely related to the deviatoric stress and the void ratio when the sample is at critical state, under which the stress variability reaches an ultimate condition. Multifractal analysis further shows that the local stress field is highly heterogeneous and exhibits self-organized patterns. © 2019 Elsevier B.V. All rights reserved.

1. Introduction A granular material is an assembly of particles interacting at their contacts. When subjected to deviatoric shearing, the particles and forces within such particulate systems often arrange and form complex and heterogeneous patterns [1]. In particular, some of the highly-stressed particles tend to aggregate and align along the deviatoric loading direction, forming chain-like structures which sustain the majority of external loads. These chain-like structures are normally referred to as force chains in the literature. The heterogeneity of force chains within a granular media when subjected to external loading has been extensively studied both experimentally [2,3] and numerically [3–6]. The pattern of a force transmission network has been shown to evolve in accordance with external loading process [7,8], and have close correlations with the mechanical properties of granular systems, including the strength and deformation [9,10] as well as acoustic propagation [11,12]. The heterogeneity of force chains was normally evaluated by considering the spatial anisotropy of contact orientation and intensity of contact forces. Specifically, quantification of anisotropies of force chains in the soil mechanics community was mostly based on the stress-forcefabric analytical solution developed by Rothenburg & Bathurst [13]. Within such a framework, the anisotropy sources of force chains were normally decomposed into a geometrical group and a mechanical

⁎ Corresponding author. E-mail address: [email protected] (H. Qin).

https://doi.org/10.1016/j.powtec.2019.08.107 0032-5910/© 2019 Elsevier B.V. All rights reserved.

group [14–17]. The former describes the spatial distribution of contacts, while the latter denotes the spatial intensity of both normal and tangential contact force magnitudes. These anisotropy parameters have been shown to evolve in accordance with external loading, and provide fundamental insights into the triggering mechanism of several important phenomena, such as static [15] and cyclic liquefactions [18], and structural features of soils at critical state [19]. Despite its effectiveness of linking meso-scale structure and macro-scale stress-strain behavior, this approach, however, considers the orientation and magnitude of contact forces separately, which cannot adequately represent the tensorial nature of particle stresses (related to local contact forces). Furthermore, the localized characteristics of stress variability are difficult to assess based on such a homogenization approach. In fact, the external loading is sustained by individual particles through interaction at the contacts. Therefore, the spatial heterogeneity of force transmission within granular media can also be elucidated by considering the variability of individual particle stress, which is directly linked to the contact forces [20,21]. The overall variability and its localized feature of particle stress can be quantified using the Euclidean approach proposed by Gao & Harrison [22,23] in which the Euclidean distance of the local stress tensor to the bulk mean stress tensor was used to quantify local stress fluctuations, while the overall dispersion of stress field was quantified using its effective variance. Through such tensor operations, the tensorial nature of stresses was retained. The effectiveness of such approach in correlating stress variability with fracture geometry as well as far-field stress conditions was demonstrated by Lei & Gao [24,25].

582

X. Huang et al. / Powder Technology 356 (2019) 581–593

In this study, the discrete element method (DEM) is used to simulate typical loading paths in soil element testing for determining the mechanical properties of soils. Stress variability within the DEM samples and the evolution of its localized feature during the entire shearing processes are assessed using the tensor-based approach proposed by Gao & Harrison [22,23]. The spatial heterogeneity of local stress fluctuations is assessed by multifractal analysis. Special attentions are paid to these correlations at critical state [26], at which the sample deforms continuously with no changes in stress and volume. The critical state is the ultimate state that a soil can attain and the concept of critical state has become the footstone of modern constitutive modeling of soil behavior. We aim at, through the current study, assessing whether the properties of stress field also become invariant at critical state and exploring how they are correlated with the mechanical properties and conditions of the soil samples. 2. Fundamentals of tensor-based quantification of stress variability We analyze the stress data using a tensor-based mathematical formulation [22], which faithfully honors the tentorial nature of stress information and overcomes the drawbacks of conventional decoupled analysis of stress magnitude and orientation [24]. We quantify the local stress fluctuation based on the Euclidean distance between a local stress tensor and the mean of the entire stress field, which measures how far away the local stress tensor is with respect to the average stress state (see Section 2.1). We also evaluate the overall stress variability of the entire stress field by calculating the effective variance, which gives a scalar-valued measure of how spread out a stress tensor group is with respect to their mean (see Section 2.2). The two parameters will together provide a comprehensive and quantitative characterization of the stress variability in granular media in both local and global perspectives. Further justification of the advantages of this tensor-based stress characterization method can be found in [23–25]. 2.1. Characterization of local stress fluctuation For two stress tensors, Si and Sj, their Euclidean distance can be quantified using the Frobenius norm between them [22–24], that is: 
  d Si ; S j ¼ Si −S j  F

ð1Þ

where ║•║F denotes the Frobenius norm, which is also known as the Euclidean norm or Hilbert-Schmidt norm [27]. Intuitively, if Sj in Eq. (1) is replaced by the mean tensor S of the entire stress field, Eq. (1) can be used to quantify the local stress fluctuation in a specified region or around a specified particle, i.e., dðSi ; SÞ ¼ kSi −Sk F . Here, dðSi ; SÞ is a tensor-based parameter reflecting the uniformity of stress distribution within a granular system: it is zero for an ideally uniform stress field but is variable for a heterogeneous stress field.

the mean value of all the column vectors within a granular system, P i.e., s ¼ 1=n∙ n1 si , where n is the total number of local stress tensors. A larger Ve(S) indicates a more dispersed stress field. Ve(S) has been successfully used by Lei & Gao [24] to evaluate the bulk stress variability of fractured rocks and is adopted herein to quantify the bulk stress variability within granular systems. 3. Overview of DEM simulations We modelled the geomechanical responses of granular media based on an open-source code, LAMMPS [28]. The particle size distribution (PSD) approximated that of Toyoura sand (Fig. 1). Initially noncontacting rigid spherical particles (20,164 in total) were generated within a cubic periodic cell. A simplified Herz-Mindlin contact model was used with the particle shear modulus and Poisson's ratio taken to be 29 GPa and 0.12, respectively, approximating the values of quartz [29]. The sample was firstly compressed isotropically until the prescribed confining pressure has been reached. Different interparticle friction coefficients were applied during this procedure to obtain samples with different void ratios (e) subjected to prescribed confining pressures. Interparticle friction coefficient was set to be 0.25 before shearing for all the samples. Three types of loading conditions were simulated, including the traditional drained, undrained and constant mean effective stress (i.e. p’ = (σ’1 + σ’2 + σ’3)/3 is kept constant) tests. Specifically, the lateral confining stress remained constant throughout the drained simulations (CDC); in the undrained simulations (CVC), the overall sample volume was maintained unchanged by continuously adjusting the positions of lateral boundaries in accordance with the movement of vertical boundaries; p’ was kept constant in the constant p’ simulations (CPC). DEM samples of selected void ratios in the CVC simulations showed very dilative to limiting contractive responses, which can well replicate the typical undrained behaviors of medium dense to dense sands. The loose region of the critical state diagram was determined by drained simulations. More simulation details can be found in [30]. In total, 14 numerical experiments were conducted and the simulation information is given in Table 1, in which e0, p’0, ecs, p'cs and qcs denote the initial void ratio, initial mean effective stress, void ratio at critical state, mean effective stress at critical state and deviatoric stress (i.e. q = σ’1 − (σ’2 + σ’3)/2) at critical state, respectively. 4. Macroscopic stress-strain behavior Fig. 2 shows the stress-strain curves for three drained simulations which have the same confining pressure but different initial void ratios. For all the three samples, the deviatoric stress increases rapidly at the beginning of shearing. The slope of the densest sample CDC-3 is the

2.2. Calculation of overall stress dispersion The overall dispersion of the stress field within a granular system (S) can be assessed by its effective variance: V e ðSÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 mðm þ 1Þ detðΩÞ 2

ð2Þ

in which m is the dimension of the stress tensor field S, and Ω is the covariance matrix of the stress vector field s given by [23]: Ω ¼ covðs; sÞ ¼

n h i 1X ðs −sÞðsi −sÞT n i¼1 i

ð3Þ

where si = [σxx,i σyx,i σyy,i σzx,i σzz,i]T is the column vector which stacks only the lower-triangular components of a local stress tensor, and s is

Fig. 1. Comparing the grading curve of DEM samples and that of Toyoura sand.

X. Huang et al. / Powder Technology 356 (2019) 581–593 Table 1 Summary of DEM simulations. Simulation type

ID

e0

p’0 (kPa)

ecs

p'cs (kPa)

qcs (kPa)

Drained

CDC-1 CDC-2 CDC-3 CDC-4 CDC-5 CDC-6 CVC-1 CVC-2 CVC-3 CVC-4 CVC-5 CPC-1 CPC-2 CPC-3

0.588 0.606 0.646 0.627 0.625 0.625 0.578 0.603 0.604 0.612 0.620 0.643 0.625 0.612

500 500 500 1000 2000 5000 1000 1000 5000 5000 5000 419.8 5000 10,000

0.630 0.630 0.629 0.628 0.627 0.618 0.578 0.603 0.604 0.612 0.620 0.632 0.619 0.612

646.7 646.9 646 1286 2580.7 6488.1 35,757 15,548 14,006 9061.3 5684.7 419.8 5000 10,000

440.1 440.7 437 865.4 1741.7 4463.8 24,265 10,491 9506.8 6169.5 3864.6 284.9 3444.6 6877.8

Undrained

Constant p’

largest, while the slope of the loosest sample CDC-1 is the smallest. A peak q value is attained at 8.06%, 8.90% and 14.82% axial strain (εa) for CDC-1, CDC-2 and CDC-3, respectively. The loosest sample CDC-3 contracts throughout the simulation, while the densest sample CDC-1 is dilative over the entire loading process. Response of the medium dense sample CDC-2 is intermediate: it initially contracts but becomes dilative immediately after around 2% axial strain. The q-εa and εv-εa curves of all three samples become approximately constant and converge at large strain levels, i.e., an identical critical state is reached. Note that within the framework of critical state soil mechanics (CSSM), critical state refers to the special feature of soils which once have been sheared to a large strain level will deform continuously with no changes in volume and stress state [26]. Fig. 3 shows the stress-strain behaviors of three constant-volume (undrained) simulations CVC-2, CVC-3 and CVC-5. Note that CVC-2 and CVC-3 have similar initial void ratios but different confining pressures, while CVC-3 and CVC-5 have the same confining pressure but different void ratios. For CVC-2, q increases continuously and becomes approximately constant after around 30% axial strain and its p’ value also increases continuously and becomes approximately constant in accordance with q. For CVC-3, the q value increases initially and reaches a local maximum at around 0.79% axial strain and increases again until it becomes roughly constant after the axial strain is beyond around 30%. The local maximum of q marks the instability state (IS) [31]. p’ of CVC-3 decreases initially due to contraction but starts to increase at around 3.47% axial strain until the critical state is reached. The transition of decreasing p’ to increasing p’ indicates that the behavior of CVC-3 changes from a contractive manner to a dilative manner, which is thereby marked as a phase transformation (PT) state [32]. q and p’ values at the critical

(a)

583

state of CVC-3 are slightly smaller than those values of CVC-2 because CVC-3 is slightly looser than CVC-2 (see Table 1). The stressstrain responses of CVC-5 are similar to those of CVC-3, sequentially experiencing an IS, a PT state and a critical state. The phase transformation phenomenon of CVC-5 is more pronounced than CVC-3 and the critical-state q and p’ values of CVC-5 are much smaller than those of CVC-3 because CVC-5 is looser than CVC-3. Fig. 4 shows the evolutions of q and volumetric strain (εv ) with the axial strain for a representative constant-p’ simulation, i.e. CPC-1. q increases from 0 to 307 kPa at 14.52% axial strain and remains approximately constant thereafter. ε v increases rapidly at the beginning of shearing, starts to decrease at around 7.26% axial strain, and becomes approximately constant beyond 24.46% axial strain, where the critical state is attained. The volumetric strain shows a sudden drop at around 52% axial strain while the deviatoric stress remains steady. The reason may be two-fold: firstly, the constant-p’ servo control was updated every 1000 timesteps. This frequency may be too sparse when the critical state was reached. Furthermore, the stress-strain responses were recorded at intervals of 0.5% axial strain. The real volumetric strain may oscillate between the two step points around 52% axial strain and may go up again if the sample is further sheared. This transient drop of volumetric strain does not violate the critical state approximation and could be compensated as the stress-strain values averaged over the last 5% axial strain are taken to be the critical-state values. The stressstrain behavior of CPC-1 is overall contractive, indicating an initially loose condition. Due to the limiting space, only the details of some representative simulation data are presented in Figs. 2 to 4. Fig. 5 compiles the general evolution paths of all the 14 simulations in q-p’ and e-p’ spaces. The values at critical state are determined by corresponding values averaged over the last 5% of axial strain and are indicated by solid circles outlined by different colors. Considering the slight oscillations in the stress-strain responses, the critical-state values are determined by averaging corresponding quantities over the last 5% axial strain. As Fig. 5a shows, the critical-state data in q-p’ space can be fitted by a linear relationship with a slope of 0.697 regardless of the initial states and loading conditions. The critical-state e values and p’ values normalized by the atmospheric pressure (p a = 101.325 kPa) can be well represented by a power-law regression curve. A power of 0.7 is selected according to [33] for quartz-based sand. Note that the selection of the power value may alter other fitting parameters but would not violate the general linear trend. The intermediate stress ratio plays an important role in mechanical behaviors of granular media and its influences on the micro and micro-scale critical-state responses have been profoundly investigated by [14,16,17] among others.

(b)

Fig. 2. Stress-strain curves of drained simulation with the same confining pressure: (a) q against εa; (b) εv against εa (positive εv indicates contraction, and vice versa).

584

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

(b)

Fig. 3. Stress-strain curves of representative undrained simulations: (a) q-εa; (b) q-p’ (IS denotes instability state and PT means phase transformation state).

5. Stress variability and localization characteristics 5.1. Local stress fluctuation Since no stress is present within a void, the local stress tensor of a granular system can be represented by the average stress tensor of individual particles. For a particle (p) in equilibrium, its average stress ðpÞ

tensor σ ij can be calculated from the contact force through the following formula [21]:  1 X  ðcÞ ðpÞ  ðc;pÞ ðcÞ σ ðpÞ ij ¼ − x −xi ni F j Vp N i

ð4Þ

c

where Vp is the particle volume, Nc is the number of contacts that the th (p) particle has, x(c) i and xi are the i coordinate components of the contact (c) and particle (p), respectively, n(c,p) is the unit-normal vector directed i th from the particle center to the contact location and F(c) j is the j component of contact force at contact (c). Fig. 6 shows the probability distribution function (PDF) of the Euclidean distance between a local stress tensor and the mean stress tensor at characteristic states of three representative simulations, CDC-1, CVC-5 and CPC-1. For all the simulation instances considered, the PDF can be fitted by a lognormal regression function. The fitting parameters of these lognormal regression functions are summarized in Table 2. All the fitting curves have a R-square value higher than 0.91, which is large enough to show that the fitting curves and the original data have good correlations. Such a lognormal correlation is not

Fig. 4. Stress-strain curves of a representative constant-p’ simulation CPC-1.

invariant but changes during the process of shearing, which may possibly be attributed to the change of stress state and packing density. We calculate the geometric mean and standard deviation of the lognormal-type PDFs of local stress fluctuation. As can be noted in Fig. 6, the degree of local stress fluctuation depends on the stress state and packing density. The PDFs at different loading stages of CDC-1, CVC-5 and CPC-1 can be fitted by different lognormal functions. The discrepancy of PDFs at different loading stages of CVC-5 is the most notable, while the discrepancy of PDFs at different loading stages of CPC-1 is the least notable and the fitting curves almost overlap. Subject to deviatoric shearing, the particles within a granular system tend to aggregate and form some distinct force transmission networks. Apart from the major force transmission network which contains the most particles and sustain the majority of deviatoric load, there are also some smaller clusters. These two sets of particle clusters can be identified using a modified approach of Huang et al. [34] by considering only strong contacts, i.e., the contacts whose normal contact force equals to or is larger than the average normal contact force. Specifically, an individual particle is selected and its contacting neighbours are identified by searching through the neighbour list which only contains strong contacts. The same procedure is conducted recursively to all of its contacting particles to identify all the particles comprising the cluster. Then this procedure is repeated for particles not already members of a cluster until all particles have been either associated with a cluster or found to be isolated particles. The cluster containing the most particles is the major force transmission network (MF). Fig. 7 shows evolutions of geometrical mean and standard deviation of local stress fluctuation for three representative simulations (CDC-1, CVC-5 and CPC-1) considering different groups of particles, i.e., particles comprising the major force transmission network (MF), all the particles (All) and particles which are not in the major force transmission network (Others). In comparison to Figs. 2 to 4, the geometric mean of local stress fluctuation generally increases as q increases and becomes approximately constant once the critical state is reached. Stress fluctuation is mostly concentrated in particles comprising the major force transmission network, which is partially compensated by that of rest clusters. MF is dominant as the overall local stress fluctuation (All) is closer to MF. In all the three simulations, the standard deviation of local stress fluctuation of MF and All increases initially to a peak value and drops thereafter until approximately constant values are reached at critical state, while the standard deviation of local stress fluctuation of particles not in MF decreases initially to a valley point and increases thereafter until the critical state is reached. Variation of local stress fluctuation is more notable in MF than in the rest clusters. Since the trend of local stress fluctuation determined from MF is similar to that considering all the particles, discussions hereafter will be focused only on the characteristics of local stress variance considering all the particles.

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

585

(b) Fig. 5. Critical state lines in (a) q-p’ space; (b) e-p’ space.

Fig. 8 illustrates the evolution of the geometric mean and standard deviation of local stress fluctuation with e and q in 3D space considering all the particles, where the solid circles mark the values at critical state which are obtained by averaging the corresponding values over the last 5% axial strain. The evolutions of both geometric mean and standard deviation do not follow the same trace between different simulations. For the CVC samples, a general trend may be that the geometric mean increases but its standard deviation decreases as q increases. No obvious correlation between the degree of local stress fluctuation and stresses or

packing density can be observed for the CDC and CPC simulations as both stresses and packing density change simultaneously in these simulations. To isolate the influence of q on the degree of local stress fluctuation from that of packing density, correlations of geometric mean and standard deviation of PDF with q for the five constant volume simulations are presented in Fig. 9a and b, respectively. For the two samples showing a consistently dilative response, i.e. CVC-1 and CVC-2, the degree of local stress fluctuation varies linearly with q in the bi-logarithmic space,

(a)

(b)

(c) Fig. 6. Probability density distribution of local stress fluctuation for three representative simulations: (a) CDC-1; (b) CVC-5; (c) CPC-1 (Peak corresponds to peak state, IS denotes instability state, PT indicates phase transformation state and CS means critical state).

586

X. Huang et al. / Powder Technology 356 (2019) 581–593

Table 2 Statistical parameters of lognormal type PDFs of local stress fluctuation for three representative simulations at different loading instances. Statistical parameters

Logarithmic geometric mean, μ (MPa) Logarithmic variance, h (MPa) R-square

CDC-1

CVC-5

CPC-1

εa = 0%

εa = 2.4%

Peak

CS

εa = 0%

IS

PT

CS

εa = 0%

εa = 2.4%

Peak

CS

0.413 0.392 0.979

0.488 0.394 0.932

0.622 0.346 0.934

0.504 0.396 0.958

2.334 0.368 0.977

2.197 0.371 0.975

1.935 0.365 0.955

1.677 0.386 0.972

0.344 0.057 0.919

0.352 0.061 0.918

0.360 0.071 0.951

0.360 0.065 0.945

which implies that dðSi ; SÞ is in a power-law relationship with q. For the rest samples which initially contract but then dilate, the degree of local stress fluctuation varies with q following a bilinear trend before and after PT state. The slope before PT state is larger than that after PT state, which indicates that degree of local stress fluctuation is more

sensitive to q before PT state than after PT state. This may be due to the fact that the mean effective stress drops before PT state but increases after PT state. The latter tends to compact the sample more isotropically, which thereby surpasses the trend of increasing heterogeneity due to increasing q. Furthermore, the slopes of these linear correlations seem

(a) CDC-1

(b) CVC-5

(c) CPC-1 Fig. 7. Evolutions of geometrical mean and standard deviation of local stress fluctuation for representative simulations: (a) CDC-1; (b) CVC-5; (c) CPC-1 (MF: major force transmission network; All: All particles; Others: particles not in MF).

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

587

(b)

Fig. 8. Evolutions of (a) geometrical mean and (b) standard deviation of local stress fluctuation with e and q (solid circles represent values at critical state).

to slightly depend on the initial packing density. The variance of local stress fluctuation does not show a correlation with q as obvious as the degree of local stress fluctuation. Despite some variability, as Fig. 9b shows, the variance of local stress fluctuation generally decreases with increasing q. Although the evolution of local stress fluctuation characteristics depends on initial state and loading condition, unique relationships between these quantities and q as well as unique relationships between these quantities and e can be identified at critical state. Specifically, as Fig. 10 shows, the geometric mean of local stress fluctuation at critical state increases linearly with increasing q but decreases linearly as e increases. In contrast, the variance of local stress fluctuation correlates with both q and e in a power-law manner with a deviation term, such that the curves in Fig. 11 are nonlinear. Opposite to the geometric mean of local stress fluctuation, the variance of local stress fluctuation declines as q increases but increases with increasing e. This is understandable as referring to Fig. 5 increasing q at critical state is associated with increasing p’, but an opposite variation is observed for increasing e. A larger p’ tends to bring a granular system to a more isotropic state, which thereby reduces the heterogeneity of stress distribution. Note that the critical-state stress ratio has been shown to vary with the intermediate stress ratio [16,17]. Since the intermediate stress ratio does not change the correlations between bulk stress and deformation at critical state, it is reasonable to anticipate that the intermediate stress ratio will also not alter the form of correlations presented in Figs. 10 and 11 but the fitting parameters may vary systematically with the intermediate stress ratio. This needs to be further investigated in the future studies.

Eq. (1) only gives the magnitude of local stress fluctuation. The relative orientation (co-axiality) between a local stress tensor Si and the bulk stress tensor SB can be assessed using the normalized quantity of their first joint invariant [15,16]: A0 i ¼

Si  SB kSi k F kSB k F

ð5Þ

The mean value of A'i, A0 , quantifies the average deviation of local stress orientation from the bulk stress orientation. As Fig. 12 shows, in all the three simulations considered, A0 initially decreases due to the abrupt change of isotropic loading mode to deviatoric loading mode. Evolution of A0 does not show apparent correlations with the evolution of stress-strain responses except for CVC5. For CVC5, evolution of A0 follows a similar trend as that of q. A0 in all simulations becomes approximately constant at critical state. Overall, A0 is maintained at values higher than 0.81, which indicates that local stress tensor is approximately coaxial with the bulk stress tensor. This is understandable as the bulk stress tensor is essentially the volumetric average of the summation of particle stress tensors weighted by solid volume of corresponding particles [21]. 5.2. Overall stress variability Fig. 13 shows the evolution of effective stress variance Ve with axial strain εa at representative simulations. For the three drained

PT

(a)

(b)

Fig. 9. Relationships of local stress fluctuation characteristics with q for CVC simulations: (a) geometric mean; and (b) standard deviation.

588

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

(b)

Fig. 10. Relationships of geometric mean of local stress fluctuation at critical state with: (a) q; and (b) e.

simulations with the same confining pressure, i.e., CDC-1, CDC-2 and CDC-3, Ve increases initially and reaches a peak value at εa = 2.96%, 7.26% and 5.08%, respectively. These strain levels at which the peak Ve values are reached are different from those at which the peak stresses are attained (see Fig. 2). This is probably because both the deviatoric stress and void ratio change during drained simulations. Both of them affect the stress distribution within the sample as discussed in Section 5.1. Ve values become approximately constant when εa exceeds 20%. Evolutions of Ve against εa in the three representative undrained simulations CVC-2, CVC-3 and CVC-5 are presented in Fig. 13b. The Ve-εa relationship is very similar to the q-εa relationship shown in Fig. 4. In particular, the Ve values at critical state in the two simulations CVC-2 and CVC-3 which have a similar e are close to each other. In the CVC-5 simulation, Ve increases initially and reaches a local maximum at IS state. It decreases thereafter and approaches a valley point at PT state, after which it starts to increase until it becomes approximately constant after about 30% axial strain. The Ve values are significantly smaller in the CVC-5 simulation than those in the CVC-2 and CVC-3 simulations. For the constant-p’ simulation CPC-1, the Ve-εa relationship can be divided into three stages in accordance with the stress-strain behavior shown in Fig. 4. In Stage I (εa is smaller than 10%) where both q and e change, Ve initially increases and then decreases during shearing. In Stage II (εa is between 10% and 22.86%) in which q becomes saturated but e continues to increase, Ve decreases accordingly. Ve decreases faster

(a)

in Stage I than in Stage II because q increases but e decreases in Stage I, both of which, as shown in Fig. 10, will intensify the decline trend of stress variability. When εa exceeds 22.86% (Stage III), Ve becomes approximately constant. This is because both q and e become approximately constant at critical state beyond such a strain level. As having been noted in Section 5.1, the local stress dispersion depends on both q and e. To isolate the influence of q, Fig. 14 plots Ve against q for all undrained simulations. Both the Ve-q relationships before and after the PT state can be represented by a power-law relationship. In particular, the slope of the Ve-q relationship before the PT state is larger than that after the PT state. This may be also because the mean effective stress p’ decreases before PT state but increases after PT state. The former exaggerates the anisotropy of force transmission network induced by deviatoric loading. Furthermore, the slope of the Ve-q relationship after PT state seems to increase with increasing e. As Fig. 15 shows, despite that the evolution of Ve depends on the initial states and loading paths, similar to the geometric mean of local stress variance, Ve at critical state increases with increasing q following a unique power-law relationship. In contrast, Ve at critical state decreases as e increases. The Ve-e relationship at critical state can be represented by a unique power-law regression curve with a deviation term regardless of initial states and loading paths. Intuitively, although the bulk stress and volume remain unchanged at critical state, the particles continue to move and the internal structure continues to evolve as the deviatoric deformation keeps developing, which will also alter the distribution of stress field within granular media. Ve reflects the integrated

(b)

Fig. 11. Relationships of variance of local stress fluctuation at critical state with: (a) q; and (b) e. X-Z projection of distribution contour of local stress fluctuation within the middle slice of samples at different stages of simulation.

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a) CDC-1

589

(b) CVC-5

(c) CPC-1 Fig. 12. Evolution of co-axiality between local stress tensor and bulk stress tensor.

effect of change of local stress variability on the overall stress dispersion. As Fig. 13 shows, the overall stress dispersion is statistically stable and constant at critical state in accordance with bulk behavior shown in Figs. 2 to 4. Fig. 15 further shows that the critical-state stress field is uniquely related to the bulk properties regardless of initial stress state, packing density and the loading path. These correlations enable a more profound understanding of critical state for granular media with reference to stress transmission. Some data points at high void ratios seem to be slightly off the fitting curve. This may be because in these simulations the fraction of rattlers is relatively high. Most of these rattlers are very lightly loaded, which makes the calculated stress variance deviate from the main trend. The spatial heterogeneity of the local stress fluctuation field can be further analyzed considering the Q-th order fractal dimension [34–36] which is defined as: DQ ¼ lim r→0

log

PN i¼1

logr

P Qi

ð6Þ

in which Pi is the proportion of stress fluctuation within the cubic box of size r in 3D and N is the total number of boxes. In particular, D0 is called the capacity dimension which equals to the dimension of the embedding space, i.e. 3 in the current problem. A larger DQ value at each positive Q-th order indicates a more dominant influence of sub-boxes where local stress fluctuation is more pronounced; whereas, a larger DQ value at each negative Q-th order indicates a more dominant influence of sub-boxes where local stress fluctuation is less pronounced. The closer the DQ values are to 3, the more homogeneous the distribution of the analyzed quantity within the analysis domain will be. Since the sample size changes during simulation, we define the domain for box-counting calculation as an extended region ranging from

the middle bottom of the sample upwards for 0.45 times the initial sample height (H0), and transversely to each side for 0.225 times the initial sample width (B0) (Fig. 16a). Variations of the fractal dimension with the order Q at characteristic instances of the three representative simulations, i.e. CDC-1, CVC-5 and CPC-1, are shown in Fig. 16b to c. In the CDC-1 case, when Q is positive, at all loading instances DQ decreases rapidly when Q is smaller than 10. It continues to drop as Q is further increased and becomes approximately constant when Q exceeds 50. The local stress variance at peak state is dominated by sub-boxes where local stress fluctuation is not very obvious as its stabilized DQ of positive Q values is the smallest and its stabilized DQ of negative Q values is the largest among all characteristic instances. This can also be inferred from Fig. 17, which shows the projection of local stress fluctuation contour in the X(horizontal)-Z(vertical) plane within the fixed analysis domain of Fig. 16a at selected loading stages. The contact force network is also overlaid in Fig. 17, in which the contacts are represented by lines connecting the centroids of contacting particles and the line thickness is proportional to the magnitude of corresponding contact normal force. As Fig. 17a shows, although there are some highly concentrated areas of Pi at peak state, a large proportion of the analysis domain is occupied by regions with Pi values lower than 0.0005. In comparison to that at peak state, the local stress variability at initial state is more dominated by sub-boxes where the local stress tensor deviates intermediately from the global stress tensor as its DQ of positive Q values is the largest and its DQ of negative Q values is medium among all characteristic instances. This can also be reflected by Fig. 17a which shows the distribution of local stress fluctuation at initial state is the most homogeneous with most Pi values approximating an intermediate value of 0.001. The number of sub-boxes with medium to large local stress fluctuations at critical state is closer to that at initial state. However, there are slightly more sub-boxes with small local stress fluctuation at critical

590

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

(b)

(c) Fig. 13. Evolution of effective stress variance Ve with axial strain εa: (a) representative drained simulations with the same confining pressure; (b) representative undrained simulations; (c) representative constant-p’ simulation CPC-1.

state than at initial state. Therefore, DQ of positive Q values at critical state is larger than that at initial state, but DQ of negative Q values at critical state is close to that at initial state. It can also be inferred from Fig. 17 that the regions where local stress fluctuations are notable coincide with those regions where the contact force network is intensified. In these regions, the contact force is relatively large and the particles bear higher stress than particles in other regions. Consequently, as Fig. 7 shows, they experience higher stress variabilities than particles in other regions. For the undrained CVC-5 simulation (Fig. 16c), DQ of positive (negative) Q values decreases (increases) rapidly when Q is small and

becomes approximately constant when Q exceeds 20 (or −20). The local stress field is the most homogeneous at initial state as its DQ of positive Q values is the largest and its DQ of negative Q values is the smallest among all characteristic instances. This can also be inferred from Fig. 17b which shows that a majority of Pi values at initial state are close to an intermediate value of 0.001. Distribution of Pi at IS is similar to that at initial state. The spatial variability of local stress fluctuation at critical state is dominated by subboxes with small Pi. Therefore, DQ at critical state deviate the most from 3 for both negative and positive Q values. For the constant-p CPC-1 simulation, Fig. 16d shows that the DQ-Q relationships are close to each other for all instances except that DQ of negative Q values at peak state is notably larger than at other instances and DQ of positive Q values at initial state is slightly larger than at other instances. This agrees with the distribution of local stress fluctuation which as shown in Fig. 17c is similar between different characteristic loading instances. There are more blue to light blue areas at peak state in comparison with other instances, and the distribution of Pi is the most homogeneous at initial state. 6. Conclusions

Fig. 14. Variation of effective stress variance Ve with q for the undrained simulations.

This paper performed a tensor-based analysis on stress variability of granular media in DEM simulations of various loading conditions that are often adopted in soil element testing, including the conventional drained (constant confining pressure), undrained (constant sample volume) and constant mean effective stress simulations. During each simulation, the local stress fluctuation was quantified based on the Euclidean distance between a local stress tensor and the bulk mean

X. Huang et al. / Powder Technology 356 (2019) 581–593

(a)

591

(b) Fig. 15. Variation of Ve at critical state with: (a) q; (b) e.

stress tensor, while the overall stress dispersion was assessed by the effective variance of the stress field within the DEM samples. Multi-fractal analysis was conducted to further characterize the spatial stress variance. Unlike conventional methodologies for assessing the heterogeneity of force transmission network which separately consider the magnitude and orientation of contact forces, the current approach honors the tensorial nature of stresses. The key observations are summarized below:

a) A critical state was reached in all simulations regardless of the initial conditions and loading paths. In particular, a unique q-p’-e relationship at critical state was identified, in agreement with experimental observations. b) The probability density functions (PDFs) of local stress fluctuation could be fitted by lognormal functions whose geometric mean and variance varied during the simulations. Both the geometric mean and variance were influenced by many factors including the loading paths, stress and packing conditions. In particular, in the undrained

(a) Box division within the analysis region

(c) CVC-5

(b) CDC-1

(d) CPC-1

Fig. 16. Multifractal analysis on the local stress fluctuation within a determined region (Peak corresponds to the peak state, IS denotes the instability state, PT indicates the phase transformation state and CS means the critical state).

592

X. Huang et al. / Powder Technology 356 (2019) 581–593

Fig. 17. X-Z projections of local stress fluctuation contour scaled by Pi values within the fixed analysis domain of Fig. 15a at characteristic stages of representative simulations. The contact force network is overlaid, in which the line connects centers of contacting particles and the line thickness is proportional to the magnitude of contact normal forces. (Peak corresponds to the peak state, IS denotes the instability state, PT indicates the phase transformation state and CS means the critical state).

simulations, the geometric mean decreased linearly with decreasing q before phase transformation state and increased linearly with increasing q after phase transformation state. The local stress fluctuation was dominated by particles comprising the major force transmission network. c) Despite its variability during shearing, the geometric mean of local stress fluctuation and its standard deviation became approximately constant at critical state. The critical-state geometric mean of local stress fluctuation increased linearly with increasing deviatoric stress but decreased linearly as void ratio increased. In contrast, the standard deviation of local stress fluctuation at critical state was uniquely correlated with deviatoric stress and void ratio by powerlaw functions with a deviation term. d) The relative orientations between the local stress tensors and bulk stress tensor evolved during shearing but became approximately

constant at critical state. The coaxiality between the local stress tensors and the bulk stress tensor was high throughout the simulations. e) The overall stress variability also showed strong dependencies on loading paths and initial conditions. Similar to the geometric mean of local stress fluctuation, the overall stress variability also varied linearly with q in the undrained simulations. Regardless of the loading paths and initial conditions, the overall stress variability at critical state had a unique power-law relationship with deviatoric stress, while it was correlated with void ratio by a unique power-law function with a deviation term. f) Stress distribution within the DEM samples was highly heterogeneous. Even at initial state in which the bulk stress was isotropic, the stress dispersion was not uniform but also localized. Multifractal analysis showed that in some loading instances, local stress variance was dominated by subboxes where the local stress tensor

X. Huang et al. / Powder Technology 356 (2019) 581–593

deviated slightly from the bulk stress tensor, while in other loading instances, local stress variance was more influenced by subboxes where the local stress tensor deviated notably from the bulk stress tensor.

To conclude, the current research demonstrated the heterogeneous nature of stress variation within a granular material through a set of tensor-based analyzes. In additional to the unique macro stress-strain relationships, the stress variability and its localization characteristics also became invariant at critical state and were uniquely correlated with the macro stress and packing density. The aforementioned observations may provide deep insights into rheology of stress transmission within granular media subjected to external loading. It should be acknowledged that the current study only considered ideally spherical particles which are different from real soil grains. Future efforts will be devoted to investigating the stress variability within non-spherical granular assemblies. Declarations of Competing Interest None. Acknowledgement The research was funded by the National Natural Science Foundation of China (No. 51509186 and 41877227) and the National Key Research and Development Program of China (2017YFC0806004). References [1] L. Papadopoulos, J. Puckett, K.E. Daniels, D.S. Bassett, Evolution of network architecture in a granular material under compression, Phys. Rev. E 94 (3–1) (2016), 032908. [2] T.S. Majmudar, R.P. Behringer, Contact force measurements and stress-induced anisotropy in granular materials, Nature 435 (1079) (2005) 1079–1082. [3] D.S. Bassett, E.T. Owens, M.A. Porter, M.L. Manning, K.E. Daniels, Extraction of forcechain network architecture in granular materials using community detection, Soft Matter 11 (14) (2015) 2731–2744. [4] F. Radjai, M. Jean, J.J. Moreau, S. Roux, Force distributions in dense two-dimensional granular systems, Phys. Rev. Lett. 77 (2) (1996) 274–277. [5] C. Voivret, F. Radjaï, J.Y. Delenne, M.S.E. Youssoufi, Multiscale force networks in highly polydisperse granular media, Phys. Rev. Lett. 102 (17) (2009), 178001. [6] X. Huang, Catherine O'Sullivan, K.J. Hanley, C.Y. Kwok, Partition of the contact force network obtained in discrete element simulations of element tests, Comput. Part. Mech. 4 (2) (2017) 145–152. [7] C. O'Sullivan, L. Cui, Micromechanics of granular material response during load reversals: combined DEM and experimental study, Powder Technol. 193 (3) (2009) 289–302. [8] Miroslav Kramár, A. Goullet, L. Kondic, K. Mischaikow, Evolution of force networks in dense particulate media, Phys. Rev. E 90 (5) (2014), 052203.

593

[9] A. Tordesillas, Force chain buckling, unjamming transitions and shear banding in dense granular assemblies, Philos. Mag. 87 (32) (2007) 4987–5016. [10] J. Zhao, T. Shan, Coupled CFD–DEM simulation of fluid–particle interaction in geomechanics, Powder Technol. 239 (2013) 248–258. [11] R.C. Hidalgo, C.U. Grosse, F. Kun, H.W. Reinhardt, H.J. Herrmann, Evolution of percolating force chains in compressed granular media, Phys. Rev. Lett. 89 (20) (2002), 205501. [12] E.T. Owens, K.E. Daniels, Sound propagation and force chains in granular materials, Eur. Phys. Lett. 94 (5) (2012) 138–161. [13] L. Rothenburg, R. Bathurst, Analytical study of induced anisotropy in idealized granular materials, Géotechnique 39 (1989) 601–614. [14] 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 (3) (2014) 317–333. [15] N. Guo, J. Zhao, The signature of shear-induced anisotropy in granular media, Comput. Geotech. 47 (2013) 1), 1–15. [16] J.D. Zhao, N. Guo, Unique critical state characteristics in granular media considering fabric anisotropy, Géotechnique 63 (8) (2013) 695–704. [17] W. Zhou, J. Liu, G. Ma, X. Chang, Three-dimensional dem investigation of critical state and dilatancy behaviors of granular materials, Acta Geotech. 12 (3) (2017) 527–540. [18] X.M. Xu, D.S. Ling, Y.P. Cheng, Y.M. Chen, Correlation between liquefaction resistance and shear wave velocity of granular soils: a micromechanical perspective, Géotechnique 65 (5) (2015) 337–348. [19] X.Q. Gu, M.S. Huang, J.G. Qian, DEM investigation on the evolution of microstructure in granular soils under shearing, Granular Matter 16 (2014) 91–106. [20] K. Bagi, Stress and strain in granular assemblies, Mech. Mater. 22 (1996) 165–177. [21] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock Mech. Min. Sci. 41 (2004) 1329–1364. [22] K. Gao, J.P. Harrison, Mean and dispersion of stress tensors using Euclidean and Riemannian approaches, Int. J. Rock Mech. Min. Sci. 85 (2016) 165–173. [23] K. Gao, J.P. Harrison, Generation of random stress tensors, Int. J. Rock Mech. Min. Sci. 94 (2017) 18–26. [24] Q.H. Lei, G. Gao, Correlation between fracture network properties and stress variability in geological media, Geophys. Res. Lett. 45 (2018) 3994–4006. [25] Q.H. Lei, G. Gao, A numerical study of stress variability in heterogeneous fractured rocks, Int. J. Rock Mech. Min. Sci. 113 (2019) 121–133. [26] A. Schofield, P. Wroth, Critical State Soil Mechanics, McGraw-Hill, London, 1968. [27] J.E. Gentle, Matrix Algebra: Theory, Computations, and Applications in Statistics, Springer, New York, 2007. [28] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1) (1995) 1–19. [29] G. Simmons, W.F. Brace, Comparison of static and dynamic measurements of compressibility of rocks, J. Geophys. Res. 70 (1965) 5649–5656. [30] X. Huang, K. Hanley, C. O'Sullivan, C.Y. Kwok, Exploring the influence of interparticle friction on critical state behaviour using DEM, Int. J. Numer. Anal. Methods Geomech. 38 (2014) 1276–1297. [31] J. Yang, Non-uniqueness of flow liquefaction line for loose sand, Géotechnique 52 (10) (2002) 757–760. [32] T.G. Murthy, D. Loukidis, J.A.H. Carraro, M. Prezzi, R. Salgado, Undrained monotonic response of clean and silty sands, Géotechnique 57 (3) (2007) 273–288. [33] X.S. Li, Y. Wang, Linear representation of steady-state line for sand, J. Geotech. Geoenviron. Eng. 124 (12) (1998) 1215–1217. [34] X. Huang, K.J. Hanley, Z.X. Zhang, C.Y. Kwok, M.Z. Xu, Jamming analysis on the behaviours of liquefied sand and virgin sand subject to monotonic undrained shearing, Comput. Geotech. 111 (2019) 112–125. [35] P.A. Cowie, D. Sornette, C. Vanneste, Multifractal scaling properties of a growing fault population, Geophys. J. Int. 122 (1995) 457–469. [36] X. Wang, Q. Lei, L. Lonergan, H. Jourde, O. Gosselin, J. Cosgrove, Heterogeneous fluid flow in fractured layered carbonates and its implication for generation of incipient karst, Adv. Water Resour. 107 (2017) 502–516.