Dependency of Dilatancy Ratio on Fabric Anisotropy in Granular Materials

The volumetric deformation behavior of granular materials under shear is a fundamental characteristic of such material. This volumetric deformation, either in contraction or expansion, can be referred to in a general sense as dilatancy. Since the establishment of the dilatancy concept for granular materials by Reynolds (1885), significant effort has been dedicated to seeking the relationship between the dilatancy ratio and the stress ratio (e.g., Oda 1975; Bolton 1986; Wood 1990; Collins and Muhunthan 2003), following the pioneering work by Taylor (1948) and Rowe (1962). The dilatancy ratio $D$ can be expressed as the plastic volumetric strain increment $\mathrm{d}{\epsilon }_v^p$ divided by the plastic deviatoric strain increment $\mathrm{d}{\epsilon }_q^p$, while the stress ratio $\eta$ can be expressed in terms of the deviatoric stress $q$ divided by the mean effective stress $p$.

Internal material characteristics such as void ratio and anisotropy are also considered to influence the dilatancy behavior of granular materials (Rowe 1962; Dafalias 1986; Houlsby 1993). The dependency of dilatancy on void ratio has been studied extensively with significant accumulation of experimental data (Been and Jefferies 1985; Li and Dafalias 2000), and its inclusion in constitutive models has become a standard practice (Manzari and Dafalias 1997; Gajo and Wood 1999; Wang et al. 2014). However, a conclusive consensus is yet to be reached regarding the influence of anisotropy on dilatancy ratio. Constitutive theories have been formulated to reflect the dependency of dilatancy on fabric anisotropy in granular materials (Nemat-Nasser 2000; Li and Dafalias 2002; Dafalias et al. 2004; Wan and Guo 2004; Wan et al. 2010; Li and Dafalias 2012; Zhao and Gao 2016). In particular, Li and Dafalias (2012) and Dafalias (2016) directly incorporated a fabric anisotropy variable within the dilatancy-stress relationship of the ACST.

Meanwhile, direct data evidence regarding these theoretical propositions for the relationship between dilatancy ratio and fabric anisotropy are scarce and often conflicting. Oda (1972) and Arthur and Menzies (1972) conducted tests under a conventional triaxial compression stress path on sand samples, with various angles between the major principal stress axis and the normal to the bedding plane on which the sand is deposited (i.e., bedding plane angle $\phi$). Based on these tests, Oda (1972) and Tatsuoka (1976) suggested that the dilatancy-stress relationship is independent of fabric anisotropy. However, the scarce data points at low stress ratio and the strong scatter of the data undermines the comparison of dilatancy ratio between samples with different bedding plane angles, especially considering that the influence of inherent fabric anisotropy is most prominent at low stress ratio. In contrast, Yoshimine et al. (1998) and Nakata et al. (1998) indicated through undrained hollow cylinder torsional tests that changes in the bedding plane angle of the sample can dramatically affect the dilatancy behavior of sand, ranging from dilative hardening to contractive static liquefaction. Being undrained tests, the dilatancy ratio could not be directly analyzed based on these test results.

Under both conventional triaxial compression and undrained torsional stress paths, the evolution of mean effective stress can be drastically different for samples with different bedding plane angles, obscuring the influence of anisotropy. Drained constant mean effective stress hollow cylinder torsional test results from Kandasami and Murthy (2015) show clear variation of dilatancy ratio for samples with different bedding plane angles, especially at low stress ratio. Unfortunately, the resolution of the test data does not allow for scrutiny of the relationship between dilatancy ratio and the bedding plane angle.

For the analysis of the relationship between dilatancy ratio and fabric anisotropy, discrete element method (DEM) (Cundall and Strack 1979) is a powerful complementary tool to physical testing by providing high resolution data, better control in sample repeatability and boundary conditions, and easier access to fabric measurements (Kruyt and Rothenburg 2006; O’Sullivan 2011); and thus, is a means to link stress and micromechanical fabric (Rothenburg and Bathurst 1989; Li and Yu 2013). The few existing DEM studies on this subject yield contradicting results and are still inconclusive (Yunus et al. 2010; Vincens and Nouguier-Lehon 2012; Hosseininia 2012; Jiang et al. 2018). Through applying load reversal at different strain levels during triaxial tests, Yunus et al. (2010) observed that the stress ratio at which the material behavior changes from contractive to dilative is dependent on fabric anisotropy. Vincens and Nouguier-Lehon (2012) discussed the dependency of dilatancy on fabric anisotropy based on observations from cyclic loading. The data from Hosseininia (2012) also shows variation of dilatancy with respect to fabric anisotropy. In contrast, Jiang et al. (2018) suggested otherwise. However, close scrutiny of the results from simulations on a loose sample in Jiang et al. (2018) does, in fact, show some dependency of dilatancy ratio on initial bedding plane angle, which the authors suggested should be further analyzed.

