## Abstract

The peak and critical state strengths of sands are linearly related to the stress level, just as the frictional resistance to sliding along an interface is related to the normal force. The analogy with frictional sliding has led to the use of a ‘friction angle’ to describe the relationship between strength and stress for soils. The term ‘friction angle’ implies that the underlying mechanism is frictional resistance at the particle contacts. However, experiments and discrete element simulations indicate that the material friction angle is not simply related to the friction angle at the particle contacts. Experiments and particle-scale simulations of model sands have also revealed the presence of strong force chains, aligned with the major principal stress. Buckling of these strong force chains has been proposed as an alternative to the frictional-sliding failure mechanism. Here, using an idealized abstraction of a strong force chain, the resistance is shown to be linearly proportional to the magnitude of the lateral forces supporting the force chain. Considering a triaxial stress state, and drawing an analogy between the lateral forces and the confining pressure in a triaxial test, a linear relationship between stress level and strength is seen to emerge from the failure-by-buckling hypothesis.

## 2. Background

The shear strength of sand is often described using a Mohr–Coulomb linear failure envelope [1,2]. This relates shear strength at failure (*τ*_{ff}) and the effective normal stress on the failure plane (*σ*′_{ff}) as
*ϕ*′ is the angle of shearing resistance, sometimes termed the friction angle. Recall that Coulomb friction relates the normal force and maximum shear force along a surface (interface) as
*T*_{f} is the maximum shear (tangential) force that can be attained along the surface, *N* is the normal force and *ϕ*_{surface} is the friction angle of the interface [3]. The similarity between equations (2.1) and (2.2) gives rise to use of the term ‘friction angle’ in soil mechanics. When triaxial stress conditions are encountered, equation (2.1) is better expressed [4] as
*σ*′_{1f} and *σ*′_{3f} are the major and minor principal stresses at the point of failure, respectively. For sands, two angles of shearing resistance are of particular importance: the critical state (*ϕ*′_{cv}) and the peak (*ϕ*′_{p}).

The critical state soil mechanics (CSSM) framework, proposed by Roscoe *et al*. [5] and documented by Schofield & Wroth [6], relates the deviator stress (*q*) and mean effective stress (*p*′) at the critical state as
*M* is a critical state parameter. For a triaxial stress state, *q*=*σ*′_{1}−*σ*′_{3} and *p*′=(*σ*′_{1}+2*σ*′_{3})/2, where *σ*′_{1} is the major principal stress and *σ*′_{3} is both the minor principal stress and the confining pressure. Equation (2.5) can also be expressed as equation (2.6) for a triaxial stress state:
*σ*′_{1} and *σ*′_{3}.

Soil strength is typically attributed to originate from both frictional resistance at the contacts between grains and the kinematics of the relative motion of grains [7]. The kinematic constraint is usually related to dilation and the initial state parameter; this is related to the difference between *ϕ*′_{p} and *ϕ*′_{cv} [8,9]. Countless laboratory experiments (e.g. triaxial tests of Verdugo & Ishihara [10] and plane strain tests of Wanatowski & Chu [11]) have confirmed that, for a given sand with a fixed mineralogy and particle size distribution and hence a constant *ϕ*′ value, *τ*_{ff} does increase linearly with stress level. However, the evidence to suggest that friction at the particle contacts is the key factor explaining the increase of *τ*_{ff} with increasing *ϕ*′ is not entirely convincing. Mitchell & Soga [7] clearly acknowledged the ample experimental evidence which indicates that *ϕ*′_{cv} is not solely a function of mineralogy and highlighted the influence of both particle size distribution and particle shape. Mitchell and Soga proposed that *ϕ*′_{cv} is determined by a combination of ‘true friction’ and particle rearrangement including fabric development.

