An evolution law for fabric anisotropy and its application in micromechanical modelling of granular materials

Micromechanical studies of granular materials have demonstrated the importance of their microstruc- ture to their behaviour. This microstructure is often characterized by fabric tensors. Experimental and computational studies have shown that the fabric can change signiﬁcantly during deformation. Therefore, the evolution of fabric is important to constitutive modelling. Current fabric evolution laws for granular materials have generally been developed for continuum-mechanical models, and use a loading index mul- tiplier associated with a yield surface. Such evolution laws can not be employed with micromechanical models that do not involve an explicit macro-scale yield surface. This study develops an evolution law for fabric anisotropy, based on observations from experiments and Discrete Element Method simulations from literature. The proposed evolution law considers the effects of inherent anisotropy, void ratio, stress ratio, loading direction and intermediate principal stress ratio. In the critical state, the value of the fabric anisotropy depends only on the Lode angle. The predicted evolution of fabric anisotropy is in good agreement with results of Discrete Element Method simulations, showing both hardening and softening behaviour and describing the inﬂuence of the initial void ratio. The proposed evolution law can be embedded into micromechanics-based constitutive relations as well as conventional continuum-mechanical models. As an example, a well-established micromechanical model (in which the fabric is considered as constant) has been extended by accounting for the variations in fabric, in combination with the proposed fabric evolution law. The performance of this enhanced micromechanical model has been demonstrated by a comparison between the predicted behaviour and ex- perimental results from literature for Toyoura sand under various loading conditions.