This study aims at investigating the dependency of dilatancy ratio on fabric anisotropy in granular materials through carefully conducted 2D and 3D DEM numerical tests, and to provide important evidence for the development of constitutive theories associating soil behavior with fabric anisotropy. Existing knowledge of dilatancy ratio dependency on void ratio allows for the numerical test program to focus mainly on fabric anisotropy. Detailed descriptions of the numerical test setup are provided. The test results are presented and analyzed, with special focus on the dependency of dilatancy ratio on fabric anisotropy. The quantitative relationship between dilatancy and fabric anisotropy is examined using the test results with ACST, connecting discrete micromechanical observations with continuum constitutive theory.

In this study, 2D DEM biaxial tests are conducted using the Polyarc Parallel-processing Discrete Element Modeling code (PPDEM) (Fu et al. 2012), which has been successfully applied to a number of studies on the mechanical behavior and fabric anisotropy of granular materials (Fu and Dafalias 2011; Wang et al. 2016, 2017a, b, 2019). 3D DEM triaxial tests are conducted using Yade (Šmilauer et al. 2015), a widely used open source code (e.g., Kozicki et al. 2013; Zhao and Guo 2015). 2D and 3D granular materials share the same basic microscopic mechanisms related to particle interactions that govern the macroscopic behavior of the material. Due to computational efficiency considerations, 2D DEM is often adopted to provide adequate qualitative assessment of granular material behavior. However, for investigations with relatively limited existing data, such as that in this paper, it is important to provide corroborating data from both 2D and 3D DEM tests.

Sixteen tests on various samples are conducted in this study to investigate the influence of fabric anisotropy on dilatancy, including eight 2D constant mean effective stress biaxial tests and eight 3D constant mean effective stress triaxial tests. In both 2D and 3D, two sets of samples with different initial void ratios are used. Within each set, four samples with the same void ratio, but different bedding plane angles of 0°, 30°, 60°, and 90°, are adopted to generate variations in initial fabric orientation. The numerical test program is listed in Table 1. Each of the 2D samples consists of approximately 20,000 randomly generated elliptical particles with an aspect ratio of 1.5:1, grain size of $D_{\min }=0.3\mathrm{mm}$, $D_{\max }=1.0\mathrm{mm}$, $D_{50}=0.72\mathrm{mm}$, and interparticle friction angle of 35°. The elliptical particles are generated using a polyarc formulation presented by Fu et al. (2012). In 2D, the interparticle contact normal and tangential stiffness parameters, $K_n$ and $K_s$, are $500\mathrm{GPa}/\mathrm{m}$ and $167\mathrm{GPa}/\mathrm{m}$, respectively, using the contact law described by Fu et al. (2012). For the 3D tests, each sample consists of approximately 40,000 clumped particles with the same aspect ratio of 1.5:1, grain size of $D_{\min }=0.115\mathrm{mm}$, $D_{\max }=0.305\mathrm{mm}$, $D_{50}=0.235\mathrm{mm}$, and interparticle friction angle of 35°. A 3D clumped particle is created by rigidly connecting two identical overlapping spherical particles. The efficacy of clumped particles in representing the behavior of sand has been justified in numerous studies (e.g., Kuhn 2017; Li et al. 2016). A traditional linear elastic–plastic contact law (Cundall and Strack 1979) is used for interparticle contact in the 3D tests, with $K_n=5\times {10}^{5}r\mathrm{kN}/{\mathrm{m}}^2$ and $K_s=0.3K_n$, where $r$ equals the radius of the particle.

