A constitutive model for granular materials with evolving contact structure and contact forces—Part I: framework

This and the companion paper present a constitutive model for granular materials with evolving contact structure and contact forces, where the contact structure and contact forces are characterised by some statistics of grain-scale entities such as contact normals and contact forces. And these statistics are actually the “fabric” or “force” terms in the “stress–force–fabric” (SFF) equation. The stress–strain response is obtained by inserting the predicted “fabric” or “force” terms from evolution equations into the SFF equation. Discrete element modelling is used to verify the slightly modified SFF equation and also to obtain the data of how the contact structure and contact forces evolve in various loading paths. It is demonstrated that a normalised contact force is a better measure of the contact forces in polydisperse granular assemblies and strong contacts should be contacts with larger normalised contact forces. The modified SFF equation is shown to predict the stress accurately. The constitutive equations regarding the response of the contact structure and contact forces are presented and they along with the SFF equation form a constitutive model, which is found capable of capturing the observed phenomena correctly and predicting the mechanical response in various loading conditions. The model is shown to be an extension to the hypoplastic models with more state variables.


List of symbols
Deviatoric tensor measuring anisotropy of contact normals Directional coordination number D Diameter of grain E Elastic modulus of grain E PDF (n) Probability density function of contact normal e Void ratio f gc Contact force on grain g at contact point c f n,gc Normal contact force f t, gc Shear contact forcẽ f gc Normalised contact force Deviatoric tensor measuring anisotropy of normal contact forces t Deviatoric tensor measuring mobilisation of contacts I Inertia number k n,t,r Elastic stiffness of contacts m g Mass of grain g M gc Contact moment on grain g at contact point c n gc Contact normal on grain g at contact point c N g Number of grains N c Number of contacts p Mean stress q Deviatoric stress in triaxial settings r g Position vector of grain g u t, gc Relative shear displacement V Volume Z Coordination number

Introduction
The constitutive modelling of granular materials is a popular topic [1][2][3] in the research community because it is important not only in understanding the material but also in the numerical investigation of various geotechnical problems [4][5][6]. Granular materials are conventionally modelled as continuum media because a body of interest in the problem-scale (such as a landslide mass) can be continually subdivided into infinitesimal elements with similar properties to those of the bulk material, which is due to the fact that there are still a great number of grains in the infinitesimal elements such that the fluctuation of macro-measurable entities is negligible. At every continuum point, physical quantities should be actually seen as statistics of the grain-scale entities over a representative volume element (RVE), which contains a large number of grains and voids. The classic continuum mechanical descriptions, such as the yield surface, flow rule and hardening rule, are summarised from observations of experiments [2,7]. With recent developments in experimental technology [8][9][10] and grain-based numerical algorithms [11,12], direct observation and quantitative measurement of grain-scale data and processes offer researchers opportunities to study and inspect granular materials at the grain-scale. Oda [8] was among the first to study the anisotropy of contact structure and grain orientations (both are the fabric). Kanatani [13] studied directional functions such as the probability density function of contact normal, approximated them with Fourier-Laplace series, defined several fabric tensors and unveiled their relationship. Through experiments [8][9][10]14], the soil has been shown to be highly anisotropic in terms of fabric entities associated with orientation of particles, voids, contact normal vectors, etc. and this anisotropy of fabric significantly influences the response of soils. These micromechanical findings have inspired and been incorporated into classic models. For example, Dafalias and Li [15] developed a model for inherently anisotropic sands. In the model, some classic ingredients such as the critical state line and plastic modulus are functions of a scalar-valued parameter measuring the inherent anisotropy. In their later models [16], a fabric tensor enters the framework as an internal variable and a rate equations of evolution is developed for it.
In terms of analytical study in micro-mechanics, Rothenburg and Bathurst [17] were the first to realize that average contact forces are also directionally distributed and they presented a stress-force-fabric (SFF) relationship for two-dimensional (2D) systems, which establishes a connection between the stress state and the grain-scale measures of contact structure ("fabric" term) and contact forecs ("force" term). A large number of related studies have been conducted in the next several decades, which are mostly about exploring how anisotropic features influence the shear resistance of granular materials and how the "force" and "fabric" terms change under various kinds of loadings [18,19] because the SFF equation is only an equation of stress and does not explicitly contain deformation. One possible approach of constitutive modelling is that, if evolution equations for both "force" and "fabric" terms under deformation are developed, the predicted "force" and "fabric" terms are then inserted into the SFF equation and a stress-strain relationship is naturally obtained. Therefore, in this framework, the SFF equation and evolution equations form the full constitutive model, which is also the primary aim of the present study. However, these "force" and "fabric" terms are defined on grain-scale entities, which are very hard to determine unless highly idealised discrete element modelling (DEM) is used. In real sands, not only the contact are complex, but also the shape of grains is irregular. The DEM simulations are far from capturing the physical picture. Thence the proposed study may be more appealing in the understanding of some grain-scale mechanism and also possibly in giving some insights to constitutive modelling, rather than in numerical investigations where model parameters are calibrated from tests of real materials.
In this part, the stress and strain obtained from DEM tests serve as the "experimental" results of our virtual granular material. Grain-scale entities are also recorded in DEM and these data are used to calculate the the "force" and "fabric" terms. These terms are firstly inserted into the SFF equation to verify the accuracy. Secondly, these terms can also serve as observations of how contact structure and contact forces change under deformation and rate equations with model parameters are proposed for them. In the constitutive model, the "force" and "fabric" terms are not from DEM results any more, but from rate equations and these terms are inserted into the SFF equation to predict the stress, which is compared to the stress from DEM to examine the performance of the model. The detailed description of the granular assembly of interest, contact model, parameters, simulation procedure is given in Sect. 2. In Sect. 3, firstly, the SFF equation is briefly recapped and related notations are introduced. and we show that for polydisperse granular assemblies, strong contacts should be the contacts with larger normalised contact force and also the benefit of using normalised contact force in SFF analysis. In Sect. 4, the normalised contact force leads to our slightly different SFF equation and in our derivation, the uncorrelated assumption between contact vector and contact force is not necessary. The constitutive model and some discussions are presented in Sect. 5.