Introduction
The mechanical behaviour of granular materials is important in many disciplines of engineering and science. For quasi-static deformation, this behaviour is described by constitutive relations (that relate increments of stress and strain) that generally are formulated in a heuristic manner within the framework of classical continuum mechanics.
Since granular materials consist of a large number of particles and voids, the alternative (and complementary) micromechanical viewpoint has been developed, where (constitutive) relationships are investigated between the behaviour at the micro-scale of particles and contacts and the continuum macro-scale. Wang et al., 2019, Yang et al., 20 07, 20 08, 2013Yoshimine et al., 1998 ). Therefore, considerations on the evolution of fabric anisotropy are important to constitutive modelling of granular materials ( Li and Dafalias, 2011;Gao et al., 2014;Gao and Zhao, 2017;Yang et al., 2018 ).
Various fabric tensors have been defined in the literature ( Oda, 1972;Cambou et al., 2013 ), that differ in the involved micro-scale quantities. The considered micro-scale quantities are: branch vectors that connect the centres of particles in contact, particle orientation vectors, contact normal vectors and so-called void normal vectors ( Yang et al., 2008;Li, 2016;Wang et al., 2017 ). In order to appropriately consider the evolution of microstructure of granular materials, as characterised by a fabric tensor, under external loading conditions, three aspects should be accounted for: (1) an initial value, (2) an ultimate value in the critical state and (3) an evolution law ( Yang et al., 2018 ).
The initial, or inherent, fabric can be isotropic or anisotropic, depending on the method of sample preparation. Dry deposition and moist damping methods have been widely used in laboratory experiments ( Yang et al.,20 07;20 08 ) for obtaining samples with different inherent anisotropy. In DEM simulations ( Fu and Dafalias, 2011;Zhao and Guo, 2013;Kruyt and Rothenburg, 2014;2016;Yang and Wu, 2016;Huang et al., 2014 ), samples with inherent anisotropy have been generated by gravity deposition or by the radius expansion method. Both experimental investigations and DEM simulations show that inherent anisotropy significantly affects the shear behaviour of granular materials ( Yoshimine et al., 1998;Yang et al., 2007;Kruyt, 2012;Fonseca et al., 2013;Fu and Dafalias, 2015;Li, 2016;Yang et al., 2013;Kruyt, 2003 ).
In Critical State Theory ( CST for short) ( Roscoe et al., 1958;Schofield and Wroth, 1968 ) the critical state of granular materials is the state in which the material keeps deforming in shear with constant volume and under constant stress. The uniqueness of fabric anisotropy in the critical state has long been controversial. Recently, observations from results of DEM simulations have shown that the fabric anisotropy approaches asymptotic values ( Fu and Dafalias, 2011;Zhao and Guo, 2013;Kruyt and Rothenburg, 2014;2016;Yang and Wu, 2016 ). These critical states have been explicitly considered in the Anisotropic Critical State Theory ( ACST for short) ( Li and Dafalias, 2011 ) by introducing a term involving fabric into the classical CST ( Roscoe et al., 1958 ). Using the framework of ACST, elasto-plastic constitutive models ( Gao et al., 2014;Gao and Zhao, 2017;Petalas et al., 2019Petalas et al., , 2020Woo and Salgado, 2015;Yang et al., 2018 ) have been developed with which non-coaxial deformation of granular materials can be described naturally.
When fabric is incorporated into constitutive models, a corresponding evolution law becomes necessary. Such evolution laws for fabric anisotropy have generally been developed for continuummechanical models ( Li and Dafalias, 2011;Gao et al., 2014;Gao and Zhao, 2017; Nemat-Nasser, 20 0 0; Sun and Sundaresan, 2011;Oboudi et al., 2016;Oda, 1993;Wan and Guo, 2004;Lashkari and Latifi, 2008;Woo and Salgado, 2015;Yuan et al., 2019;Fang et al., 2019 ). One type of these evolution laws involves the loading index associated with the macro-scale yield surface (see for instance ( Li and Dafalias, 2011;Yang et al., 2018;Gao et al., 2014;Gao and Zhao, 2017;Nemat-Nasser, 20 0 0;Woo and Salgado, 2015;Petalas et al., 2019b )). Although these evolution laws have demonstrated their effectiveness in continuum-mechanical modelling, they can not be employed with micromechanical constitutive models that do not possess an explicit macro-scale yield surface. Other types of evolution laws ( Wan and Guo, 2004;Tian and Yao, 2018;Sun and Sundaresan, 2011;Lashkari and Latifi, 2008;Yuan et al., 2019 ) do not involve the loading index, but they can only describe fabric evolution for limited loading conditions. For instance, the fabric evolution law incorporating only the stress ratio ( Wan and Guo, 2004 ) cannot describe the evolution of the microstructure under loading paths with constant stress ratio. The evolution law that only involves the Lode angle ( Lashkari and Latifi, 2008 ) fails to describe the non-coaxial deformation behaviour of granular materials under proportional loading. Evolution laws from literature are described in more detail in Section 2.3 . Therefore, it is desirable to construct a more general evolution law for fabric anisotropy, with as few parameters as possible.
Micromechanics-based constitutive models for granular materials have recently been developed with homogenization methods ( Chang and Hicher, 2005;Kruyt and Rothenburg, 2004;Misra and Singh, 2014;Nicot and Darve, 2011;Xiong et al., 2017;Yin and Chang, 2009;Zhao et al., , 2019. In this micromechanical approach, the force-displacement relationship is specified at the micro-scale level of interparticle contacts and the stress-strain behaviour at the continuum, macro-scale level is then obtained by averaging over all contact orientations. Such models have demonstrated their effectiveness in describing the behaviour of granular materials under various loading conditions with only few parameters that have physical meaning, and they have been successfully extended to describe the influence of inherent fabric anisotropy ( Hicher and Chang, 2006;Yin and Chang, 2009;Chang and Yin, 2009;Yin et al., 2010;Chang and Bennett, 2015 ). However, the evolution of the fabric tensor in these models still has not been accounted for.
Since micromechanical models do not have explicit macro-scale yield surfaces, current fabric evolution laws that generally have been formulated for elasto-plastic continuum-mechanical models can not be employed in micromechanical models. The objectives of the current study are therefore to: • analyse the parameters that influence the evolution of fabric anisotropy, based on experimental and DEM simulation results from literature; • develop a novel evolution law that accounts for these parameters of influence and that can be embedded in micromechanical as well as continuum-mechanical models; • extend the well-established Chang-Hicher micromechanical model ( CH -model for short) ( Chang and Hicher, 2005;Zhao et al., 2018e ) (in which the fabric is considered as constant) by accounting for the variations in fabric, in combination with the proposed fabric evolution law; • demonstrate the effectiveness of the resulting enhanced CHmicromechanical model for predicting the behaviour of granular materials under proportional loading.
The following notations are adopted here. The sign convention from soil mechanics is followed, according to which compression is considered positive. The strain tensor is denoted by , with corresponding plastic strain tensor p and deviatoric strain tensor q . The principal values of the stress tensor σ are denoted by σ 1 , σ 2 , σ 3 , with σ 1 ≥ σ 2 ≥ σ 3 . The intermediate principal The mean pressure is given by where I is the second-order identity tensor. The stress deviator q is defined by q = 3 2 s : s . The stress ratio tensor η and the (scalar) stress ratio η are given by The rate of change of some quantity A (such as the stress, strain and the fabric tensors) is denoted by ˙ A . The Lode angle θ A of a symmetric second-order A is defined in terms of its (ordered) principle values A 1 , A 2 , A 3 by The overview of this study is as follows. In Section 2 the parameters are identified that have a significant effect on the evolution of fabric anisotropy. A corresponding evolution law is proposed in Section 3 . It is also shown here that the proposed fabric evolution law well reproduces results from three-dimensional DEM simulations from literature. In Section 4 an overview of the CH-micromechanical model is given and it is explained how this model has been extended by accounting for variations in the fabric tensor. In Section 5 , the capability of the enhanced model is demonstrated by comparing numerical predictions with experimental data from literature by Yoshimine et al. (1998) for Toyoura sand under various loading conditions. Section 6 gives the main findings of this study and discusses further developments.

Observations on fabric anisotropy and its evolution
Definitions of fabric anisotropy for granular materials are reviewed in Section 2.1 . Based on observations from experimental studies and DEM simulations, the most important parameters affecting fabric anisotropy are identified in Section 2.2 . Evolution laws from literature are discussed in Section 2.3 .

Definition of fabric tensor
To describe the microstructure of granular materials, many fabric tensors have been defined ( Oda, 1972;Subhash et al., 1991;Cambou et al., 2013 ). The description of the microstructure is based on unit vectors u c associated with the interparticle contacts c . Examples of such vectors u c are the contact normal vector ( Oda, 1982 ) and the normalised void vector  ). The exact nature of the vector u c associated with a contact is not important in the following.
Based on these vectors u c , the symmetric second-order fabric tensor G is defined as the average of (the dyadic product of) the vectors u c that are associated with interparticle contacts where N c is the total number of contacts and w ( u c ) are scalar weight factors associated with contacts c ( Cambou et al., 20 0 0;. According to Li and Dafalias (2015) the fabric tensor should be volumetrically normalized in order to satisfy thermodynamic principles. The dimensionless fabric anisotropy tensor F , that is studied in the following and that has been widely used in the literature ( Zhao and Guo, 2013;Yang et al., 2018 ), is proportional to the deviatoric part of the fabric tensor G F = 5 2 The (scalar) norm F and the direction tensor n F of the tensor F are defined by Note that the direction tensor n F satisfies n F : n F = 1 and tr ( n F ) = 0 . Since F is a symmetric second-order tensor, its principal values F 1 , F 2 and F 3 (not necessary implying that F 1 ≥ F 2 ≥ F 3 ) are real.
To quantify fabric anisotropy of granular materials, experimental techniques (such as photo-elastic tests , X-ray tomography ( Andò et al., 2013 ) and scanning electron microscope tests ( Yang et al., 2008 )) have been successfully employed. However, it is impractical to experimentally quantify the evolution of fabric anisotropy along all loading stages.
Alternatively, using DEM simulations, the value of fabric anisotropy at the initial stage and its evolution can be quantified in a more convenient manner. Series of DEM simulations have been conducted with various initial void ratios, different particle shapes, and different loading schemes ( Chantawarangul, 1993;Li and Yu, 2009;Fu and Dafalias, 2011;Yang and Wu, 2016;Xie et al., 2017 ). Results of these DEM simulations show that the evolution of fabric anisotropy strongly depends on sample density, inherent anisotropy, loading path and stress level.
Since F is a deviatoric tensor, F 1 = −(F 2 + F 3 ) . For transversely isotropic granular materials, such as granular soils deposited under gravity, F 2 = F 3 . Therefore, for such cases the value of F can be quantified by a single scalar quantity. Based on this idea, several definitions for the quantification of fabric anisotropy have been proposed, see for instance ( Oda and Nakayama, 1988;Yang et al., 2008;Gao and Zhao, 2017 ). Among these, the degree of fabric anisotropy suggested by Oda and Nakayama (1988) has been widely adopted ( Yang et al., 2008;Yang and Wu, 2016 ). For transversely isotropic granular materials, the fabric tensor can be expressed (in the coordinate system whose principle directions are parallel and perpendicular to the direction of gravity) as This expression has been employed in the quantification of DEM simulation results in Yang and Wu (2016) , which are analyzed in Section 2.2 . Fig. 1 shows some typical results of three-dimensional DEM simulations from Chantawarangul (1993) , Yang and Wu (2016) , in which the second-order fabric tensor has been based on the contact normal vectors. Chantawarangul (1993) investigated fabric evolution in samples of spherical particles under shear for different loading directions, i.e. for different intermediate principal stress ratios b . Using three-dimensional samples consisting of so-called "clumped" particles, Yang and Wu (2016) investigated the influence of void ratio, inherent fabric anisotropy and loading direction on the evolution of fabric anisotropy.

Parameters affecting the evolution of fabric anisotropy
In Fig. 1 a the influence of the initial void ratio e 0 of (initially) isotropic samples on the evolution of fabric anisotropy is shown. The fabric anisotropy increases from the initial value of zero for an isotropic sample to a peak value, after which softening occurs. The dense sample has a higher peak fabric anisotropy than the loose sample. The peak fabric anisotropy is reached at a smaller strain for the dense sample than for the loose sample. At large deformations, all samples reach the critical state, with identical fabric anisotropy.
In Fig. 1 b the influence of the inherent (i.e. initial) fabric anisotropy (characterised by in Eq. (6) ) on the evolution of fabric anisotropy is shown. The peak fabric anisotropy increases with increasing inherent fabric anisotropy 0 . Again, all samples reach the critical state, with identical fabric anisotropy, at large deformations.
In Fig. 1 c and d the influence of the loading direction, characterised by the intermediate principal stress ratio b , on the evolution of the fabric anisotropy is shown. The peak fabric anisotropy and the critical-state value increase with increasing intermediate principle stress ratio b -value, see also Zhao and Guo (2013) . However, these important observations have not been adequately considered in fabric evolution laws in the literature.
The influence of the confining pressure on the critical state fabric anisotropy has been investigated in Zhao and Guo (2013) , Yang and Wu (2016) , Huang et al. (2014) , using threedimensional DEM simulations. The results in Zhao and Guo (2013) , Huang et al. (2014) show that the critical state fabric anisotropy is weakly dependent on the confining pressure p , but strongly de- pendent on the b -value, i.e. on the loading direction. The results by Yang and Wu (2016) for the fabric tensor that has been normalised by volume (i.e. the fabric tensor F / (1 + e ) is considered) show that the fabric anisotropy is dependent on the b -value, but not on the confining pressure p .
The norm of F at the critical state is denoted by F c and M is the value of the stress ratio η at the critical state. The values of F c and M can be described by where θ F and θ s are the Lode angle (see Eq. (2) ) of the fabric and the stress tensors, respectively; F c 0 and M c are the fabric anisotropy and the stress ratio in the critical state under triaxial compression.
The function g s ( θ s ) describes the effect of loading path and is given by Dafalias et al. (2004) where M e and M c are the values of η in the critical state under triaxial extension and triaxial compression conditions, respectively. In Eq. (7) , θ s = 0 • corresponds to the loading path of triaxial compression, while θ s = 60 • for triaxial extension. Other expressions (for instance, as described in Sheng et al. (20 0 0) ) can also be adopted instead of Eq. (7) to describe the effect of loading path on M . According to Zhao and Guo (2013) , Li and Dafalias (2015) , the first joint invariant of the deviatoric fabric tensor and the stress tensor has the same value for all loading directions. In other words, the product F c · M is independent of the intermediate principal stress ratio b , which shows that This expression for F c is consistent with that in Li and Dafalias (2015) . As shown in Fig. 2 , the values of fabric anisotropy predicted by Eqs. (7) and (8) agree well (with values of the critical state fabric anisotropy in triaxial compression F c 0 and c given in Zhao and Guo (2013) ) with those obtained from DEM simulations in Zhao and Guo (2013) . This important observation has not been adequately considered in current fabric evolution laws.

Overview of fabric evolution laws from literature
For describing the strongly nonlinear stress-strain behaviour of granular materials, constitutive models are generally developed in rate form. For constitutive models that account for the variations of fabric, it is therefore preferable to also develop a fabric evolution law in rate form.
The fabric tensor G in Eq. (3) is determined by the vectors u c associated with contacts. These vectors in turn are determined by the particle positions, X p for particle p . Since changes in particle positions are determined by the (plastic) strain increment, it is expected that the rate of fabric change ˙ F depends on the plastic strain rate ˙ p .
This argument can be formalized using some (strong) assumptions. If all particles in contact move according to the uniform strain assumption (i.e. ˙ X p = ˙ · X p ; see Rothenburg and Kruyt (2001) , Kruyt and Rothenburg (2001) for limitations of this assumption) and the contact topology is unchanged (i.e. no contacts are created and no contacts are disrupted), then it follows that the rate of fabric change ˙ F is proportional to plastic strain deviator rate ˙ p q , where the proportionality factor depends on the current fabric anisotropy F .
With a fabric tensor based on the contact normal vectors, three mechanisms of fabric change have been identified ( Kruyt, 2012 ): contact creation, contact disruption and contact reorientation. Based on analyses of two-dimensional DEM simulations of isobaric tests, it has been found that contact creation is important at all levels of deformation. Contact disruption is dominant in the small strain range, especially for dense samples. Contact reorientation is only important at large strain.
According to Yuan et al. (2019) , contact creation occurs primarily in the initial phase of shearing, which they consider to be mainly controlled by the change in the stress ratio tensor. They consider the contribution of contact reorientation to be also important, which is driven by the plastic strain rate ˙ p . Based on these considerations, their fabric evolution law involves both rate of the stress ratio tensor, ˙ η, and the plastic strain rate ˙ p .
Focussing on two-dimensional systems, Wan (2016, 2017) have developed separate models for the rates of contact creation and contact disruption, as function of the strain rate and stress rate, respectively. These models are based on a statistical description of the Dirichlet tessellation. These analytical models have been validated by results of two-dimensional DEM simulations of granular assemblies with various initial void ratios.
Based on experimental results with photo-elastic materials Mehrabadi et al., 1988 ), the evolution law for the fabric anisotropy proposed by Nemat-Nasser (20 0 0) incorporates coaxiality between the rates of fabric change and plastic strain. Wan and Guo (2004) argue that the rates of fabric and the stress ratio tensor η are coaxial (referring to the same experiments by Oda et al. (1982) ). Sun and Sundaresan (2011) formulated a heuristic evolution law (involving strain rate ˙ , fabric anisotropy F and void ratio e ), with different terms that counteract and allow for the existence of a critical state where the fabric anisotropy is constant.
The results from the literature analysed in Section 2.2 , combined with the previous descriptions of fabric change mechanisms, show that the rate of change of fabric anisotropy, ˙ F , depends on the deviatoric plastic strain rate ˙ p q , the rate of the stress ratio tensor ˙ η, loading direction ( b -value, determined by the current stress σ), (inherent) fabric anisotropy F and the void ratio e . Therefore, a general evolution law should have the following functional form In the critical state with ˙ p q = 0 , e = e c (p) and ˙ η = 0 , the fabric anisotropy does not change, i.e. ˙ F = 0 . Hence, for consistency with the ACST, H must satisfy The value of the fabric anisotropy F crit in the critical state depends only on its Lode angle, as described in Section 2.2 .
Fabric evolution laws can be categorised with respect to the variables considered to be relevant, i.e. where the fabric rate depends on: (i) stress ratio rate, current fabric tensor and void ratio: ( Wan and Guo, 2004;Lashkari and Latifi, 2008 ); (ii) (plastic) strain rate, current fabric tensor and void ratio: ( Sun and Sundaresan, 2011;Fang et al., 2019 ); (iii) (plastic) strain rate as well as stress ratio rate, current stress and fabric tensors and void ratio: ˙ The most commonly-used fabric evolution laws in the literature are summarised in Table 1 . However, none of these evolution laws completely considers all the described aspects of the evolution of the microstructure of granular materials. Therefore, a novel evolution law for fabric anisotropy is developed in the following Section that is based on the observations from DEM simulations in Section 2.2 and that more comprehensively considers the parameters that influence fabric evolution.