In both 2D and 3D, the particles are first pluviated under gravity. The pluviated assembly is rotated by the desired bedding plane angle $\phi$, and then each sample (consisting of approximately 20,000 and 40,000 particles in 2D and 3D, respectively) is obtained from within the assembly and consolidated under 100 kPa isotropic stress. The stress after pluviation prior to consolidation is within 5 kPa, far smaller than the consolidation stress. Following this procedure, which has been presented in detail in Fu and Dafalias (2011), two sets of samples with initial void ratios $e_{in}$ of $0.2150\pm 0.0004$ and $0.1820\pm 0.0004$ are generated in 2D, and two sets of samples with initial $e_{in}$ of $0.690\pm 0.005$ and $0.656\pm 0.005$ are generated in 3D, as listed in Table 1. The unit of stress in 2D is presented in meters for consistency with 3D. The homogeneity of the samples are checked rigorously following the procedures proposed by Fu and Dafalias (2011). The damping parameters and load rates are chosen to achieve quasi-static loading with the simulation results being insensitive to moderate variation of their values. Two pairs of loading walls in the vertical and horizontal directions are used to apply the biaxial load in 2D, while three pairs of loading walls in the vertical and two in the horizontal directions are used to apply the triaxial load in 3D. Strain controlled loading is applied, during which constant $p$ is achieved through a servomechanism. The loading walls are frictionless to allow uniform application of stress. Diagrams for the samples, the loading schemes, and typical initial fabric in 2D and 3D are illustrated in Fig. 1.

The stress within the samples are calculated following Bagi’s (1996) formulation. In 2D and 3D, the deviatoric stress is the difference of the major and minor principal stresses, expressed as $q={\sigma }_1–{\sigma }_2$ and $q={\sigma }_1–{\sigma }_3$, respectively, while the mean effective stress is expressed as $p=({\sigma }_1+{\sigma }_2)/2$ and $p=({\sigma }_1+{\sigma }_2+{\sigma }_3)/3$, respectively. For the triaxial tests in this study, ${\sigma }_2={\sigma }_3$. The deviatoric stress ratio $\eta =q/p$. The strain of the samples are calculated following the procedure by Fu and Dafalias (2012). The volumetric strain is given by ${\epsilon }_v={\epsilon }_1+{\epsilon }_2$ and ${\epsilon }_v={\epsilon }_1+{\epsilon }_2+{\epsilon }_3$ for 2D and 3D, respectively, while the deviatoric strain is given by ${\epsilon }_q=1/2({\epsilon }_1-{\epsilon }_2)$ and ${\epsilon }_q=2/3({\epsilon }_1-{\epsilon }_3)$ for 2D and 3D, respectively. In this paper, compression is considered positive for both stress and strain, following the sign conventions in soil mechanics. It should be noted for samples with oblique bedding planes in 3D, ${\epsilon }_2$ and ${\epsilon }_3$ may be different due to fabric anisotropy, and ${\epsilon }_q=2/3({\epsilon }_1-{\epsilon }_3)$ is an approximation considering that only a small difference in ${\epsilon }_2$ and ${\epsilon }_3$ is observed.

The plastic strain increment $\mathrm{d}{\mathbf{\epsilon }}^p$ generated during loading can be calculated by subtracting the elastic strain increment $\mathrm{d}{\mathbf{\epsilon }}^e$ from the total strain increment $\mathrm{d}\mathbf{\epsilon }$; $\mathrm{d}{\mathbf{\epsilon }}^p=\mathrm{d}\mathbf{\epsilon }-\mathrm{d}{\mathbf{\epsilon }}^e$. The elastic strain increment $\mathrm{d}{\mathbf{\epsilon }}^e$ is obtained using the procedure suggested by Calvetti et al. (2003) and Wan and Pinheiro (2014). Numerical snapshots of the DEM sample are taken at designated intervals during loading, and a probe in the loading direction with no contact sliding is conducted using each snapshot to obtain the elastic strain increment. This process is computationally costly, but is important when trying to bridge the gap between DEM and continuum models. For the data reported in this study, the deviatoric strain increment is set to 0.0005 during the strain hardening loading stage before the occurrence of peak stress, in which drastic changes in stress and dilatancy occur. After peak stress, during the strain softening stage, the deviatoric strain develops significantly while stress and dilatancy changes slowly, and the deviatoric strain interval is set to a larger value of 0.02.