The discrete element method (DEM) facilitates parametric studies that have examined the influence of particle–particle contact friction on *ϕ*′_{cv}. The DEM simulations of Thornton [12], Peña *et al.* [13], Yang *et al.* [14] and Huang *et al.* [15], among others, have shown that *ϕ*′_{cv} is sensitive to the inter-particle friction angle (*ϕ*_{surface}) only at low values of *ϕ*_{surface}, as illustrated in figure 1. These numerical results indicate that frictional sliding along the interfaces defined by particle contacts is not the key mechanism determining the critical state strength findings. These numerical studies are in agreement with the earlier experimental observations of Skinner [16], which are also presented in figure 1.

It is now broadly accepted that the stress distribution within a sand is heterogeneous; discrete subnetworks or chains of contacting relatively highly stressed particles form upon loading. These ‘strong force chains’ are aligned with the major principal stress orientation and they have been observed in DEM simulations [17], in photoelastic experiments using analogue sands [18], and in real sands using both image analysis of photographs [19] and micro-computed tomography [20]. Figure 2 illustrates these force chains using data from a two-dimensional DEM simulation of biaxial compression. There are clear hypotheses in the literature linking soil failure to the buckling or collapse of these strong force chains [21]. At the critical state, or in a shear band, it seems that there is continuous formation of new force chains concurrent with failure of pre-existing force chains by buckling [22].

Mitchell & Soga [7] explicitly attribute the non-proportional relationship between *ϕ*_{surface} and *ϕ*′_{cv} illustrated in figure 1 to the formation of the strong force network, and state that friction (i.e. *ϕ*_{surface}) acts to stabilize this network but is not ‘the direct source of macroscopic resistance to shear’. Mitchell and Soga's statement is really speculation that force chain buckling may provide a valid conceptual framework for soil failure; proof of this hypothesis requires a rational explanation for the linear stress–strength relationship expressed in equations (2.1), (2.4) and (2.6). Here, following the earlier contributions of Hunt *et al.* [23] and O'Sullivan *et al.* [24], several simple models of isolated force chains are used to examine the link between buckling and stress-dependent strength. In each of these models, the particles are represented by nodes that are connected by rigid links to model the strong force chain and lateral forces are applied to represent the weak force network orthogonal to *σ*′_{1} that transmits the confining pressure (i.e. *σ*′_{3}). The relationship between the lateral forces, which simulate *σ*′_{3}, and the post-buckling load following collapse of the force chain, which is analogous to *σ*′_{1}, is assessed. The main difference between this work and the earlier work of Hunt *et al.* [23] is the development of an explicit relationship between the principal stresses within a real sample subjected to triaxial loading and the output of these analogue models. In the prior contribution of O'Sullivan *et al.* [24], data obtained from DEM simulations of true triaxial tests were used as input for a spring-and-link model to show that the relative support offered to strong force chains in the minor and intermediate principal stress directions explains the observed dependency of sample strength on the intermediate stress ratio,

First, a two-dimensional model containing a single node is developed (§3). This model is sufficiently simple to permit an analytical solution for the post-critical equilibrium path to be obtained in closed form without resorting to numerical bifurcation analysis. This model is then extended into three dimensions and multiple nodes (§4).

## 3. Development of a two-dimensional single-node model

This study neglects explicit consideration of the supporting role played by inter-particle friction to increase the stability of strong force chains and considers only the contribution of the lateral supports. Previously, Tordesillas *et al.* [22] acknowledged that lateral supporting contacts play a key role in supporting the force chains. Using DEM, Barreto & O'Sullivan [25] showed that both inter-particle friction and the weaker network of contacts oriented orthogonal to the major principal stress direction contribute to the stability of these force chains.

The analogue models used in this paper have a similar form to the three-dimensional model developed by O'Sullivan *et al.* [24], which was based on a two-dimensional model proposed by Hunt *et al.* [23]. Consider the very simple abstraction of a force chain shown in figure 3 which contains one node and two rigid links of equal length *L*. Two springs, a lateral spring of stiffness *k*_{f} and a rotational spring of stiffness *k*_{r}, resist motion of the node from its equilibrium position at which both springs are unstressed. The lateral spring represents the lateral force chains or the weak contact network oriented in the direction of *σ*′_{3} that acts to support the strong force chains. The rotational spring represents the rotational resistance present at the contacts owing to particle geometry effects and friction.