An evolution law for fabric anisotropy
A novel evolution law for fabric anisotropy is formulated in Section 3.1 , based on the observations from experiments and results of DEM simulations that have been summarised in Section 2.2 . In Section 3.2 the developed evolution law is calibrated and compared with results of three-dimensional DEM simulations from literature ( Zhao and Guo, 2013;Yang and Wu, 2016 ).

The proposed rate-form fabric evolution law
The evolution of fabric anisotropy has been quantified in threedimensional DEM simulations by Yang and Wu (2016) , while the value of fabric anisotropy in the critical state has been investigated with DEM simulations by Fu and Dafalias (2011) , Zhao and Guo (2013) , Kruyt and Rothenburg (2014) . From these DEM observations it follows that the fabric rate must depend on (initial) fabric anisotropy, deviatoric plastic strain tensor rate, Lode angle of both stress and fabric and sample density (void ratio).
The fabric evolution law that is proposed here is based on the following assumptions: Table 1 Evolution laws for fabric anisotropy from literature.

Category
Evolution laws for fabric anisotropy References Wan and Guo (2004)
• The fabric rate is proportional to the deviatoric plastic strain rate ˙ p q ; • The dependence of the fabric rate on mean pressure p and void ratio e is only via the state function ψ ( Been and Jefferies, 1985;Li and Wang, 1998 ) where the void ratio in the critical state, e c , depends only on the confining pressure p ; • The stress ratio tensor is important via η and its Lode angle θ s , where the latter is only important for the value of the stress ratio M at the critical state, i.e. M = M(θ s ) ; • The critical state value for the magnitude of the fabric anisotropy, F c , conforms to Eq. (5) , and hence is dependent on • The loading direction n with respect to the fabric tensor is characterised by F : n , as also adopted in Li and Dafalias (2011)  An evolution law that incorporates these assumptions is where h is a material constant that accounts for the influence of the rate of the plastic deviatoric strain tensor; f ( η, ψ) is a function of the stress ratio η and the state parameter ψ that accounts for the influence of void ratio. In the critical state, F c = F : n and the state parameter ψ = 0 , and for the rate of fabric anisotropy to be equal zero, it is necessary that f (η = M, ψ = 0) = 1 where the critical state stress ratio M depends on the Lode angle θ s . Such an evolution law belongs to category (iii) described in Section 2.3 . The essential step for developing an evolution law following the framework in Eq. (12) is to further define the loading direction n and the function f ( η, ψ). As discussed in Li andDafalias (2011, 2015) , the loading direction n , a unit-norm deviatoric tensor, can be related to the direction of the deviatoric plastic flow or to the normal to the yield surface in the deviatoric stress space. Here the loading direction is based on the direction of the deviatoric plastic flow. More specifically, the loading direction n is defined as In the critical state, the value of η equals the critical stress ratio M and the state function ψ becomes zero. Given these considerations, a multiplicative form for the function f ( η, ψ) in Eq. (12) is proposed, leading to the fabric evolution law where h and m are two model parameters.
The term F :  Zhao and Guo (2013) , in the critical state the direction of the fabric and the loading direction are coaxial, thus n F : n = 1 . In addition, the stress ratio η reaches its critical value M and the state parameter ψ becomes zero, hence the fabric anisotropy F reaches the critical state F c , irrespective of the development of the deviatoric plastic strain ˙ p q . Thus, in the critical state F c = F c : n , which closely resembles the relationship A c = 1 defined by Li and Dafalias (2011) (who normalize their fabric anisotropy to the value in the critical state). Therefore, the ultimate state of the fabric anisotropy given by Eq. (13) satisfies the requirement for the fabric anisotropy described by the ACST ( Li and Dafalias, 2011 ).