Discrete element modelling
Discrete element simulations have been used to verify the modified stress-force-fabric relationship and also to obtain the data of contact structure and contact forces such that their response under deformation could be observed, summarised and a constitutive model could be built.
The commercial software PFC 3D is used in this study. The granular material of interest is an assembly of spherical grains, which have a mean diameter D of 0.17 mm with ± 20% dispersion (uniform distribution). The mean diameter is chosen similar to that of Toyoura sand [20]. A rolling resistance liner model is used for the contacts. Given a grain g and one of its neighbouring grain b. They have position vectors ( | is positive, they are in contact (the contact point is denoted as c) and grain g experiences a force f gc = f n,gc + f t,gc and a rolling resistance moment M gc from grain b. The contact forces are also divided into an elastic part and a viscous part ( f * ,gc = f * e,gc + f * v,gc ). They are calculated as: where n gc is the contact normal on grain g at contact point c, m * = m g m b ∕(m g + m b ) is the effective mass, k n,t and n,t are the elastic stiffness and viscous damping constants, Δv t,c is the relative shear velocity at the contact c and Δ b,c is the relative bend rotation. Equations 2 and 3 are in an incremental style. The normal stiffness k n is calculated as k n = [ ∕2]ED * , where E is the elastic modulus and is the effective diameter. The ratio between the normal and tangential stiffness = k n ∕k t is fixed for all contacts. The rolling stiffness is calculated as k r = k t (D * ∕2) 2 . Both the tangential contact force and the rolling resistance moment have limits, which are imposed as |f t,gc | ⩽ |f n,gc | and |M gc | ⩽ r D * |f n,gc |∕2 , respectively. The (2) rolling moment is introduced to model some physical phenomena existing at contacts [21] such as the uneven distribution of contact pressure, plastic deformation around contacts, surface adhesion, etc. Most importantly, it is included to compensate for the lost rolling resistance due to the use of spherical grains in DEM.
The density of quartz (a very common sand mineral) is 2650 kg/m 3 and this value is used for grain . The viscous damping constants ( n and s ) are chosen as 0.2, which is equivalent to a restitution coefficient of about 0.5. A relatively small Young's modulus ( E = 0.15 GPa ) [22] is used such that a greater time step ( Δt ∼ √ m E ) can be used and simulations are speeded up slightly. It is not recommended to have greater than 0.5 [23], therefore, 0.4 is used in this study. Regarding the rolling resistance coefficient ( r ), although a higher value can lead to higher critical state strength that is comparable to real sand, a small value 0.1 is used here because the higher critical state strength of real sands is due to other unconsidered physical effects. In summary, the DEM model parameters are listed in Table 1.
All the simulations are performed in 3D periodic domains (explained by Thornton [24]) without gravity. Periodic simulations have an advantage over wall-controlled simulations in that specimens are homogeneous over large strain scales. In wall-controlled simulations, close to the rigid walls, the void ratio can be larger than that far away from walls. Additionally, the contact normals between grains and rigid walls constitute a large portion of whole contacts due to the limited number of grains used in DEM and they will always be perpendicular to the wall.
Specimens are generated by randomly inserting grains within a cuboidal domain (each side is 4 mm long) with the possibility of overlap until a target void ratio is achieved. Then, the domain is enforced to have no deformation, contacts are created and specimens are left to reach a stable state in two steps. In the first step, different combination of and r is used to have various specimens and in the second step, and r are fixed as in Table 1. Specimens with a variety of initial densities (density is defined relative to the critical state line in this paper) can be obtained. Depending on the void ratio, the number of grains in a specimen ranges from 13,800 to 15,000.
Several axisymmetric loading paths considered in this study are: (a) constant volume triaxial compression (CV), which approximates the conventional undrained triaxial tests in laboratory; (b) constant radial stress triaxial compression (CR), which approximates the conventional drained triaxial tests; (c) constant mean stress triaxial compression (CP), which keeps the mean stress constant and compresses the specimen in the axial direction; (d) isotropic compression (ISOC), which compresses the specimen at identical rates in all three directions and (e) isotropic dilation (ISOD). During simulations, the strain ( ) and strain rate ( ̇ ) are estimated from the size and deformation of the periodic domain and the stress is estimated from the grains by the Love's equation as in the literature [17,25]: Here, V is the volume of RVE, the contact vector l gc is the vector pointing from the centre of gravity of grain g to the contact point c. The summation is performed over all contacts ( c ∈ g ) of every grain ( g ∈ V ) in RVE. Wensrich [26] argued that when contact moment is introduced, the contact vector l gc should be shifted by an eccentricity vector calculated from the moment and the normal contact force. However, his analysis is based on the assumption that contact moment arises solely from uneven distribution of contact pressure. However, in most DEM simulations, the contribution of this source to the contact moment is negligible. And contact moment is introduced mostly to model the missing rolling resistance due to the use of spherical grains. We have also conducted several wall-controlled triaxial simulations where the stress could be estimated directly from the wall, it is found that the estimation of Eq. 4 is very accurate even with contact moment.
In the present study, compressive stress and strain are defined as positive. The z axis is in the axial loading direction, therefore, z is the axial strain. the deviatoric stress is q = 1 − 3 , the mean stress is p = ( 1 + 2 + 3 )∕3 and the void ratio is denoted as e.
In stress-controlled simulations (e.g. CR and CP tests), a servo-control mechanism is used and the strain rate is estimated from the following equation if a target stress t is expected.
where * can be x, y or z and indicates the direction, c * is the current stress in direction * , A * is the face area of the cuboidal domain. The stiffness K * is estimated from the stiffness of all contacts.
where N c is the number of contacts in RVE and n * is the directional unit vector. For the smooth control of the specimen, the strain rate is actually a weighted average of the value calculated from Eq. 5 at current step and at previous step (i.e. ̇ * = 0.25[̇ * ] t + 0.75[̇ * ] t−Δt ). All simulations are conducted with p below 3 MPa such that the average penetration rate of contacts ( c ∕D ) is smaller than 2% and the small-overlap assumption of DEM is not violated. Also, the grains are far from breakage such that crushability does not enter into play [27]. da Cruz et al. [28] suggested that the quasi-static regime is achieved when the inertia number I = |̇|D( grain ∕p) 1∕2 is smaller than 10 −3 and the norm of a tensor is its Euclidean norm. In the present shear tests, the z direction is displaced at a constant strain rate of ̇z = 10 s −1 , therefore, even at the lowest p of 0.1 MPa, I is smaller than 3 × 10 −4 and quasi-static deformation condition is met. Figure 1 presents the e-p path of some shear tests. It could be seen that the initial density of specimens spans a wide range, they all converge towards a unique critical state e-p line after undergoing either kind (CV, CR or CP) of shearing, which agrees with laboratory observations of real sands [20]. However, because spherical grains are used in DEM, the void ratios are as expected smaller than that of real sands. One thing to note is that in CV simulations, the specimens are sheared at strain rate with zero volumetric term. For some CV tests where p changes dramatically (For example, tests with e 0 = 0.606 in Fig. 1), the volume of the specimen does not change, but the void ratio changes slightly because when the void ratio is estimated in this study, we account for the volume of the small overlap between contacting grains. And this volume undergoes noticeable changes when p changes dramatically, which leads to the small variation of e. The critical state strength is found to be M = 0.92 . However, most real sands have M greater than 1. For example, M = 1.25 for Toyoura sand [20]. This is because some other physical effects are not considered in the highly idealised DEM model such as the irregular shape effect of grains, surface adhesion at contacts, etc. and these effects can contribute to the critical state strength as well. It is also possible that M in DEM is comparable to that of real sand if higher values of and r are used. But in this case, the higher values of and r are "unrealistic" and we prefer to use the "realistic" DEM parameters here.

Normalised contact force
The SFF equation relates the stress of granular material to some statistics of grain-scale entities which have clear physical meanings. The derivation of the SFF equation is extensively discussed in the literature [17,25]. Firstly, The summation in Eq. 4 can be grouped within different contact normal ( n ) directions as in the following equation. Please note that two grains are in contact at one point, l gc and f gc are different on the two grains and therefore two different contacts are considered in the present study.
where Ω is the solid angle and it represents all possible orientations of the contact normal n . Denote E PDF (n) as the probability density function of contact normal, it has ∫ n∈Ω E PDF (n)dΩ = 1 and there are N c (n) = N c E PDF (n)dΩ contacts whose normal is in direction n . Denote the average of l gc ⊗ f gc for contacts in direction n as l ⊗ f (n) , Eq. 7 can be written as Subsequently, distributions of l(n) and f (n) are assumed uncorrelated and therefore l ⊗ f (n) = l(n) ⊗ f (n) . In the next step, the directional functions ( E PDF (n) , l(n) and f (n) ) are approximated by the Fourier-Laplace series [17]. For example, E PDF (n) can be approximated by the following second-order approximation in the present study where the fabric tensor is a second-order deviatoric tensor measuring the anisotropy of contact normals. It is calculated as a statistic of the contact normals in the RVE using the following equation.
where is the isotropic unit tensor. f (n) can be approximated by the following second-order approximation [17].
where f n is the average normal force over all directions, n and t are deviatoric tensors for normal and tangential force, respectively. l(n) ≈ Dn∕2 . The last step is to substitute all the Fourier-Laplace approximations into Eq. 8 and the following SFF equation is obtained if only lower-order terms are considered.
For a polydisperse granular RVE whose diameter of grains ranges from D 0 to D N d . The grains can then be grouped into N d groups, where the diameter for grains in group k is between D k−1 and D k . The summation in Eq. 4 can then be grouped not only based on n , but also based on the diameter D.
The inner is the summation of l gc ⊗ f gc for the N c (n, k) contacts whose corresponding grain has diameter between D k−1 and D k and the contact normal is in direction n . Taking a similar procedure as the one from Eqs. 7 to 8, one has where N c (k) is number of contacts on grains in group k. Therefore, the contact network is partitioned into N d distinct subsets, where the union of all N d subsets form the complete contact network. The stress tensor is contributed by Similarly, an SFF equation for the polydisperse granular RVE is obtained: However, Eq. 14 can be re-arranged and re-interpreted.
is the total volume of grains. For each subset, V g (k) = N g (k) D 3 k ∕6 . Define Z = N c ∕N g as the coordination number and also Z(k) = N c (k)∕N g (k) . With these definitions, Eq. 14 is expressed as: Different from Eq. 15, this equation means that the total stress tensor is contributed by V (k) of each subset with a weight of Z(k) Z is a ratio of coordination numbers, Here, some properties of spherical grains are used such as the contact vector is always in the direction of contact normal l gc = � (D∕2)n gc . Because two contacting grains have a small overlap, the magnitude of contact vector is slightly smaller than the radius and ′ here (slightly smaller than 1) is to correct this. A normalised contact force can be defined f gc = 3 � f gc ∕( D 2 ) and Eq. 17 is expressed as A DEM simulation is conducted to investigate N (k) and V (k) of each subset. During a CV test of a dense specimen ( p 0 = 0.21 MPa and e = 0.606 ), the grains are divided into 5 groups based on their diameters and both N (k) and V (k) are calculated and presented in Fig. 2. For comparison, N (k) and V (k) are normalised with . It can be seen that coarse grains (larger D) have larger N (k) and their N (k) could be several times larger than that of fine grains. However, Fig. 2 shows that V (k) of different groups is approximately the same. To better explain this, a 2D example is given in Fig. 3. A specimen (Fig. 3a) is sustaining an external deviatoric stress with zz ∕ xx = 2.0 . The specimen is stratified and the grains in upper stratum have greater diameter than those V g ∫ n∈Ω n ⊗f (n, k)E PDF (n, k)dΩ below. Figure 3b gives the contact forces between grains and warmer colour means larger contact forces. Therefore, the average contact force in the upper stratum is larger than that in the bottom stratum and strong contacts concentrate in the upper part of specimen, which corresponds to the fact that coarse grains have larger N (k) . Figure 3c gives the normalised contact forces and it can be seen that contacts with larger f scatter around the whole specimen, but primarily oriented along z axis (the direction of major principal stress), which corresponds to the fact that grains of different groups have similar V (k) . It is extensively addressed in the literature that strong contacts are the dominant source in sustaining the deviatoric stress [25,29,30]. For the case in Fig. 3, strong contacts concentrate in the upper part, does it mean that the stress in the bottom part is less deviatoric? This is clearly not the case and the evaluated stress tensor in both the top and bottom parts is actually the same with zz ∕ xx = 2.0 . Therefore, in the analysis of polydisperse granular assemblies, strong contacts should be the contacts with larger normalised contact force.
Another benefit of using normalised contact force is that because V (k) of different groups is approximately the same , the summation and weight in Eq. 19 can be ignored Therefore, there is no need to study the SFF of each subset as in Eq. 16.
The normalised contact force is actually closely related to the penetration rate of contacts. The normal force modelled in this study is f n,gc = − ED * gc n gc ∕2 if the viscous term is ignored. Therefore, the normalised normal contact force is f n,gc = −(3 ∕2)E( gc ∕D * )n gc , which means that it is proportional to the elastic modulus of the material and the penetration rate of contact ( gc ∕D * ). Generally, = −(4 ∕ )E( gc ∕D * ) 3∕2 n gc . The only difference is that this relationship reveals a power proportionality. The coarse and fine grains have similar normalised contact forces and therefore, they have similar penetration rate ( gc ∕D * ). In the present study, a new notation ( gc =f gc ∕E ) is introduced and it is also a normalised contact force. However, it is dimensionless and it implies penetration rate of contacts in linear contact models ( gc ∕D g is smaller than 2% in the present study). Accordingly, Eq. 20 is

Modified stress-force-fabric relationship
Subsequently, the techniques used in deriving the SFF equation [17,19] can be similarly applied, but here a slightly different procedure is used and the uncorrelated assumption between l(n) and f (n) is not necessary any more. Firstly, the normalised contact force can be decomposed into a normal and tangential component ( = −Λ n n + t =f n ∕E +f t ∕E ) and the following equation is obtained.
Here, Λ n (n) is the average normal component of gc for contacts in direction n and (n) = − n ⊗ t (n) is the average of the outer product between n and the tangential component of gc . Λ n (n) is approximated with a second-order series.
where Λ is a scalar measuring the average of normalised normal contact force and is a second-order deviatoric tensor measuring the anisotropy of normalised normal contact force. They are both calculated as statistics of contact normals and normalised normal contact forces in the RVE as (n) is approximated with a zero-order series.
As shown by He et al. [19], t measures the overall mobilisation of contacts. It is also a deviatoric tensor and is calculated as a statistic of contact normals and normalised tangential contact forces.
The last step is to substitute all the Fourier-Laplace approximations into Eq. 22 and the following SFF equation is obtained.
This SFF equation relates the stress to some statistics of grain-scale entities (e.g. Z, Λ , , , t , etc.) which have clear physical meanings. Ignore the higher-order coupled terms Fig. 3 An example to illustrate the difference between contact force and normalised contact force involving and (relatively smaller than other terms), the mean stress p, the deviatoric stress tensor = − p and the stress ratio ∕p can be written as This means that the mean stress p of granular material is proportional linearly to the elastic modulus E of the grain, the coordination number Z, a statistic Λ measuring the average of normalised normal contact force and the solid fraction = 1∕(1 + e) . The stress ratio ∕p is related to three dimensionless deviatoric tensors ( , and t ). In other words, the strength is contributed by three sources: the anisotropy of contact structure , the anisotropy of normal contact force and the mobilisation of contacts t . The SFF equation is from rigorous derivation and some "edge" cases can also be inferred form it. Firstly, = if (a) E = 0 , which means that the grains are so soft to sustain any external load, the stress of the assembly can therefore only be zero; or (b) the coordination number Z is zero, which is only possible when the grains in an RVE are not in contact at all. In this case, the RVE is not sustaining any external load and the stress is zero; or (c) the void ratio is infinity, which means that in a RVE, there is no grains, the stress is of course zero; or (d) the average normalised contact force Λ is zero, which means that although the grains in an RVE can have geometric contacts, the contact forces are zero. In this case, due to the balance of forces for boundary grains, the external forces are also zero and the stress is zero. The stress is isotropic if 2 5 + 2 5 + 3 t = . One trivial case is that all the deviatoric stress tensors are zero. Another possibility is that the contact structure and the normal contact forces are both anisotropic and t is zero, but 2 5 + 2 5 = . This corresponds to the consolidation of specimen with initial fabric anisotropy.
To verify the SFF equation, the "force" and "fabric" terms are calculated from DEM simulations and inserted into the SFF equation and results are shown in Figs. 4, 5, 6, 7 and 8 with dashed lines. It could be seen that all the SFF predictions are extremely close to the results from DEM simulations. Better results can also be obtained if higherorder terms in the Fourier-Laplace are used. For example, Sufian et al. [25] used fourth-order approximation for E PDF (n) instead of two.

Constitutive model and discussion
The major aim of this paper and the companion paper is to build a constitutive model around the SFF equation. If all the variables in the right hand side of Eq. C0 could be appropriately modelled by several equations of incremental type, a constitutive model is immediately obtained. In this approach, more state variables are introduced and also more equations are required. However, as a result, the state of granular materials is also more precisely defined. The constitutive equations regarding the response of the contact structure and contact forces are explained in detail in the companion paper and they are listed here  Table 2. There are overwhelming number of parameters in the model. This is partly due to the complexity of granular materials, which also makes the study of them constantly a popular topic. For example, at least five model parameters in the evolution equations of Z or Λ are to account for the rate under volumetric deformations. However, if we are interested only in the response under shear deformation other than the ISOC or ISOD deformations, a single parameter is found enough. Additionally, it is highly possible that some parameters are not specific for a single granular material, but applicable to all granular materials and therefore, it is not necessary to calibrate them. For example, the equation for has several model parameters. The mechanism of contact number change is the same for all granular materials, so do we really need to have different parameters for this equation for different granular materials? These questions are not within the primary aim of the present study and therefore are saved for later research.
Here, Z com (e) = 4 + c com (e cmax − e) cz characterise a compression line in the e − Z space.

Hypoplastic nature
This model could be essentially categorised into the family of hypoplastic models [1]. Firstly, within this model there is no need to define the loading-unloading criteria or to decompose the total strain rate into elastic and plastic parts. Additionally, taking a total differential of the SFF equation and ignoring the higher-order coupled terms, the following incremental equation regarding the stress is obtained.
Considering the format of the constitutive equations (Eqs. C1-C5), Eq. 33 could be conceptually expressed as Therefore, this model can be seen as an extension to the hypoplastic models [1] by adding more state variables. e, Z and Λ appear in both linear and non-linear terms and they are mostly to measure the compaction of the materials. In the literature, a number of similar parameters (such as the state parameter in Dafalias [2]) have been used for the same purpose. The anisotropy of contact structure and contact force only enter into the non-linear term, which suggests that for the extension of existing hypoplastic models in literature to account for the fabric effect, a fabric tensor should be regarded as an internal variable and added into

Performance
Model predictions are compared with virtual experimental results in Figs. 4, 5, 6, 7 and 8. The stress and strain obtained from DEM tests are plotted in solid lines. For both the dashed lines and dotted lines, the stress is obtained by inserting the "force" and "fabric" terms into the SFF equation. The difference is that in the verification of SFF equation, both these terms are calculated as statistics of grainscale entities, but in the constitutive model, they are from evolution equations. Figure 4 gives the results of some CV tests at different void ratios. Some key phenomena observed in undrained triaxial tests of real sand [20] can also be found in DEM simulations. For example, specimens ultimately reach the same p-q line ( M = 0.93 ) irrespective of the void ratio or initial confining pressure. The stiffness and the peak strength is largely affected by the initial state ( Denser specimens have higher stiffness and peak strength). Additionally, in the p-q diagram, the hook-type response of medium dense samples (p decreases initially before q∕p = M and then, increase along the q∕p = M line) can also be observed as shown in the inset of Fig. 4a and the rightmost line in Fig. 4e. However, this effect is far less profound than that in real sands. As shown in the companion paper, the coordination number Z has initially a decreasing trend in CV tests and this decrease of coordination number can lead to the decrease of p. Due to the use of spherical grains in DEM, Z is smaller in DEM than in real sands and the decreasing trend is also less profound, which can probably explain the insignificance of the hook-type response. Figure 5 presents the results of some CR tests at different confining pressure. Similar to real sands [20], depending on the initial state, specimens may have different stiffness and peak strength. They all reach the critical state at large shear strain. Figure 6 presents the results of some CP tests, which have similar phenomena to the CR tests. However, the initial drop of void ratio for dense specimens are not clearly observed in CP tests. Because the critical state is a reference state in all the equations (Eqs. C1-C5), as long as the critical state lines are correctly captured by the fittings (e.g. Z c (e) = 4 + c cZ (e cmax − e) cZ and Λ c (e) = c cΛ (e cmax − e) ), the model is able to correctly predict the critical state for all monotonous shear tests (Figs. 4, 5, 6). Figure 7 presents some ISOC or ISOD tests in e-p plots. Similar to findings in laboratory oedemeter tests, the granular material is found to be stiffer in compressibility after been compressed. The compressbility depends on the void ratio and the confining pressure. Figure 8 presents some CP cyclic tests. In the cyclic program, p is kept constant at 0.5 MPa and the shear direction is reversed at a greater axial strain than that of the last loop (e.g. at 1%, − 1%, 2%, − 2%, 3%, − 3% and 4% axial strain). Figure 8b is the response of q for a loose specimen ( e 0 = 0.77 ) and Fig. 8d is for a dense specimen ( e 0 = 0.60 ). Figure 8f is the response of e for them. For comparison, the monotonous shear at different initial densities is also illustrated. The phenomena observed in laboratory cyclic shear tests [31] are also found in DEM simulations such as the different response of specimens with different initial void ratios, the large compressive trend when shear direction is reversed, etc.
From Eq. 30, the stress ratio ∕p is only related to the deviatoric tensors. Because these tensors are reasonably modelled and predicted by Eqs. C1, C3 and C5 as shown in the companion paper, the prediction of the stress ratio is also good in both monotonous shear tests (Figs. 4,5,6) and cyclic shear tests (Fig. 8).
The mean stress p under isotropic deformations is correctly predicted (Fig. 7) except for a relatively larger discrepancy for the compression test of a very dense specimen and the dilation test of a very loose specimen. In the stress-controlled tests (CR and CP tests in Figs. 5, 6), the specimen is predicted to contract or dilate correctly to the critical state void ratio, but the predicted void ratio does not fit the DEM results very well. This is even worse in cyclic tests. Although the model can predict the contraction of a loose specimen and the dilation of a dense specimen under cyclic loading, the predicted void ratio deviates the DEM results considerably. The reason could be that when there is accompanying shear deformations, the variation of Z and Λ under volumetric deformation does not follow exactly the observations found in ISOC and ISOD tests. An important feature of granular material is the dilatancy. Although dilatancy does not appear in the present constituve eqautions, from the simulations above, the model is able to reproduce some effects related to it. For example, the hook-type response in CV tests. This is because, for the averge terms (e.g. Z and Λ ), their response to volumetric and deviatoric deformations is modelled seperatelly.

Implications and limitations
It can be shown that, to some extend, this model is able to show the different response of specimens with initial contact structure anisotropy. For example, an initial fabric anisotropy with a triaxial structure can be denoted as 0 = a 1 n 1 ⊗ n 1 − (a 1 ∕2)n 2 ⊗ n 2 − (a 1 ∕2)n 2 ⊗ n 2 , where a 0 1 is not zero and it measures the magnitude of the initial anisotropy. n 1 , n 2 and n 3 are the eigenvectors. Specimens are consolidated to an isotropic stress condition before shearing. t is initially zero. From Eq. C0, the normal contact forces must initially be 0 = − 0 [ + 4 7 0 ] −1 to have an isotropic stress tensor for the specimen. Considering a specimen with 0 and a 1 = 0.3 but different loading directions are applied. Denote the angle between n 1 and z axis as ( cos = n 1 ⋅ n z ). Figure 9 gives some model simulations with variable for a loose specimen. The simulations show that when the axial loading direction is along the major principal direction of | 0 ( = 0 • ), the response is predicted to the stiffest, which also agrees with experimental observations [8,9]. However, the predicted influence of fabric is less profound than that observed in real sands. For example, it is shown in CR tests [8] that, the specimen can behave like very dense sand when = 0 • , but like very loose sand when = 90 • . Also, in Yoshimine et al. [14]'s undrained experiments of the same specimen, liquefaction may be observed when sheared in some direction, but in other directions, the effective stress ends up with a very big value. As mentioned several times in the manuscript, the DEM model is highly idealised and it is still far from capturing the true response of real sands. For the case of fabric anisotropy, the real sand grains are irregular and the orientation of grains is a major source of fabric anisotropy in influence the response, but is not considered in either DEM model or the constitutive model. Another problem is with the constitutive equation of . In the analysis in the companion paper, we conclude that the contacts forces are important in modelling , but only a isotropic term Λ is included. However, in contact creation, disintegration or rotation, the local contacts forces may not be fully characterised by the average Λ . Therefore, may be the missing element and the extension of the present model to account for the fabric anisotropy requires further investigation.
One interesting thing to note is that even for initially isotropic specimens ( 0 = , 0 = , t0 = ), the state is not completely known by defining only the confining pressure and the void ratio in our model. In other words, different specimens can be made by varying the combination of Z and Λ to have the same confining pressure and void ratio. Figure 10 shows some results, where different Z 0 is chosen and Λ 0 is changed accordingly to have specimens at the the same initial confining pressure and void ratio. The range of Z 0 is between 4.2 and 6, which is well in the reasonable range. The difference is negligible, but also noticeable. Therefore, in geotechnical application, where only the shear deformation is important, it is not necessary to know the coordination number due to its limited influence. However, in terms of the grain-scale mechanism, it suggests a method to test the hypothesis whether the evolution of Z and Λ follows different equations by observe the response of specimen at the same confining pressure and void ratio, but with varying Z 0 .

Towards experimental validation
The present evolution equations for all the "force" and "fabric" terms are summarised from observations of DEM simulations and parameters are also calibrated from DEM data.
With the development of laboratory technique, the data may be directly obtained on real granular materials in the future and the models maybe be verified thereafter. The "fabric" terms are relatively easier to measure. For example, as early as 1970s, Oda [8] was able to measure the fabric by freezing the soil specimen and examine it under microscopes. Recent advancement includes the X-ray computed tomography (CT) technology [32]. In terms of the "force" terms, photoelastic techniques [33] are used to make measurements of the forces within idealised granular materials. In terms of parameter calibration, in addition to the traditional measurement of stress and void ratio in experiments, the minimum additional requirement is the "fabric" terms such as Z and . Then, an equivalent"force" term can be inverted from the SFF equation and the model parameters are calibrated.

Conclusions
The primary aim of this paper and the companion paper is to build a constitutive model for granular materials with evolving contact structure and contact forces, where the contact structure and contact forces are characterised by some statistics of contact normals and contact forces. And these statistics are actually the "fabric" or "force" terms in a modified SFF equation.
The verification of the modified SFF equation and the acquire of data regarding the evolving of these statistics under various loading conditions are through DEM simulations. In the present DEM model, the granular material is modelled as assemblies of spherical grains, and a rolling resistance linear contact model is adopted. The DEM simulations are conducted under quasi-static condition which is checked by the inertia number. Also, the simulations are all in a stress level where the small overlap assumption of contacts is not violated. The axisymmetric loading paths considered in the present study include constant volume triaxial compression, constant radial stress triaxial compression, constant mean stress triaxial compression, isotropic compression and isotropic dilation.
In the analysis of the SFF equation, we have addressed that it is more appropriate to use a normalised contact force for polydisperse granular assemblies. As been demonstrated, in a randomly-mixed polydisperse granular assembly sustaining an external load, coarse grains have greater average contacts forces than fine grains. But their average normalised contact forces are similar. Because both coarse and fine grains are equally sustaining the external deviatoric stress, this average normalised contact force is a better indicator of the contact forces of the whole assembly. Also, in deriving the SFF equation, this normalised contact force should be used.
This paper has demonstrated that the modified SFF equation is able to predict the stress accurately in various tests. The constitutive equations regarding the response of the contact structure and contact forces are explained in detail in the companion paper. They along with the SFF equation compose a constitutive model, which is found capable of capturing the observed phenomena correctly and predicting the mechanical response in various loading conditions. In the discussion, the model is found to be an extension to the hypoplastic models but with more state variables.