O'Sullivan *et al.* [24] considered a true triaxial stress state when developing a three-dimensional spring-and-link model, i.e. *σ*′_{3}≤*σ*′_{2}≤*σ*′_{1}, and the relative values of the spring stiffnesses were chosen by considering the axial stiffnesses in the *σ*′_{3} and *σ*′_{2} directions measured in DEM simulations of true triaxial tests [26]. In this paper, the earlier work of O'Sullivan *et al.* has been advanced, so that a constant lateral confining force (*F*_{L}) can be applied to mimic a constant value of *σ*′_{3}, i.e. the lateral support is force-controlled rather than stiffness-controlled.

Two different types of lateral spring were used: a linear spring for which the force is directly proportional to deflection (i.e. *F*_{L}=*k*_{f}*δ*) and a nonlinear Hertzian spring for which *F*_{L}=*K*_{f}*δ*^{3/2}. Here, *k*_{f} corresponds to the stiffness of the linear spring with SI units of N m^{−1}, whereas *K*_{f} is the stiffness coefficient of the Hertzian spring (N m^{−3/2}). For both types of spring, the force increases until a deflection *δ*_{0} is attained; thereafter, the force becomes constant and independent of compressive deflection of the spring, *δ*, as shown in figure 4. It is necessary to have a piecewise description to avoid singularities at the critical state.

Referring to figure 3, the bottom link is pin-jointed, and a load *P* is applied to the uppermost link, increasing from zero until buckling occurs. A total potential energy function, *V* , is developed for the system, which has the same form as that described by O'Sullivan *et al.* [24], i.e.
*U*_{L} and *U*_{L} represent the total elastic potential energy of the lateral and rotational springs, respectively, and *P*Δ is the work done by the load *P*. The *U*_{L} term was developed by dividing the response into two parts. The strain energy stored in the linear lateral spring, i.e. the area enclosed by figure 4*a*, is
*b*, are

Of course, this assumes that *δ*>0 during the loading history. Following the approach of O'Sullivan *et al.* [24], *P*Δ and *U*_{L} can be expressed as
*U*_{L} and *δ*=*qL*. At equilibrium, the potential energy is stationary, i.e. energy dissipation due to friction or particle breakage is neglected and so d*V*/d*q*=0. Thus, solving for *P* for the four distinct cases discussed gives the following relationships