Comparison of predictions by fabric evolution law and results of DEM simulations
In order to demonstrate the capabilities of the fabric evolution law in Eq. (13) , its predictions are compared to results of threedimensional DEM simulations from literature ( Yang and Wu, 2016 ). These DEM simulations involve cases with different initial void ratios e 0 and different loading paths (different b -values).
The results from Yang and Wu (2016) have been selected, as Yang and Wu (2016) provides information on the evolution of the fabric (based on the contact normal vectors) as well as the information required in the fabric evolution law, Eq. (13) , on the stresses and the void ratio e along all loading stages. The numerical samples in these DEM simulations are composed of threedimensional clumped particles, with an interparticle friction coefficient μ = 0 . 5 . The deviatoric plastic strain rate tensor ˙ p q is approximated here by the deviatoric strain rate tensor ˙ q . This approximation is considered appropriate, since the elastic deformations in granular materials are generally very small (and other data are lacking).
The fabric evolution law, Eq. (13) Table 2 .  In the comparison, the values for the stress ratio η and the state parameter ψ from the DEM simulations of Yang and Wu (2016) have been employed in Eq. (13) . According to Yang and Wu (2016) , the direction of fabric anisotropy changes rapidly (for axial strains up to around 2%) to the loading direction n . Therefore, the value of F : n in Eq. (13) is equal to the norm of F , i.e. F = F : n . Fig. 3 a shows the prediction for the evolution of fabric anisotropy of the medium-to-dense sample ( e 0 = 0 . 576 ) under triaxial compression with an initial confining pressure of 10 0 0 k Pa. The hardening and softening behaviour, as well as the peak values of the fabric anisotropy, are well predicted by the proposed fabric evolution law. In the critical state, the values of the fabric anisotropy for the different initial void ratios reach the specified value F c0 = 0 . 48 . Fig. 3 b compares the predictions of Eq. (13) with DEM simulations under triaxial compression and extension in Yang and Wu (2016) , which shows that the fabric evolution law is capable of describing the effect of loading path on the fabric evolution.
This comparison demonstrates that the proposed evolution law adequately reproduces the results from the DEM simulations, based on a modest number of (additional) model parameters ( h, m and F c 0 ; M c and M e are already required to describe the criticalstate yield surface).
As an example for a micromechanical model, the developed fabric evolution law will be employed to extend the CHmicromechanical model (that is based on a constant fabric and that does not involve an explicit yield surface) in the following Section.