The second order fabric tensor $\mathbf{N}$ reflecting the anisotropy of granular materials is formulated as (Satake 1982)where $N$ = count of the entity (e.g., contacts, particles) being quantified in the domain; $\mathbf{v}$ = unit norm vector of the directional entity that the fabric tensor is based on; $\otimes$ = tensor product; and superscript $k=k$th entity in the domain. Use of the contact normal vector ${\mathbf{v}}_{\mathrm{c}}$ in Eq. (1) results in the contact normal fabric tensor ${\mathbf{N}}_{\mathrm{c}}$, while the particle orientation vector ${\mathbf{v}}_{\mathrm{p}}$, which is the long axis of the particle, yields the particle orientation fabric tensor ${\mathbf{N}}_{\mathrm{p}}$. The fabric anisotropy intensity $\alpha$ can be depicted in 2D by the difference between the major and minor principal fabric tensor components, $N_1-N_2$, and in 3D by $\sqrt{[{(N_1-N_2)}^2+{(N_2-N_3)}^2+{(N_1-N_3)}^2]/2}$, while the orientation of the fabric tensor $\theta$ is referred to as the angle between the major principal fabric tensor axis and the horizontal plane. The initial fabric anisotropy intensity and orientation of all the samples are listed in Table 1. In general, the looser samples have greater initial fabric anisotropy intensities compared to denser samples for both 2D and 3D.

In ACST proposed by Li and Dafalias (2012), to establish the connection between fabric anisotropy of granular materials with continuum constitutive theory, the effects of fabric anisotropy on the dilatancy ratio is incorporated via

In Eq. (2), $M$ = critical state stress ratio, and $d$ and $m$ are often assumed to be material constants in constitutive modeling (e.g., Dafalias and Manzari 2004; Wang et al. 2014). $\zeta$ = the dilatancy state parameter expressed aswhere $e_c$ = critical state void ratio that is a function of the mean effective stress $p$; and ${\widehat{e}}_A$ = a function of the void ratio $e$ and/or the mean effective stress $p$, in general. For the 2D and 3D materials in this study, the critical state void ratios $e_c$ at $p=100\mathrm{kPa}$ are 0.240 and 0.761, respectively. $e_c$ is determined as the mean of all eight samples of the same material after the void ratio reaches a steady state, at 0.8 shear strain. The fabric anisotropy variable $A=\mathbf{F}:\mathbf{n}$. $\mathbf{F}$ is a second order deviatoric fabric tensor normalized at the critical state, that is:where $r=2$ and 3 for 2D and 3D cases, respectively; the fabric tensor calculated by means of DEM, is here normalized by the specific volume $1+e$ to become a per volume measure, achieving thermodynamic consistency with the continuum definition of fabric (Li and Dafalias 2015); and $\mathbf{n}$ = the unit-norm deviatoric tensor-valued loading direction. For the 2D biaxial loading in this study, the unit-norm deviatoric tensor-valued loading direction is constant at

For the 3D triaxial loading, the unit-norm deviatoric tensor-valued loading direction is simplified by neglecting the difference in the two lateral axial strain components caused by fabric anisotropy

According to the ACST, $A=\mathbf{F}:\mathbf{n}=1$ should be true at the critical state, indicating that $\mathbf{F}$ and $\mathbf{n}$ have the same orientation.

The evolution of the fabric anisotropy variable $A$ under 2D biaxial load and 3D triaxial load are presented in Figs. 10 and 11, respectively, based on both contact normal and particle orientation fabric measured from the numerical tests. Note that when the particle orientation fabric is used, in order for the particle orientation-based fabric anisotropy variable $A_{\mathrm{p}}$ to be positive at the critical state, $A_p=-{\mathbf{F}}_p:\mathbf{n}$. For samples with different initial bedding plane angle $\phi$ values, the initial fabric anisotropy variable $A$ values are distinctly different. During loading, contact normal-based fabric anisotropy variable $A_{\mathrm{c}}$ evolves gradually, while $A_{\mathrm{p}}$ remains almost constant until high $\eta$ values. At the critical state, both $A_{\mathrm{c}}$ and $A_{\mathrm{p}}$ reach unity, as is suggested in the ACST, regardless of the sample’s initial state. As shown in Figs. 10 and 11, the difference of the initial $A$ values between samples with different bedding plane angles is greater within the looser sets, especially for the 2D samples, because the looser samples have greater initial fabric anisotropy intensity, echoing the greater difference of the initial dilatancy ratio within the looser sets of samples shown in Figs. 4 and 5.