Any load beyond the transition point (which occurs at *δ*=*δ*_{0} in this simple model) can be identified as a suitable resistance because it is beyond that point at which the constant force confinement becomes valid. This post-buckling load represents the vertical load that maintains equilibrium; for a given displacement, *δ*, any load greater than this resistance would cause collapse. Note that it is essential to compare the model at the same deflection. For the single-node model, this was arbitrarily chosen as *q*=0.003, because this value is always higher than *δ*_{0}/*L*.

## 4. Extension to a three-dimensional multiple-node model

The two-dimensional single-node model developed in §3 can be extended to three dimensions and multiple nodes. A vertical force chain is idealized as *N* nodes connected by *N*+1 rigid links, each of length *L*, as shown in figure 5 for the *N*=3 case. The displacements of any node *i* in the *x* and *y* directions from the undeflected state with unstressed springs are given by *q*_{ix}*L* and *q*_{iy}*L*, respectively. Each node is supported by four springs: two linear lateral springs of equal stiffness (*k*_{f}) and two rotational springs.

One of the rotational spring stiffnesses, *k*_{ry}, was specified by fixing the ratio *k*_{ry}:*k*_{rx} at 1.1 to avoid numerical difficulties with concurrent zero eigenvalues arising during the matrix inversion required for the solution that is outlined below. Two approaches were adopted to estimate suitable values of *k*_{rx}. One approach was to set the rotational stiffnesses at arbitrarily chosen constant values. Alternatively, the stiffnesses were taken to be proportional to the load *P*; this is comparable to some of the rolling resistance models proposed in the literature where the rolling stiffness is proportional to the normal contact force [22,27–29].

The relevant equations to describe both *U*_{L}, the total elastic potential energy of the rotational springs, and *P*Δ, the work done by the load *P*, were developed by O'Sullivan *et al.* [24]:

For the three-dimensional multiple-node model, linear lateral springs were used which have the force–deflection characteristics shown in figure 4*a*. For this multiple-node model, negative deflections are permitted for which a symmetric relationship between *F*_{L} and *δ* is assumed, i.e. *F*_{L}=−*k*_{f}|*δ*_{0}| if *δ*≤−|*δ*_{0}|. The deflections of the two lateral springs are
*U*_{L} term was developed by dividing the response into two parts for each node *i*. When *δ*_{ix} or *δ*_{iy}≤*δ*_{0}, the strain energy stored in the lateral spring is *et al.* [24]. When *δ*_{ix} or *δ*_{iy}>*δ*_{0}, the respective strain energy terms become
*N* linear springs is obtained by summing the energy terms for both lateral springs over all *N* nodes.

The expression for the total potential energy function, *V* , contains five constants: *k*_{f}, *k*_{rx}, *k*_{ry}, *L* and *δ*_{0}; *V* can be non-dimensionalized as *k*_{f}*L*^{2} to leave three dimensionless groups of constants:
*i*=1,2,…,*N*. These partial derivatives were calculated symbolically using Maple 17.00 [30]. The CodeGeneration package within Maple was used to generate a Fortran equations file suitable for input to AUTO-07p [31] following the approach of O'Sullivan *et al.* [24].

For each model evaluation, non-dimensional resistances (*q*=0.003 for the two-dimensional single-node model. *σ*′_{1}. Thus, any pattern in the variation of *F*_{L} is analogous to a variation in *σ*′_{1} with *σ*′_{3}. In the majority of the simulations, the initial displacements of all nodes were zero (i.e. *q*_{xi}=*q*_{yi}=0 for all *i*). As real materials contain flaws, an imperfection was included in a subset of the simulations by giving the middle node an initial perturbation of equal size in both the *x* and *y* directions. The simulations were run using two different numbers of nodes: three and seven. The confining force was varied by using a range of values for

## 5. Results and discussion

### 5.1 Two-dimensional single-node model

Four separate cases were considered for the two-dimensional model: both linear and Hertzian lateral springs were used, and for each spring type, the confining force *F*_{L} was varied by adjusting either the spring cut-off, *δ*_{0}=*q*_{0}*L*, or the spring stiffness parameter (i.e. *k*_{f} or *K*_{f} for linear or Hertzian springs, respectively). It is noted that the confining force at *δ*>*δ*_{0} is linearly proportional to the spring stiffness parameter for both linear and Hertzian lateral springs (i.e. *F*_{L}=*k*_{f}*δ* or *K*_{f}*δ*^{3/2}), whereas *F*_{L} is nonlinearly related to deflection for the Hertzian model.

Figure 6 shows model evaluations for the four cases considered. For those two cases in which *δ*_{0} was varied, *q*_{0} values of 3×10^{−4}, 6×10^{−4} or 9×10^{−4} were used with a fixed *k*_{f} of 200 N m^{−1} (linear springs) or 200 N m^{−3/2} (Hertzian springs). In the other two cases, a fixed *q*_{0} of 6×10^{−4} was used with *k*_{f} or *K*_{f} values of 100, 200 or 300 (N m^{−1} or N m^{−3/2}). Regardless of spring type or the level of the confining force, there is a clear transition point at *q*=*q*_{0} characterized by a sharp change of slope beyond which the load *P* decreases monotonically with increasing deflection. The loads may be compared at any arbitrarily chosen deflection along these decreasing paths; the small inset subfigures plot load against confining force at *q*=0.003. The load is linearly related to confining force even if the initial lateral spring behaviour is nonlinear (i.e. Hertzian) and this linear relationship holds for all *q* values that exceed *q*_{0}. The latter point is important as soil is inherently a nonlinear material. Because the column resistance and confining force in this simple model represent *σ*′_{1} and *σ*′_{3}, a linear relationship emerges between *σ*′_{1} and *σ*′_{3}, irrespective of the manner in which the confining force is applied, i.e. equations (2.4) and (2.6) emerge from this simple abstraction.

### 5.2 Three-dimensional multiple-node model

The two-dimensional single-node model is convenient for analysis because of its simplicity. However, real force chains exist in three dimensions and can comprise many contacting particles. For this reason, it is necessary to confirm that the linear relationship between confining force and resistance observed for the two-dimensional single-node model is also true for the three-dimensional multiple-node model developed in §4.

Figure 7*a* shows that the column resistance increases linearly with confining force for a three-node chain using two definitions of *α* and *β* are arbitrarily chosen constants. Figure 7*b* is similar to figure 7*a* except that a chain length of seven nodes is considered. Because the aim of this study is to give semi-qualitative trends, the non-dimensional confining and post-buckling forces for each set of multiple-node simulations were normalized by their respective maximum values.

The critical buckling force increases linearly with confining force for both chain lengths, regardless of whether *b* are due to the inability of AUTO to find the lowest eigenvalue in a limited number of cases owing to close proximity with other eigenmodes causing numerical problems, as discussed by Wadee *et al.* [32]. Because the post-buckling force is taken to be analogous to *σ*′_{1} and the confining force is analogous to *σ*′_{3}, the trends expressed in equations (2.4) and (2.6) are again captured by this model which does not contain an inter-particle friction parameter. It is noted that when the rotational stiffness is independent of confining force, the rotational springs are sufficient to prevent buckling at negligibly small forces. Hence,

Figure 7 is for a column in which the nodes are initially collinear (vertical). In a real material, the particles comprising the force chains are not perfectly aligned (as illustrated in figure 2). The three-node model was modified by giving the middle node a small initial perturbation. Figure 8*a* shows the sharp reduction of *b* shows that when the force chain is initially perturbed, in this case by 0.00001% of link length, the linear relationship between resistance and confining force observed in figure 7 still fundamentally exists.

## 6. Conclusion

The dependency of soil strength upon stress level is well known and it is commonly described as a frictional relationship. Prior experiments and DEM simulations have shown that *ϕ*′_{p} and *ϕ*′_{cv}, which are often referred to as friction angles, are not simply related to particle surface friction, *ϕ*_{surface}, suggesting that a framework which considers the critical state strength to be a purely frictional strength has a tenuous scientific basis. There is by now ample evidence to indicate that, when subject to a non-isotropic stress state, sand particles and their contacts are oriented in the direction of *σ*′_{1} to form columns of preferentially stressed particles. Abstractions of these columns can be created using simple spring-and-link models. By using these abstract models and controlling the applied lateral forces, it has been shown that the actual post-buckling strength of these columns is approximately linearly related to the lateral supporting force when the confining force is basically constant. This model, therefore, indicates that the linear relationships between *τ*_{ff} and *σ*′_{f} and between *q* and *p*′ emerge from the fundamental mechanics of force chain buckling. The model explains the existence of these relationships despite the absence of an explicit link between inter-particle friction (*ϕ*_{surface}) and the critical state angle of shearing resistance (*ϕ*′_{cv}).

## Data accessibility

All of the data used to plot figures 7 and 8 are deposited in Dryad. Sets of AUTO input files for these simulations (*.f and c.*) and the associated AUTO output files (fort.7, fort.8, fort.9) are also available at http://dx.doi.org/10.5061/dryad.05478.

## Funding statement

K.J.H. was supported as a post-doctoral scholar by the Royal Commission for the Exhibition of 1851 during the period of this research study.

## Author contributions

C.O.S. conceived of the main idea and wrote most of §§2 and 6. X.H. contributed to §2 and conceived of introducing a perturbation to the model (figure 8). K.J.H. and M.A.W. developed the analytical method and designed the simulations. K.J.H. ran the simulations and wrote most of §§3–5 with the assistance of M.A.W. All authors contributed to the data analysis and reviewed the results. All authors reviewed and approved the final manuscript.

## Conflict of interests

None of the authors have any competing interests.

- Received January 22, 2015.
- Accepted March 23, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.