Application of the evolution law in granular micromechanical modelling
In this Section, the developed fabric evolution law is combined with a micromechanical model for granular materials. The wellestablished CH micromechanical model ( Chang and Hicher, 2005 ) has been selected as an example to demonstrate how the evolution law can be combined with micromechanical models, in order to incorporate the important effect of the varying microstructure into such models. The CH micromechanical model ( Chang and Hicher, 2005 ) was initially developed for sand. Further developments by Yin and Chang (2009) , Yin et al. (2010, 2012 , Zhao et al. (2018e) demonstrated its good performance in modelling the mechanical behaviour of sand and clay. Recently, this model has been successfully implemented into finite element code to solve geotechnical boundary value problems ( Zhao et al., 2018a( Zhao et al., , 2018b. However, in these developments the fabric tensor is taken as constant , and thus only isotropic and inherently anisotropic samples have been considered, see Hicher and Chang (2006) , Yin et al. (2010) , Chang and Bennett (2015) . Therefore, an evolution law for fabric anisotropy should be considered in this family of micromechanical models in order to make their predictions more physically realistic.
In the CH-model the interaction between particles is described at the micro-scale of interparticle contacts. This description involves: elastic behaviour, yield criterion and dilatancy relation. Forces and displacements are decomposed into directions normal and tangential to interparticle contacts c , denoted by n and t respectively. The model expressions of the CH micromechanical model are summarized in Table 3 .
At the micro-scale, the Coulomb friction law is adopted to describe a yield function for the contact forces (with components f c n and f c t at contact c ), and a non-associated flow rule is defined for the plastic component of the interparticle deformations (with components δ cp n and δ cp t at contact c ). To account for the inter-locking effect on interparticle friction and dilatancy, a dependence of the contact friction angle φ c p and the dilatancy angle φ c d on the void ratio and the critical state void ratio (conforming to the classical CST) has been adopted. The expression for the stress tensor ( Chang and Hicher, 2005 ) in Table 3 is employed to determine the average macro-scale stress from the interparticle forces. The expression for the strain tensor based on the best-fit hypothesis ( Liao et al., 1997 ) in Table 3 is adopted to determine the average macro-scale strain from the micro-scale interparticle deformations. By using energy conservation considerations at the micro-and macro-scales, an expression for force increments in terms of involving stress increments has been obtained, which is regarded as a so-called "static localisation operator", see also Cambou et al. (1995) , Zhao et al. (2018e,c) . For the numerical implementation of the CH micromechanical model, an implicit integration algorithm has been developed in Zhao et al. (2018c) . Table 3 Formulations of the CH micromechanical model ( Chang and Hicher, 2005;Zhao et al., 2018e ).

Relations
Components Equations Critical state line μ is inter-particle friction angle; k c n 0 is a reference normal stiffness; k tR and k pR are material parameters; f ref is a reference force; r is half of the particle diameter.
The expressions for the average stress and the strain tensors involve averages for sets of contacts. For a large number of interparticle contacts, these averages can be converted to integrals over contact orientations ( Bathurst and Rothenburg, 1988;Rothenburg and Bathurst, 1989 ). For an arbitrary quantity c associated with contacts c , this conversion for the average is given by Here E ( n ) is the contact distribution ( Horne, 1965 ) defined over the unit sphere and ¯ is the group average ( Bathurst and Rothenburg, 1988;Rothenburg and Bathurst, 1989 ) of c over contacts with similar orientations, n c . The latter equality indicates that the integral over the unit sphere is approximated by the Gauss numerical integration method in Bažant and Oh (1986) , with Gauss points (number of Gauss points is NG) that are indicated by α with corresponding weights w ( α).
The contact distribution function E ( n ) present in Eq. (14) can be approximated, with good accuracy, by Kanatani (1984) In the CH-model, all particles are considered as equal-sized spherical particles for which the branch vector l c can be expressed as l c = 2 rn c where n c is the contact normal and r is the average radius of particles. Accordingly, the tensor A ij (present in Table 3 ) is related to the fabric tensor F ij by (see also Kanatani (1984 ), Li and Yu (2011) ) where m V is the contact density, i.e. the number of contacts per volume. If non-spherical particles are to be directly considered in the CH-model, then Eq. (16) should be further verified, for instance by using DEM simulations. Since the expressions for the average stress and strain and for the "localisation operator" involve the actual fabric tensor F , the CH micromechanical model can be extended in a straightforward manner by employing the actual fabric (rather than some constant, initial fabric), where the actual fabric is given by the proposed fabric evolution law, Eq. (13) . This shows the importance of employing the current fabric (rather than an initial fabric) in these expressions in the CH micromechanical model.
As described in Chang and Hicher (2005) , the original CH micromechanical model satisfies the classical CST, with respect to the influence of confining pressure and void ratio. At large deformations, the stress ratio η = η c = M is obtained from the evolution of interparticle friction angle φ c p that is a function of the void ratio e . The void ratio e = e c = e c (p) is controlled directly by the critical state line given in Table 3 . With the introduction of the evolution of the fabric tensor, the effect of loading direction on the evolution of fabric anisotropy is described by the term F : n of the fabric evolution law. In the critical state, this term becomes equal to F c , that is dependent only on the Lode angle. Therefore, the enhanced CH micromechanical model satisfies the ACST developed by Li and Dafalias (2011) .

Capabilities of the enhanced CH micromechanical model
In this section, the new features of the enhanced CH micromechanical model will be demonstrated. In Section 5.1 experimental results from literature ( Yoshimine et al., 1998 ) for Toyoura sand are described, which are subsequently used in Section 5.2 to calibrate the enhanced CH model. The new capabilities of the enhanced CH model will then be demonstrated in Section 5.3 by simulating the behaviour of Toyoura sand under various loading directions. In the Appendix, predictions by the enhanced CH-model (for the stressstrain response and for the fabric evolution) are shown for triaxial compression and compared to results of DEM simulations from Yang and Wu (2016) (see also Section 3 ).

Hollow cylinder tests on granular soils
Hollow cylinder tests (see for instance Yoshimine et al., 1998;Yang et al., 2007 ) have been extensively used in experimental soil mechanics to investigate the effect of fabric anisotropy and loading direction on the behaviour of granular soils. More specifically, the influence has been investigated of the angle α (as indicated in Fig. 4 ) between the directions of the major principal values of the initial fabric tensor and of the stress tensor, at various stress levels (stress deviator q ). This is crucial for understanding the role of fabric on the behaviour of granular soils.
To investigate the mechanical behaviour of anisotropic granular soils, three types of experiments have been generally conducted: proportional loading ( ˙ q > 0 ) and the triaxial extension ( α = 90 • , ˙ q > 0 ) tests are two special cases of proportional loading tests. In order to demonstrate the importance of considering fabric evolution in a micromechanical model for granular materials along all loading stages, proportional loading tests are considered here in order to evaluate the enhanced CH micromechanical model. Rotational shearing and combined loading tests can also be reasonably simulated by the enhanced CH model. However, this will be considered in future studies.
The comprehensive proportional loading tests performed by Yoshimine et al. (1998) on Toyoura sand are considered here. According to Yoshimine et al. (1998) , Verdugo and Ishihara (1996) , Toyoura sand has a mean particle diameter d 50 = 0 . 17 mm, a mini-mum void ratio e min = 0 . 597 , a maximum void ratio e max = 0 . 977 , a uniformity coefficient U c = 1 . 7 and a specific gravity of G s = 2 . 65 . Toyoura sand is a uniform fine sand consisting of subrounded to subangular particles and composed of 75% quartz, 22% feldspar and 3% magnetite.
The techniques for performing the proportional loading tests are described in detail in Yoshimine et al. (1998) . For each value of the intermediate principal stress ratio b = 0, 0.25, 0.5, 0.75 and 1, different values of α = 0 • , 15 • , 30 • and 45 • were experimentally studied. In Section 5.2 the experimental results for b = 0 and α = 0 • have been used to calibrate the enhanced CH micromechanical model, from which the effect of sample density and confining pressure can be well considered. Subsequently, the calibrated model will be used in Section 5.3 to predict the behaviour of granular soils for other values of α, thus demonstrating the new features provided by the incorporation of the fabric evolution law.

Model calibration
In total, 13 parameters are required for the enhanced CH micromechanical model: 3 parameters for the critical state line, 5  parameters for the interparticle contact law and 5 parameters for the fabric evolution law, as listed in Table 4 . The parameters in the critical state line relationship (given in Table 3 ) for Toyoura sand have been calibrated by Li and Wang (1998) , Li and Dafalias (20 0 0) and are widely adopted in literature (see for instance Gao et al., 2014;Gao and Zhao, 2017 ). These parameters, i.e. e ref = 0 . 934 , λ c = 0 . 019 and ξ = 0 . 7 , are also used in the current study.
The particle radius r is taken as half of the mean particle diameter d 50 of the tested Toyoura sand. Based on experimental measurements reported in Yoshimine et al. (1998) Gao et al. (2014) , Gao and Zhao (2017) .
The values of fabric anisotropy in the initial state and in the critical state for the tested Toyoura sand have not been experimentally measured. In the simulations of Gao et al. (2014) , Gao and Zhao (2017) , the norm of critical state fabric anisotropy was taken equal to one, as their fabric tensor has been normalized by its norm. In the current study, this critical fabric anisotropy F c0 = 0 . 48 as obtained from DEM simulations in   (2016) has been adopted. The initial norm of the fabric anisotropy of Toyoura sand, prepared by dry deposition method in Yoshimine et al. (1998) , has been assumed to be 0.45 in Gao et al. (2014) ; Gao and Zhao (2017) . In this study, the initial norm of the fabric anisotropy, F , is taken equal to 0.208.

Yang and Wu
The values of k c n 0 , k tR , k pR , h and m have been determined by simulating the stress-strain behaviour of Toyoura sand under undrained triaxial compression ( Yoshimine et al., 1998 ). The values of the calibrated model parameters are given in Table 4 , which are slightly different from those given in Zhao et al., 2018aZhao et al., , 2018bZhao et al., , 2018c for Hostun sand for the original CH-micromechanical model.
For undrained triaxial compression, experimental and simulated results are presented in Fig. 5 for samples of different densities e 0 under various confining pressures p . For the same confining pressure, the simulation results are softer than the experimental results at small deformations, while they are stiffer at large deformations for lower confining pressures. In all cases, a unique stress ratio has been obtained at large deformations. Data are shown here only when the shear strain is smaller than 15%, corresponding to the range in the experiments.

Model predictions
With the calibrated parameters shown in Table 4 , the enhanced CH-micromechanical model has been used to predict the stressstrain behaviour of Toyoura sand for various values of α, as measured by Yoshimine et al. (1998) . Here typical simulation results are shown for the intermediate principal stress ratio b = 0 , with α varying from 0 • to 45 • .
The comparison between the stress-strain relations in the experiments from Yoshimine et al. (1998) and those predicted by the enhanced CH-micromechanical model is shown in Fig. 6 , for a confining pressure of 100 kPa. The influence of the initial fabric anisotropy, represented by α, on the behaviour of granular soils is well captured by the enhanced CH-micromechanical model. With increasing value of α, the experimentally observed smaller shear strength has been reproduced by the model. It is clear that the samples are far from the critical state, as the deformations are very small in comparison to the large deformation (normally when γ ≥ 50%) required for the critical state. Note that at large defor-mations, the enhanced CH-micromechanical model shows a unique critical state.
In order to highlight the importance of considering an evolution law for fabric anisotropy in the CH-micromechanical model, simulation results for the original CH-micromechanical model (see Chang and Hicher, 2005;Zhao et al., 2018e;Zhao et al., 2018c ), with a constant value of the fabric anisotropy (obtained by setting the parameter h = 0 ; see also Eq. (13) ) have also been obtained. In Figs. 6 and 7 it is shown that, when the evolution of fabric is not accounted for, the stress ratio η at critical state is not unique, and thus it is not consistent with the ACST developed by Li and Dafalias (2011) .
In comparison to the original CH-micromechanical model in Chang and Hicher (2005) , Zhao et al. (2018e,c) , the advantages of the enhanced CH-micromechnical model are: (i) the evolution of the microstructure of the granular material under loading is accounted for; (ii) the unique critical state of the ACST framework is retained; (iii) the enhanced model is more accurate for large deformation problems, such as in modelling pile installation and slope stability where strain localization occurs.
Overall, the enhanced CH-micromechanical model is capable of describing the behaviour of granular soils under proportional loadings. As shown, it is crucial to account for fabric evolution in constitutive modelling of granular materials subjected to large deformations, preferably consistent with the framework of ACST such that the unique critical state of granular materials (observed from DEM simulations) is retained.

Conclusions
The microstructure of granular materials, as quantified by fabric tensors, is important to their behaviour. Therefore, fabric is taken into consideration in some continuum-mechanical and micromechanical models. Results of experimental and DEM simulations show that fabric changes significantly during deformation. Hence, an evolution law for fabric anisotropy of granular materials has been proposed, that is based on observations from results of experiments and three-dimensional DEM simulations from literature. This fabric evolution law can be used in combination with both continuum-mechanical and micromechanical models.
The main contributions of this study are: • The evolution of fabric anisotropy depends on the void ratio, the initial fabric anisotropy and the loading directions. These factors have been considered in the developed fabric evolution law, Eq. (13) . This evolution law satisfies the framework of ACST developed by Li and Dafalias (2011) . • The developed fabric evolution law, that involves only few model parameters, has been applied to simulate the results of three-dimensional DEM simulations from literature. The predicted evolution of fabric anisotropy is in good agreement with results of these DEM simulations, showing both hardening and softening behaviour and describing the influence of the initial void ratio. • Considering that the fabric tensors in many micromechanical models in the literature are taken as constant, the developed fabric evolution law has been used to extend the well-established CH micromechanical model for granular soils by accounting for the variations in fabric. The enhanced CHmicromechanical model also satisfies the framework of ACST developed by Li and Dafalias (2011) . • The capabilities of the extended CH micromechanical model have been demonstrated by comparison with hollow cylinder tests from literature. With this extension, the CH micromechanical model is capable of predicting the dependence of the behaviour of granular materials (with evolving microstructure) from various initial states (with differences in fabric, void ratio and confining pressure) to a unique critical state.
For future studies it is recommended to study fabric evolution and the predictions by the proposed evolution law in rotational shearing ( ˙ α > 0 , ˙ q = 0 ) and combined loading tests ( ˙ α > 0 , ˙ q > 0 ) and also to implement the proposed fabric evolution law into continuum-mechanical models. In particular, DEM simulations should be performed to verify the assumption present in Eq. (13) of proportionality between the rates of fabric change and of plastic strain increment under non-proportional loading conditions. It is also recommended to consider other influences, such as cyclic loading, particle shape and particle breakage, on fabric evolution.
The enhanced CH-micromechanical model can be extended to unsaturated granular materials by introducing the capillary cohesion (for instance Kruyt and Millet, 2017;Zhao et al., 2018a;Hicher and Chang, 2007 ), and can also be further applied to geotechnical engineering boundary value problems, using the algorithms in Zhao et al. (2018c) .

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Polytechnic University), who provided essential guidance during the implementation of the enhanced CH-micromechanical model, as well as to Professor Zhongxuan Yang (Zhejiang University) for stimulating discussions and sharing the DEM simulation results reported in Yang and Wu (2016) .

Appendix: Comparison between predictions of the enhanced CH-model and DEM simulation results
Here the predictions by the enhanced CH-model are compared with DEM simulation results of Yang and Wu (2016) in order to demonstrate that the enhanced CH-model has the capability of simulating both the fabric evolution and the stress-strain curve for these data. The cases considered are for isotropic medium-dense samples in triaxial compression under various confining pressures ( Fig. 4 of Yang and Wu, 2016 ).
In the simulations with the enhanced CH-model, the parameters for the critical state line were obtained from Yang and Wu (2016) , and the parameters for the fabric evolution law were the same as those used in Section 3 ( Table 2 ). d 50 = 0 . 535 mm was obtained from Yang and Wu (2016) . The calibrated values of the remaining four CH-model parameters are: k c n 0 = 18 N/mm, φ c μ = 25 • , k tR = 0 . 6 , and k pR = 0 . 5 . Fig. 8 (a,b) show that the enhanced CH-model is capable of reproducing the stress-strain relation obtained from DEM simulations reported in Yang and Wu (2016) . Fig. 8 (c) gives the predictions of the enhanced CH-model (by involving Eq. (13) ) for the evolution of fabric anisotropy (for which no data have been published in Yang and Wu (2016) ). This figure also demonstrates that the fabric evolution is not sensitive to the confining pressure, which is also shown for a dense sample in Fig. 11 of Yang and Wu (2016) .  ( Yang and Wu, 2016 ) (symbols) and predicted by the enhanced CH model (solid lines); initial medium-dense sample is isotropic, under confining pressures of 100 kPa, 500 kPa and 1000 kPa in triaxial compression respectively: (a) deviatoric stress vs. deviatoric strain; (b) deviatoric stress vs. mean stress; (c) fabric anisotropy vs. deviatoric strain.