The authors have presented some interesting experimental studies on strength anisotropy of sand using direct shear tests. It is indeed surprising to see that the minimum peak friction angle is not observed when the sample is sheared along the bedding plane orientation. The authors also mentioned that it is desirable to develop theoretical modelling of the observed soil response. It is emphasized by the authors in this paper and a previous one [1] that it is practically more convenient to formulate the anisotropic failure criterion of sand in terms of the shear plane (or failure plane) orientation, as some existing geotechnical analysis methods which rely on shear strength of soils such as slope stability analysis and foundation failure analysis explicitly assume a potential failure plane without referring to the loading direction. The discusser thinks that a failure criterion formulated this way would be useful for special cases in which the failure planes are known but may not be suitable for general cases in which the failure planes cannot be easily determined. This discussion will show that the observed strength anisotropy of sand can be successfully modelled by a fabric-tensor-based sand model, in the formulations of which the potential failure plane orientation is not required.

1 Finite element simulation of the strength anisotropy of sand in direct shear tests

Direct shear tests should not be treated as single element tests as the stress and strain distribution inside the samples is highly non-uniform. In this discussion, the direct shear test is simulated as a boundary value problem using a newly developed fabric-based sand model [2]. This model employs a void-based fabric tensor and a physically based fabric evolution law to account for the influence of void sizes and orientations and their change during shear on the sand behaviour. Details of the model formulations can be found in Gao et al. [2]. The model has been implemented in Abaqus using the explicit stress integration method [3]. In the simulations presented here, the parameters for Toyoura sand are used (see [2]).

The test setup is shown in Fig. 1. The sample length l is 60 mm and height h is 20 mm. In addition, there is a 1.5-mm-thick sand layer between the upper and lower boxes (Fig. 1). Definition of the initial bedding plane orientation α is shown in Fig. 1. For all the tests simulated here, the same boundary conditions are used and uniform distribution of initial void ratio e 0 is assumed. An initially K 0 stress state is assumed for all the tests, and K 0 = 1 − sin φ c = 0.48 is used, where φ c (= 31.1° for Toyoura sand) is the friction angle at the critical state in triaxial compression. During the tests, uniform horizontal displacement u is applied on the vertical sides of the upper box and constant confining pressure σ n is applied on the top of the sample (Fig. 1a). Inclined horizontal deformation field is applied to the left and right vertical sides of the shear zone (Fig. 1a). The top surface of the sample is assumed to remain horizontal throughout the tests. Rough boundary condition is assumed for the interaction between sand and top cap as well as sand and the bottom of the lower box. Linear four-node plane strain element is used, and Fig. 1b shows the mesh sizes. All the meshes have the same width of 1.2 mm. The maximum and minimum mesh height in the two boxes is 2.5 and 1 mm, respectively. Table 1 shows the summary of the simulated tests. F 0 is the initial degree of anisotropy.

Fig. 1
figure 1

Illustration of a the sample dimension, boundary conditions and initial bedding plane orientation for direct shear tests and b mesh size for direct shear tests

Table 1 Summary of the simulated tests

1.1 General observation on the strength anisotropy

Typical global stress and strain relations are shown in Fig. 2a, in which γ (=u/h) is the global shear strain and τ is the global shear stress measured on the boundary (based on the horizontal reaction force at the bottom of the lower box). Strain softening response is observed for all the tests, and the peak global stress ratio (τ/σ n )max is found to be strongly dependent on α. At large strain level (γ ≈ 11 %), approximately the same residual stress ratio around 0.65 (or equivalently the same global residual friction angle) is observed for all tests.

Fig. 2
figure 2

Simulated a global stress and strain relations, b vertical displacements of the cap (global dilatancy) and c global dilation angle evolution for two samples (F 0 = 0.5)

Figure 3a shows the simulated variation of peak friction angle φ p [=tan−1(τ/σ n )max] with α for Toyoura sand. To facilitate the comparison between the numerical simulations with the experimental results presented by the authors, the data shown in Fig. 6b–d of the paper are reproduced in Fig. 3b–d. It can be seen that the simulated variation of φ p with α is similar to the experimental observations for Fujian sand and glass beads. First, φ p at α = 0° (and α = 180°) is neither the minimum nor the maximum. Secondly, two minimum φ p values are observed for these three materials. Such similarity is probably attributable to that all these three materials have similar particle shapes. The strength anisotropy observed for Mica sand is different from that for the other materials as this sand has much more angular particle shape.

Fig. 3
figure 3

Variation of peak friction angle φ p with initial bedding plane orientation α = for a Toyoura sand with σ n  = 400 kPa (numerical simulation in this study), b Fujian sand with σ n  = 200 kPa, c Mica sand with σ n  = 300 kPa and d glass beads with σ n  = 150 kPa (test data from the paper)

1.2 Shear dilatancy behaviour

The shear dilatancy response of sand samples with different initial bedding plane orientations is shown in Fig. 2b. The vertical displacement v is measured at the top of the upper shear box. Note that positive and negative v indicates volumetric expansion and contraction of the sample, respectively. It is evident that all the samples show volumetric contraction at the initial loading stage and extensive volumetric expansion afterwards. At large strain level of γ ≈ 11 %, the vertical displacement reaches the maximum and ceases to increase.

Figure 2c shows the evolution of global dilation angle Ψ [=tan−1(δv/δh] with the global shear strain for samples with α = 30° and α = 120°, where δv and δh denote the vertical and horizontal displacement increments, respectively. It is evident that the peak global dilation angle Ψ p and friction angle φ p are observed at approximately the same global strain level for both samples (Fig. 2a, c). Similar observations are also found for cases with other values of α. This is indeed in agreement with the experimental observations shown in the paper. Figure 4 shows the variation of global peak dilation angle Ψ p with initial bedding plane orientation α. It can be found that the variation of Ψ p with α is similar to that of φ p with α shown in Fig. 3a, which is also observed in the experiments by the authors (Fig. 16 of the paper).

Fig. 4
figure 4

Variation of peak dilation angle Ψ p with initial bedding plane orientation

2 Conclusion

The strength anisotropy of sand observed in direct shear tests by the authors can be captured by the numerical simulation using a fabric-based sand model. Therefore, it may be useful to develop an anisotropic failure criterion for sand in terms the orientation of failure planes for special cases in which the failure planes are known, the fabric-based model is more general and can be used to describe the strength anisotropy of sand under more general loading conditions.