The 3D tests of $e_{in}=0.690\pm 0.005$ with $\phi =0°$ and $\phi =90°$ are used as examples to analyze the influence of $A$ on dilatancy. According to the ACST, at the initial state, the only difference for the two samples in terms of Eqs. (2) and (3) is the fabric anisotropy variable $A$. For the test with $\phi =0°$, at the initiation of loading, $A_{c,\phi =0°}={\mathbf{F}}_c:\mathbf{n}=0.37$ when ${\mathbf{F}}_c$ is based on the contact normal fabric. For the test with $\phi =90°$, the initial contact normal fabric yields an initial $A_{c,\phi =90°}=-0.18$. Thus, at the initial state, ${\zeta }_{\phi =0°}<{\zeta }_{\phi =90°}$, resulting in $D_{\phi =0°}<D_{\phi =90°}$, matching the relationship observed in the triaxial tests (Fig. 5). Based on the fabric measurements, Eq. (2) would also predict a greater phase transformation deviatoric stress ratio for $\phi =90°$ than for $\phi =0°$, agreeing with the test results in Fig. 5. Once the contact normal fabric tensor major principal axis transitions from horizontal to vertical, the difference of the fabric anisotropy variable value between the two tests becomes smaller. At $\eta =1.1$, $A_{c,\phi =0°}=0.87$ and $A_{c,\phi =90°}=0.68$. According to Eq. (2), the difference in dilatancy ratio value of the two samples would become much smaller compared to the initial state. This agrees with the observation that at high deviatoric stress ratio, the dependency of dilatancy ratio on initial bedding plane angle diminishes (Fig. 5). However, if the particle orientation fabric is used for $\mathbf{F}$, $A_p$ would still be distinctly different even at $\eta =1.1$, with the difference close to that at the initial state. This further suggests that the contact normal fabric serves as a better indication for the influence of fabric anisotropy in the dilatancy-stress relationship in Eq. (2).

Based on the $\eta$, $M$, $e$, $e_c$, and contact normal based $A_c$ data obtained from the DEM tests, the dilatancy ratio can be theoretically calculated with the ACST framework based on Eqs. (2) and (3), once $d$, $m$, and ${\widehat{e}}_A$ are determined. With the simplified assumption that $d$, $m$, and ${\widehat{e}}_A$ are constants independent of other variables, taking a single set of values of 0.4, 35, and 0.03 for the eight samples in 2D and a single set of values of 0.5, 8, and 0.12 for the eight samples in 3D, the calculated dilatancy ratio for the 2D and 3D tests are plotted against $\eta$ in Figs. 12 and 13, respectively.

Comparisons between the theoretical results in Figs. 12 and 13 with the DEM results in Figs. 4 and 5 show that even under highly simplified assumptions, the ACST is able to capture the overall dependency of dilatancy ratio on fabric orientation in granular materials exceptionally well. The theoretical dilatancy ratio results not only show that tests with greater $\phi$ have greater initial $D$ and phase transformation $\eta$, but they also exhibit the convergence of $D$ as fabric anisotropy evolve toward the same critical state at high $\eta$ in the tests with different $\phi$, as observed in the DEM tests. Comparisons of the calculated dilatancy ratio between samples with different initial void ratio and fabric anisotropy intensity in Figs. 12(a and b) and 13(a and b) further highlight the ACST framework’s ability to depict the dependency of dilatancy ratio on not only fabric orientation, but also fabric anisotropy intensity and void ratio with a single set of parameters. Note that the void ratio, deviatoric stress ratio, and fabric used in the calculation of the theoretical dilatancy ratio are directly obtained from DEM, and are intended to validate the core concept of ACST’s incorporation of fabric anisotropy in the dilatancy ratio.

If the influence of fabric anisotropy is not considered in calculating the dilatancy ratio, the samples with the same void ratio and different bedding plane angles should have the same dilatancy ratio; and looser samples should display greater contraction compared to denser samples. Due to the considerations for the influence of fabric anisotropy within the ACST framework, the denser 2D sample with $\phi =90°$ and $e_{in}=0.1821$ actually results in a greater initial dilatancy ratio in contraction of 0.14, compared to the 0.13 for the looser sample with $\phi =0°$ and $e_{in}=0.2150$, agreeing with the observed DEM results (Fig. 4). The smaller initial fabric anisotropy intensity in the denser sets of samples (around 0.05 and 0.06 for 2D and 3D, respectively) compared with the looser sets (around 0.11 and 0.07 for 2D and 3D, respectively) also leads to the smaller difference in fabric anisotropy variable for samples with different bedding plane angles, which in turn results in the smaller difference of the dilatancy ratio within the denser sets of samples in Figs. 12(b) and 13(b) compared to that within the looser sets in Figs. 12(a) and 13(a).

