Partition of the contact force network obtained in discrete element simulations of element tests

The transmission of stress within a granular material composed of rigid spheres is explored using the discrete element method. The contribution of contacts to both deviatoric stress and structural anisotropy is investigated. The influences of five factors are considered: inter-particle friction coefficient, loading regime, packing density, contact model, and boundary conditions. The data generated indicate that using the above-average normal contact force criterion to decompose the contact force network into two subsets with distinct contributions to stress transmission and structural anisotropy is not robust. The characteristic normal contact forces marking the transition from negative to positive contribution to the overall deviatoric stress and structural anisotropy are not unique values but vary during shearing. Once the critical state is attained (i.e., once shearing continues at a constant deviator stress and solid fraction), the characteristic normal contact force remains approximately constant and this critical state characteristic normal force is observed to decrease with increasing inter-particle friction. The characteristic normal contact force considering the contribution to deviatoric stress has a power-law relationship with the mean effective stress at the critical state.


Introduction
Forces in granular materials are transmitted through interparticle contacts. Physical experiments and numerical studies have revealed that the contact force network is highly heterogeneous [1][2][3]. Under non-isotropic stress conditions, the contact force network exhibits strong geometrical and mechanical anisotropy [4][5][6]. The heterogeneity of the contact force network underlies the complexity of the mechanical behavior of granular assemblies. Therefore, characterization of the contact force network is of crucial importance to understand the behavior of granular media. Radjai et al. [4] proposed that the contacts carrying above-average normal forces and the contacts bearing below-average normal forces have distinct roles in force transmission and work dissipation. Contacts bearing above-average normal forces form a 'strong' force network which transmits the majority of the applied loads. The remaining contacts form a 'weak' contact force network which stabilizes the 'strong' force network and accounts for most of the frictional energy dissipation in the system. These two sub-networks, which collectively comprise the entire contact force network, have different statistical distribution features [2]. This approach to partitioning the contact force network based on average normal force has been adopted in a number of subsequent research studies [5,[7][8][9][10]. Despite its widespread adoption, there is evidence to suggest that partitioning the contact force network in this manner may not be justified. Peters et al. [11] showed that only approximately half of the 'strong' contacts belong to the force transmission chains, and Huang et al. [12] found that the contribution of the 'strong' contacts to frictional dissipation is considerable when inter-particle friction is high. These findings indicate that the use of the above-average criterion to partition contact force networks requires further examination.
This study uses three-dimensional discrete element modeling [13] to study the force transmission mechanism within a granular medium. The contact force characteristics are investigated following the approach of Radjai et al. [4] considering the contribution to deviatoric stress and structural anisotropy. The applicability of the above-average criterion, which was initially proposed based on the results of 2D granular simulations, to more realistic 3D scenarios is explored.

DEM simulations
Most of the simulations were run using a modified version of the LAMMPS code [14]. These three-dimensional simulations contained 20,164 idealized spherical particles. As Fig. 1a shows, the particle size distribution (PSD) of Toyoura sand was approximated with the particle diameters lying between 0.115 mm and 0.408 mm [12]. Thus the size span is around 0.56, where d max and d min are the largest and smallest particle diameters, respectively. This is indicative of a moderate polydisperse system [3]. The samples were enclosed within a cuboidal periodic cell (Fig. 1b) and were initially isotropically compressed until the target stress state and solid fraction had been reached. They were then subjected to triaxial shearing. A simplified Hertz-Mindlin (SHM) contact model was used with a shear modulus of 29 GPa and a Poisson's ratio of 0.12: parameters which approximate the properties of quartz sand [15]. Some periodic-cell simulations were performed using a linear elastic contact model with equal normal and tangential contact stiffnesses of 10 5 N/m. An additional set of simulations used the same PSD and the SHM contact model, but 16,024 spheres enclosed within cylindrical rigid-wall boundaries and these simulations were run using PFC3D [16]. Gravitational body forces were neglected in all simulations. Loading was achieved by moving the upper boundary while the position of the lower boundary was fixed. A constant shear rate, v/L 0 , of 1 s −1 was used for all simulations where v is the velocity of moving boundary and L 0 is the corresponding initial sample vertical dimension. The inertia number was below 10 −4 throughout all the simulations, indicative of quasi-static flow [17,18]. Three loading scenarios were applied: • type-a The sample was sheared while the lateral stresses were maintained constant; In all the three loading regimes, the lateral stresses were kept approximately identical, i.e., σ 2 = σ 3 . Therefore, the deviatoric stress can be simplified as q = σ 1 − σ 3 and the mean stress p equals to 1/3(σ 1 + 2σ 3 ). Five different inter-particle friction coefficients (μ) were used for the type-a loading scenario to investigate the influence of μ on the stress transmission characteristics, while only μ = 0.25 was applied to the type-b and type-c simulations. In all cases shearing continued to a sufficiently high strain level that critical state conditions were attained, i.e., the combination of stress level and packing density beyond which shearing continues at a constant deviator stress, constant mean stress, and constant packing density [19,20].

Results and discussion
The stress and deformation responses of a representative subset of the simulations subject to the type-a loading condition using periodic boundaries and the SHM contact model are illustrated in Fig. 2. All these simulations considered the same initial isotropic conditions prior to shearing (initial solid fraction, s, of 0.652 and an initial isotropic mean stress, p 0 , of 100 kPa) and the confining pressure was maintained at 100 kPa during shearing. Referring to Fig. 2a, for all samples the deviatoric stress q initially increases rapidly to a peak value and subsequently decreases to become almost constant when the axial strain ε a exceeds 20 %. The deviatoric stress initially increases with increasing μ but becomes insensitive to

Fig. 2
Effect of μ on the evolution of q (a) and solid fraction (b) during simulations using periodic boundaries and the type-a loading protocol μ when μ exceeds 0.5 due to the interplay between rotational behavior and sliding behavior at the contacts [12]. Figure 2b shows that all samples dilate from their initially high density and attain nearly constant solid fractions that decrease with increasing μ. The constant deviatoric stress and solid fraction at large ε a levels represent the typical critical state (steady-state) response which is the ultimate state that a soil can approach [19,20].
The average stress within a granular assembly can be calculated directly from the contact forces [21] in which V is the volume of the domain considered, f i is the ith component of the contact force vector f , l j is the jth component of the branch vector l that joins the centers of two touching particles and N c is the number of contacts. The particle velocities and particle inertias are sufficiently small that the dynamic component of the stress tensor (as described in Claudin [22]) can be neglected. Following the approach of Radjai et al. [4], Fig. 3 presents the cumulative contribution of contact forces to the overall deviatoric stress as a function of f n / f n for the type-a simulations considered in Fig. 2 at ε a = 30 % (at this axial strain level the critical state has been attained in all the simulations); f n is the mean normal contact force over the entire assembly and the corresponding q value for each f n / f n considers the contribution of the subset of contacts with f n / f n below the current value. To create Fig. 3, the contacts are firstly sorted in an ascending order according to the magnitude of f n / f n (i.e., the ξ parameter in Radjai et al. [4]). The contribution of each contact to q is then calculated based on Eq. 1. Each data point in Fig. 3 represents the cumulative sum of the contributions to q considering contacts with f n / f n below the current value. There is a range of f n / f n values where the cumulative curves give q < 0. When the stress tensor is calculated considering only contacts with force values lying in this range, the resultant principal stress orientation will differ from the overall principal stress orientation by more than 45 • . An enlarged view of the range 0 < f n / f n < 1 is superimposed on the figure to enable a closer examination of the transition from negative to positive contributions to q. Contacts that make a positive contribution to q are on average considered part of the force transmission network following Radjai et al. [4]. The data show that at the strain level considered the characteristic normal contact force f * which separates the positive and negative contributions to q is smaller than unity for all the μ values considered. Figure 4 presents the projection of the normal contact force onto the X(σ 3 ) − Z(σ 1 ) plane at ε a = 30 %.
For clarity, only contacts within the middle part of the samples (10 % of the entire sample) are considered. The contacts are represented by lines connecting the centers of the particles in contact, and the thickness of each line is proportional to the magnitude of f n / f n . The overall contact force network is divided into two subsets according to the magnitude of f n / f n with respect to f * , i.e., f n / f n ≥ f * (on the left in black) and f n / f n < f * (on the right in gray). For the μ values considered, the configurations of contact force networks are quite similar. The contacts within the f n / f n ≥ f * subset are oriented along the major principal stress direction, while the orientations of the contacts belonging to the f n / f n < f * have a random orientation.This reflects the distinct roles of the two subsets: the f n / f n ≥ f * subset acts as force-bearing backbones carrying the majority of deviatoric loading, while the f n / f n < f * subset is the supporting network which provides lateral support to the force-bearing network. Figure 5 shows the evolution of f * for q during shearing for the simulations presented in Fig. 2. For all the μ values considered, f * for q is relatively high at the isotropic state and varies during subsequent shearing, tending to nearly a constant value at the critical state (ε a ≥ 30 %). Moreover, at all strain levels f * increases with decreasing μ. The structural anisotropy also increases as μ increases [12]. The non-uniqueness of f * for q is also evident for the type-b (Fig. 6a) and the type-c (Fig. 6b) loading scenarios, for the type-a loading using the linear elastic model (Fig.  6c) and for the type-a loading using the rigid-wall boundaries (Fig. 6d). Figures 3, 4, 5, and 6 indicate that many of the contacts with 0 < f n / f n < 1 which have been c for a representative simulation using the type-a loading protocol and the linear elastic contact model; d for a representative simulation using rigid-wall boundaries and the type-a loading protocol hitherto considered 'weak' are not only contributing to the secondary supporting force network but also to the deviatoric load and so their classification as "weak" may not be appropriate. The contact anisotropy is usually quantified using the fabric tensor [23]: where n k i denotes the unit contact normal. The degree of contact anisotropy is evaluated using deviatoric fabric defined , where Φ 1 , Φ 2 , and Φ 3 are the major, intermediate, and minor principal fabrics, respectively, and these are the eigenvalues for the fabric tensor. It was found that Φ 2 and Φ 3 are effectively identical in the simulations as a σ 2 = σ 3 condition was imposed. Therefore, the deviatoric fabric can be simplified as Figure 7 shows Φ d as a function of f n / f n at the critical state for each μ value considered. To create Fig. 7, the contacts are firstly sorted in an ascending order according to the magnitude of f n / f n , as was required for Fig. 3. The contribution of each contact to Φ d is then calculated based on Eq. 2. Each data point in Fig. 7 represents the cumulative sum of the contribution to Φ d considering the contacts with f n / f n below the current value. Analogous to the stress calculations described above, there is a range of f n / f n values that give Φ d < 0 on the cumulative curves, and the fabric tensor calculated considering only contacts with force values lying in this range gives a principal fabric orientation that differs from the overall principal fabric orientation by more than 45 • . It can be seen from Fig. 7 that Φ d increases as μ increases. This is because with increasing μ, the inherent stability of the force-bearing skeleton increases,  Fig. 9 Variation of f * at the critical state with μ for representative simulations using the type-a loading protocol the lateral supporting contacts orthogonal to the strong force chains are not activated, and hence the anisotropy increases. The inset in Fig. 7 shows that for all μ values considered, the f * value which marks the transition from a negative contribution to a positive contribution to Φ d generally decreases with increasing μ value. Figure 8 presents the variation of f * considering the contribution to Φ d with axial strain (ε a ) during shearing for a representative set of the type-a simulations. It is clear that just as f * for q varied, f * for Φ d is also not a constant; it depends on μ and the strain level. Similar to f * for q, f * for Φ d also remains approximately constant when the axial strain ε a exceeds 30 % and critical state shearing conditions are attained. Figure 9 illustrates the variation of f * at the critical state with μ considering both contributions to q and Φ d . f * values defined for q differ from those defined for Φ d ; the former are smaller at small μ values (< 0.75) but the two become almost identical at μ values ≥0.75. Similar observations can also be made for other types of loading conditions for which data are not shown herein for conciseness. This is because apart from geometrical anisotropy Φ d , q is  Fig. 10 Correlations between the characteristic force f * with mean stress p at the critical state: a critical state f * for Φ d with p; b variation of power-law regression parameters with μ; c critical state f * for q with p also a function of mechanical anisotropies [4,24]. Both f * for q and f * for Φ d decrease with increasing μ.
As shown in Figures 5, 6, and 8, just as is the case for the stress and volumetric responses at large strain levels, both f * for q and f * for Φ d become almost constant at the critical state. Figure 10a shows the relationship between the critical state f * for Φ d and the mean stress p at the critical state. Due to the presence of oscillations, the critical state data presented in Fig. 10a are taken as the mean values of the corresponding quantities for the last 10 % of axial strain. f * for Φ d at the critical state increases as p increases but decreases with increasing μ, in line with Figs. 7 and 8. This is the opposite of the correlations between Φ d at the critical state and p and between Φ d and μ noted in [12]. The relationships between f * for Φ d for different μ values can be represented by a series of power-law functions, i.e., in a form of f * = ηp ξ . The fitting parameters η and ξ for each μ case are given in Table  1. η decreases with increasing μ while the opposite trend is evident for ξ . Figure 10b presents the variation of η and ξ with μ. Both of them follow a power-law relationship with μ. η decreases with increasing μ and thus the fitting function has a negative power, while ξ increases with increasing μ and the fitting function has a positive power. Figure 10c plots f * for q at the critical state against p. In comparison to f * for Φ d , the relationship between the critical state f * for q and p is more complicated; for all the μ values considered, it initially decreases with increasing p, reaches a local minimum and then continues to increase as p is further increased.

Conclusions
This study presents a systematic investigation on the force transmission characteristics within granular media. The influences of multiple factors were considered including inter-particle friction coefficient, loading regime, packing density, contact model, and boundary conditions. The evidence presented in this study clearly suggests that some of the weak contacts carrying below-average normal contact forces also participate in stress transmission and contribute to structural anisotropy. It is shown that both the characteristic normal contact force which marks the transition from a negative to a positive contribution to the overall deviatoric stress and that marking the transition from a negative to a positive contribution to structural anisotropy vary until the critical state is reached and remain approximately constant thereafter. The characteristic normal contact force for the deviatoric stress and that for the structural anisotropy are not unique and differ from each other. Thus, using the average normal contact force f n as the transition point to sepa-rate a load-bearing network which is anisotropic and aligned along the external loading direction from a supporting network which is isotropic and orthogonal to the load-bearing network within a granular packing is not robust for threedimensional systems.