In regards to the simplified assumptions made in the previous calculations, Li and Dafalias (2012) pointed out that, in general, ${\widehat{e}}_A$ should be a function of $e$ and/or $p$, rather than being constant. Although the state dependency of Eq. (2) naturally introduces nonlinearity in $D\text{-}\eta$ relationship, more significant nonlinearity of the $D\text{-}\eta$ curves is observed for DEM results in Figs. 4 and 5 compared with that for theoretical results in Figs. 12 and 13. This suggests that the parameter $d$ may also be a function of the fabric anisotropy variable $A$, rather than being a material constant, as is assumed in the theoretical calculations. Nonetheless, these results can overall serve as direct proof that the formulation for dilatancy ratio in the ACST can provide an excellent mathematical framework depicting the relationship between dilatancy ratio and fabric anisotropy with an appropriate choice of fabric definition.

DEM numerical tests in both 2D and 3D are conducted in this study to directly reveal the relationship between dilatancy ratio and fabric anisotropy for granular materials. The 2D and 3D test results are shown to be qualitatively alike, providing mutual corroboration. The dilatancy ratio is shown to be significantly influenced by the initial fabric orientation and anisotropy intensity. The dependency of dilatancy ratio on initial fabric orientation is analogous to the dependency of dilatancy ratio on void ratio. The difference between samples with greater and smaller bedding plane angles, that is, the angle between the major principal stress and the normal to the bedding plane, is similar to that between looser and denser samples. Samples with greater bedding plane angles exhibit greater dilatancy ratio (positive in contraction and negative in dilation) at the initiation of loading, greater deviatoric stress ratio at phase transformation, and smaller absolute value of peak negative dilatancy ratio, compared with samples with smaller bedding plane angles. For samples with the same void ratio and fabric anisotropy intensity but different bedding plane angles, greater initial fabric anisotropy intensity results in more significant differences in initial dilatancy ratio.

The dilatancy ratio is observed to evolve along with fabric evolution during loading, instead of only being dependent on the initial fabric anisotropy. For samples with distinctly different initial fabric orientations, as the deviatoric stress ratio, void ratio, and fabric tensor evolve towards the same critical state, the difference in dilatancy ratio between the samples also diminishes, the dilatancy ratio becomes independent of the initial fabric anisotropy, and the data between peak dilatancy ratio and the critical state converges towards the same path. In this study, the dilatancy ratio is shown to be dominantly influenced by the contact normal fabric, while the particle orientation fabric plays a less significant role. The relationship between various fabric tensors, defined based on different entities such as particle, contact, and void, with the dilatancy ratio of granular materials should be further investigated. The evolution of these fabric tensors should also be studied in detail in future studies. Such studies will further bridge the gap between micromechanical fabric observations and constitutive modeling efforts in the macro scale.

By addressing the value of dilatancy using both DEM and ACST, direct proof for the cornerstone of the theory is achieved in this study. Based on the numerical test data, ACST is shown to be able to capture the phenomenon related to the dependency of dilatancy ratio on fabric anisotropy with an appropriate choice of fabric tensor definition. Using a single set of parameters, the ACST framework is able to correctly reflect the influence of fabric orientation, fabric anisotropy intensity, and void ratio on dilatancy ratio for the same granular material. The ACST framework bridges the gap between micromechanical fabric anisotropy and macromechanical behavior of granular materials, and provides a powerful framework for incorporating the dilatancy anisotropy in constitutive modeling efforts. The observation from the tests suggests that the slope of the $D\text{-}\eta$ relationship may also be affected by fabric anisotropy, which is often considered constant in existing constitutive models. With further accumulation of laboratory and numerical test data, the currently assumed linear relationship between the dilatancy state parameter $\zeta$ and the fabric anisotropy variable $A$ can also be discussed. A fully autonomous anisotropic constitutive framework with formulations for general stress–strain relationship and fabric evolution should be explored in future studies.

The authors would like to acknowledge the National Natural Science Foundation of China (Nos. 51708332 and 51678346) for funding the work. The authors would also like to thank the reviewers for their constructive comments and suggestions.