A constitutive model for granular materials with evolving contact structure and contact forces—part II: constitutive equations

Granular Matter, Feb 2019

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. In the model, the critical state is characterised by two fitting equations and three critical state parameters. A semi-mechanistic analysis is conducted about the change of the contact number and the obtained results are combined with observed phenomena in DEM virtual experiments to give the constitutive equations for the “fabric” terms. The change of fabric anisotropy is related to the strain rate, current fabric anisotropy and also contact forces. The change of coordination number is induced by two terms related to volumetric and shear deformations, and also an additional term related to the change of fabric anisotropy. The constitutive equations regarding the “force” terms are also proposed. All the “fabric” or “force” terms are modelled to tend toward their critial state value, which agrees with Li and Dafalias’s (J Eng Mech 138(3):263–275, 2012.  https://doi.org/10.1061/(ASCE)EM.1943-7889.0000324) basic philosophy in their evolution equation for the fabric tensor. These equations along with the SFF equation form a constitutive model.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10035-019-0869-7.pdf

A constitutive model for granular materials with evolving contact structure and contact forces—part II: constitutive equations

Granular Matter May 2019, 21:20 | Cite as A constitutive model for granular materials with evolving contact structure and contact forces—part II: constitutive equations AuthorsAuthors and affiliations Xuzhen HeWei WuShun Wang Open Access Original Paper First Online: 28 February 2019 163 Downloads Abstract 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. In the model, the critical state is characterised by two fitting equations and three critical state parameters. A semi-mechanistic analysis is conducted about the change of the contact number and the obtained results are combined with observed phenomena in DEM virtual experiments to give the constitutive equations for the “fabric” terms. The change of fabric anisotropy is related to the strain rate, current fabric anisotropy and also contact forces. The change of coordination number is induced by two terms related to volumetric and shear deformations, and also an additional term related to the change of fabric anisotropy. The constitutive equations regarding the “force” terms are also proposed. All the “fabric” or “force” terms are modelled to tend toward their critial state value, which agrees with Li and Dafalias’s (J Eng Mech 138(3):263–275, 2012.  https://doi.org/10.1061/(ASCE)EM.1943-7889.0000324) basic philosophy in their evolution equation for the fabric tensor. These equations along with the SFF equation form a constitutive model. KeywordsConstitutive model Stress–force–fabric relationship Hypoplasticity DEM Granular materials  List of symbols \({{\mathbf {\mathsf{{A}}}}}\) Deviatoric tensor measuring anisotropy of contact normals \({{\mathbf {\mathsf{{C}}}}}\) Directional coordination number D Diameter of grain E Elastic modulus of grain \(E^\text {PDF}(\varvec{n})\) Probability density function of contact normal e Void ratio \(\varvec{f}^{gc}\) Contact force on grain g at contact point c \(\varvec{f}^{n,gc}\) Normal contact force \(\varvec{f}^{t,gc}\) Shear contact force \(\varvec{{\widetilde{f}}}^{gc}\) Normalised contact force \({{\mathbf {\mathsf{{F}}}}}\) Deviatoric tensor measuring anisotropy of normal contact forces \({{\mathbf {\mathsf{{F}}}}}^t\) Deviatoric tensor measuring mobilisation of contacts I Inertia number \(k_{n,t,r}\) Elastic stiffness of contacts \(m^{g}\) Mass of grain g \(\varvec{M}^{gc}\) Contact moment on grain g at contact point c \(\varvec{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 \(\varvec{r}^{g}\) Position vector of grain g \(\varvec{u}^{t,gc}\) Relative shear displacement V Volume Z Coordination number \(\beta _{n,t}\) Viscous damping constants \(\delta ^{gc}\) Penetration depth at contact c \(\varepsilon _z\) Axial strain in triaxial settings \({\dot{\varepsilon }}_v\) Volumetric strain rate \(\dot{\varvec{\varepsilon }}\) Strain rate tensor \(\dot{\varvec{\epsilon }}\) Deviatoric strain rate tensor \(\theta ^{b,gc}\) Relative bend rotation \(\kappa \) \(k_n/k_s\) \(\Lambda \) A statistic measuring the average normalised normal contact forces \({\overline{\Lambda }}^n(\varvec{n})\) Average normalised normal contact force in direction \(\varvec{n}\) \(\mu \) Friction coefficient at contacts \(\mu _r\) Rolling resistance coefficient at contact \(\rho _\text {grain}\) Density of grains \(\varvec{\sigma }\) Stress tensor \(\Omega \) Solid angle 1 Introduction In micro-mechanics, Rothenburg and Bathurst [15] 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). Equation  relating the stress of granular material to some state variables with clear physical meanings is explained in detail in the companion paper [8]. $$\begin{aligned}&\varvec{\sigma }= E\frac{Z \Lambda }{1+e}\left[ \frac{1}{3} {{\mathbf {\mathsf{{I}}}}} + \frac{2}{15}{{\mathbf {\mathsf{{A}}}}} + \frac{2}{15}{{\mathbf {\mathsf{{F}}}}}+\frac{2}{105}\text {tr}({{\mathbf {\mathsf{{F}}}}}{{\mathbf {\mathsf{{A}}}}}){{\mathbf {\mathsf{{I}}}}}+\frac{8}{105}{{\mathbf {\mathsf{{F}}}}}{{\mathbf {\mathsf{{A}}}}}+{{\mathbf {\mathsf{{F}}}}}^t \right] \end{aligned}$$ (C0) Z is the coordination number. The fabric tensor \({{\mathbf {\mathsf{{A}}}}}\) is a second-order deviatoric tensor measuring the anisotropy of contact normals and the probability density function of contact normal \(E^{\text {PDF}}(\varvec{n})\) can be approximately calculated as $$\begin{aligned} E^{\text {PDF}}(\varvec{n}) = \frac{1}{4\pi }[1 + (\varvec{n})^T{{\mathbf {\mathsf{{A}}}}}\varvec{n}] \end{aligned}$$ (1) where \(\varvec{n}\) is the directional unit vector. \(\Lambda \) is a scalar measuring the average normalised normal contact force. \({{\mathbf {\mathsf{{F}}}}}\) is a second-order deviatoric tensor measuring the anisotropy of normalised normal contact force and \({{\mathbf {\mathsf{{F}}}}}^t\) is a second-order deviatoric tensor measuring the mobilisation of contacts. Z, \(\Lambda \), \({{\mathbf {\mathsf{{A}}}}}\) , \({{\mathbf {\mathsf{{F}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}^t\) are all calculated as statistics of grain-scale entities as shown in the companion paper. E is the elastic modulus of grains and \({{\mathbf {\mathsf{{I}}}}}\) is the isotropic unit tensor. It is demonstrated that this SFF equation is able to predict the stress accurately. Ignore the higher-order coupled terms involving \({{\mathbf {\mathsf{{A}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}\) (relatively smaller than other terms), the mean stress p, the deviatoric stress tensor \({{\mathbf {\mathsf{{s}}}}} = \varvec{\sigma } - p{{\mathbf {\mathsf{{I}}}}}\) and the stress ratio \({{\mathbf {\mathsf{{s}}}}}/p\) can be written as $$\begin{aligned} p= & {} \frac{1}{3}E\frac{Z \Lambda }{1+e} \end{aligned}$$ (2) $$\begin{aligned} {{\mathbf {\mathsf{{s}}}}}= & {} E\frac{Z \Lambda }{1+e}\left[ \frac{2}{15}{{\mathbf {\mathsf{{A}}}}} + \frac{2}{15}{{\mathbf {\mathsf{{F}}}}}+ {{\mathbf {\mathsf{{F}}}}}^t \right] \end{aligned}$$ (3) $$\begin{aligned} \frac{{{\mathbf {\mathsf{{s}}}}}}{p}= & {} \frac{2}{5}{{\mathbf {\mathsf{{A}}}}} + \frac{2}{5}{{\mathbf {\mathsf{{F}}}}}+ 3{{\mathbf {\mathsf{{F}}}}}^t \end{aligned}$$ (4) This means that the mean stress p of granular material is proportional to the elastic modulus E of the grain, the coordination number Z, a statistic \(\Lambda \) measuring the average normalised normal contact force and the solid fraction \(\nu = 1/(1+e)\). The strength is contributed by three sources: the anisotropy of contact structure \({{\mathbf {\mathsf{{A}}}}}\), the anisotropy of normal contact force \({{\mathbf {\mathsf{{F}}}}}\) and the mobilisation of contacts \({{\mathbf {\mathsf{{F}}}}}^t\). A number of studies about the SFF equation has been conducted in the recent 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 [4, 5, 10, 13, 18, 21, 23, 25]. Because the SFF equation is only an equation of stress and does not explicitly contains 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 better calibrated from tests of real materials. Some research has been conducted about the evolution equations. Li and Dafalias [12]’s model, a normalised fabric tensor enters the framework of elasto-plasticity as an internal variable and a rate equations of evolution is developed for it, which is based mainly on observations that this fabric tensor evolves towards and eventually is coaxial with the loading direction under shear. Therefore, the rate is proportional to the distance between the tensor and loading direction. Sun and Sundaresan [20] proposed evolution equations for both the coordination number and the fabric tensor, which are devised to yield the desired observations. Their general technique is to use both a linear and a non-linear term of strain rate in the equation and the non-linear term is introduced to ensure the achieve of critical state values. Some studies combine semi-mechanistic analyses with some assumptions and observations to give evolution equations. Rothenburg and Kruyt [16] proposed an equation for the coordination number. Any macroscopic change in the coordination number is the overall effects of two mechanisms, i.e., contact creation and contact disintegration. They assumed that there is a ceiling value for contact creation and a floor value for contact disintegration and the rate of each effect is proportional to the distance between the current value and the limiting value. In their recent study [9], a dilatancy law is studied and proposed. In their semi-mechanistic analysis, the dilatancy is attributed to two mechanisms, i.e., loop (related to the void space) deformation and topological change. Some analytical study can be studied on the deformation of example loops such that a mathematical expression is found for this effect. The other effect due to topological change is found to related to the change of coordination number. These findings are then combined to give the equation that better fits the results. Kumar and Luding [11] revealed the importance of a history-dependent jamming density in volumetric and shear deformations and proposed evolution equations for both the stress and fabric, in which the jamming density is a variable. DEM simulations are conducted to obtain the grain-scale data and the detailed description of the granular assembly of interest, contact model, parameters, simulation procedure is given the companion paper. Some tests conducted are: (a) constant volume triaxial compression (CV); (b) constant radial stress triaxial compression(CR); (c) constant mean stress triaxial compression (CP); (d) isotropic compression (ISOC) and isotropic dilation (ISOD); (e) CP cyclic shear tests. The z axis is in the axial loading direction, therefore, \(\varepsilon _z\) is the axial strain. the deviatoric stress is \(q = \sigma _1 - \sigma _3\), the mean stress is \(p = (\sigma _1 + \sigma _2 + \sigma _3)/3\) and the void ratio is denoted as e. 2 Critical state Open image in new window Fig. 1 Contact structure and contact forces at the critical state. a \(e-Z\), b \(e-\Lambda \), c \(A_{zz}-F_{zz}\) The critical state [17] is the flow state of granular materials after been continuously sheared. This state is an important reference state in most constitutive models [3, 6, 7, 24]. Li and Dafalias [12] also proposed the anisotropic critical state theory, which account for the role of anisotropic fabric at critical state. In addition to the conventional requirements in the classical critical state theory for critical values of the stress and void ratio, it includes another requirement for the critical value of the fabric anisotropy. They also proved the uniqueness and attainability of the critical state \(e-p\) line on the basis of the Gibbs condition. There are more state variables (e.g. Z, \(\Lambda \), \({{\mathbf {\mathsf{{F}}}}}\), etc.) other than the fabric tensor in our constitutive model, and their values and relative relationships at the critical state should also be found. They are also important reference state in our model. Figure 1a presents the \(e-Z\) path for shear tests of all three types. Especially, the critical states are marked with filled solid circles and a unique relationship between e and Z at the critical state is suggested. Song et al. [19] derived a phase diagram based on theoretical analyses. They showed that possible states of granular materials lie in a region in the \(e-Z\) space because another variable compactivity is missing in this pace. Our results also show that the possible states are in a region and the critical state \(e-Z\) line is a line going through this region. Song et al. [19] also showed that the possible minimal coordination number is 4 by considering the number of balance equations and unknowns. In this study, the critical state \(e-Z\) relationship is fitted by the following power equation and the theoretical result of the minimal possible number of Z is considered (Table 1). $$\begin{aligned} Z^{c}(e) = 4 + c_{cZ}(e_{cmax} - e)^{\gamma _{cZ}} \end{aligned}$$ (5) Table 1 Parameters of the constitutive model Critical state \(e_{cmax}\) \(c_{cZ}\) \(\gamma _{cZ}\) \(c_{c\Lambda }\) \(|{{\mathbf {\mathsf{{A}}}}}^c|\) \(|{{\mathbf {\mathsf{{F}}}}}^c|\) \(|{{\mathbf {\mathsf{{F}}}}}^{tc}|\) 0.76 5.57 0.7 0.098 0.45\(\sqrt{3/2}\) 0.77\(\sqrt{3/2}\) 0.045\(\sqrt{3/2}\) Contact structure \(c_{eA}\) \(\gamma _{eA}\) \(\beta _{dA}\) \(\gamma _{A}\) \(c_{eZ}\) \(\beta _{eZ}\) \(c_{Z}\) \(c_{com}\) \(c_{vZ}\) \(\beta _{vZ}\) 1.0 0.65 35 0.4 4.96 70 0.125 7.2 0.3 3 Contact force \(\gamma _{F}\) \(\gamma _{\Lambda }\) \(c_{e\Lambda }\) \(c_{v\Lambda d}\) \(c_{M}\) 0.2 0.65 100 0.075 100 Figure 1b presents the \(e-\Lambda \) path and it is suggested that in the range of stress level considered in this study, a linear relationship is a good approximation for the critical state \(e-\Lambda \) line. This line also intersects the y axis at \(e_{cmax}\), which means that \(e_{cmax}\) is a limiting void ratio. At this void ratio, the critical state flow occurs at the minimal possible coordination number 4 and the average contact forces are zero (actually impossible). The following linear equation is used to fit the line. $$\begin{aligned} \Lambda ^{c}(e) = c_{c\Lambda }(e_{cmax} - e) \end{aligned}$$ (6) During shearing, the contact structure becomes anisotropic [12, 14] and the normal contact forces become also statistically anisotropic [15]. Non-zero \({{\mathbf {\mathsf{{A}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}\) indicate anisotropic contact structure and anisotropic normal contact forces, respectively and the norms (\(|{{\mathbf {\mathsf{{A}}}}}|\) and \(|{{\mathbf {\mathsf{{F}}}}}|\)) measure the magnitude of anisotropy. Zhao and Guo [26] conducted DEM simulations and examined the fabric tensor of contact normals (i.e. \({{\mathbf {\mathsf{{A}}}}}\) in the present paper) at the axial strain between 40% and 50% (assumed to reach the critical state). They found that \(A^c = |{{\mathbf {\mathsf{{A}}}}}^c|\) is related not only to the mean stress p, but also to the loading condition (for example, the Lode angle). Dafalias [1] suggested that if suitable fabric tensor is defined, its normalised norm may have unique critical state value. Li and Dafalias  [12] also argued that the norm of a fabric tensor can depend on the Lode angle, therefore, they defined a normalised fabric tensor, which is unique at the critical state. In the present study, all the specimens are prepared to have almost isotropic inherent contact structure (i.e. \({{\mathbf {\mathsf{{A}}}}}|_0 \approx {{\mathbf {\mathsf{{0}}}}}\)). During shearing (triaxial compression or extension only in this paper), \({{\mathbf {\mathsf{{A}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}\) are both coaxial with the stress tensor. Figure 1c gives the development of \(A_{zz}\) and \(F_{zz}\) for selected shear tests and they all approach a positive value at about 25% axial strain. At about 25% axial strain, the stress and void ratio seem to have met the requirement of critical state. However, it is possible that the critical state is not actually reached because the requirement for the fabric anisotropy is not met [12]. Our results show that \(A_{zz}\) at this strain depends on the void ratio or the mean stress slightly and ranges from 0.38 to 0.52, which agrees with the studies mentioned before. However, because the range of stress level is limited and only triaxial compression or extension is considered in the presented study, we use constant \(A^c\) and \(F^c\) for simplicity. If the range of stress level increases or more complex loading path is involved, it is not appropriate to use constant value for them. The critical state \(e-p\) relationship and the soil constant \(M = (q/p)|_{critical\ state}\) are an immediate outcome of the SFF equation and Eqs. 5–6. The predicted critical state \(e-p\) relationship is found to fit the DEM simulations very well in the companion paper and M is found to be 0.93, which also agrees with DEM simulations. However, most real sands have M greater than 1. For example, \(M=1.25\) for Toyoura sand [22]. 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 \(\mu \) and \(\mu _r\) are used. But in this case, the higher values of \(\mu \) and \(\mu _r\) are “unrealistic” and we prefer to use the “realistic” DEM parameters. 3 Contact structure There have been some constitutive equations for Z and \({{\mathbf {\mathsf{{A}}}}}\) in the literature. For example, in Li and Dafalias [12]’s model, a rate equation of evolution for a normalised fabric tensor is developed. Sun and Sundaresan [20] proposed evolution equations for both the coordination number and the fabric tensor, both of the equations include a linear and a non-linear term of the strain rate. In this paper, a semi-mechanistic analysis is firstly conducted and the results are combined with observed phenomena in DEM virtual expereiments to give the constitutive equations, which are expected to model the response in various kinds of loading conditions. Similar approach has also been used in the literature of granular materials. For example, in Rothenburg and Kruyt [16]’s study for the evolution equation of the coordination number, their later study for the dilatancy law [9, 11]’s macroscopic constitutive model, in which a history-dependent jamming density is a variable. In the present study, both Z and \({{\mathbf {\mathsf{{A}}}}}\) characterise the contact structure and they can be combined together to have a new tensor \({{\mathbf {\mathsf{{C}}}}} = Z({{\mathbf {\mathsf{{I}}}}}+{{\mathbf {\mathsf{{A}}}}})\), which is a directional coordination number. Z is the mean part of \({{\mathbf {\mathsf{{C}}}}}\) and \({{\mathbf {\mathsf{{A}}}}}\) is the scaled deviatoric part of \({{\mathbf {\mathsf{{C}}}}}\). In a RVE, denote the number of grains as \(N_g\). If there are two contacts with opposite directions at a contact point, there are then totally \(ZN_g\) contacts. Given any differential solid angle \(d\Omega \) whose directional unite vector is \(\varvec{n}\), from Eq. 1, there are $$\begin{aligned}ZN_gE^{\text {PDF}}(\varvec{n}) d\Omega &= ZN_g\frac{1 + (\varvec{n})^T{{\mathbf {\mathsf{{A}}}}}\varvec{n}}{4\pi }d\Omega \nonumber \\& = N_g\frac{(\varvec{n})^T{{\mathbf {\mathsf{{C}}}}}\varvec{n}}{4\pi }d\Omega = N_gC(\varvec{n})d\Omega \end{aligned}$$ (7) contacts whose contact normals are in this direction, where \(C(\varvec{n}) = [(\varvec{n})^T{{\mathbf {\mathsf{{C}}}}}\varvec{n}]/[4\pi ]\). After the RVE undergoes a small strain \(\dot{\varvec{\varepsilon }}dt\), this counting of contacts is changed and it can be caused by three reasons: The disintegration of existing contacts. Denote the possibility of disintegration for every contact as \(P^\text {disintegration}\) and this source has a change of \(-N_gC(\varvec{n})d\Omega P^\text {disintegration}\). The creation of contacts between close but not contacting grains. For every grain in an assembly sustaining external load, it must be in contact with several neighbouring grains if it serves to transfer load. It also has some neighbours such that there are void gaps between them and there is a possibility that contacts could created between them. It could be roughly estimated that the number of neighbouring grains per grain per unit solid angle \(C^\text {neigh}(e)\) depends only on the void ratio. Therefore, the number of neighbours with a potential to create contact is \(N_g[C^\text {neigh}(e) -C(\varvec{n})]d\Omega \). Similarly, denote the possibility of creation for every of these neighbours as \(P^\text {create}\) and the total creation is \(N_g[C^\text {neigh}(e)-C(\varvec{n})]d\Omega P^\text {create}.\) The rotation of contacts. These contacts can rotate towards other directions and also contacts in other directions can rotate to direction \(\varvec{n}\). Denote the increase rate as \(P^\text {rotation}\) and this source has a change of \(N_gC(\varvec{n})d\Omega P^\text {rotation}\). In summary, a balance equation is obtained as: $$ \begin{aligned}d[N_gC(\varvec{n})d\Omega ]&= {\dot{N}}_gC(\varvec{n})d\Omega dt + N_g{\dot{C}}(\varvec{n})d\Omega dt\\&= -N_gC(\varvec{n})d\Omega P^\text {disintegration}\\&\quad +\,\,N_g[C^\text {neigh}(e)-C(\varvec{n})]d\Omega P^\text {create} +N_gC(\varvec{n})d\Omega P^\text {rotation} \end{aligned}$$ (8) After rearrangement, the rate of \(C(\varvec{n})\) has $$\begin{aligned} \begin{aligned}&{\dot{C}}(\varvec{n}) = -\frac{\dot{N_g}}{N_g}C(\varvec{n}) + C^\text {neigh}(e)P^\text {create}/dt\\&\qquad\qquad -[P^\text {disintegration}/dt +P^\text {create}/dt - P^\text {rotation}/dt]C(\varvec{n})\\ \end{aligned} \end{aligned}$$ (9) It must be firstly recognised that the change of contact number is driven by the deformation of RVE for rate-independent system. Therefore, if \(\dot{\varvec{\varepsilon }} ={{\mathbf {\mathsf{{0}}}}}\), all three sources are zero and \(P^* = 0\) (* can be create, disintegration or rotation). Also, the possibilities are highly related to the contact forces. For example, a weak contact definitely has a higher possibility to disintegrate than a strong one. Imagining the creation of contacts as pushing grains into the bores, a stronger contact network could facilitate this process and lead to a higher possibility to create contacts. \(\Lambda \) is a statistic of normalised normal contact forces and it measures the average, which is a good indicator for the contact forces. Consequently, we decompose \(P^* = \text {tr}({{\mathbf {\mathsf{{L}}}}}^*(\Lambda )\dot{\varvec{\varepsilon }}) dt+ N^*(\Lambda )|\dot{\varvec{\varepsilon }}|dt\) into two parts, where the first part is linear in \(\dot{\varvec{\varepsilon }}\) and the second part is non-linear in \(\dot{\varvec{\varepsilon }}\). Both terms are functions of \(\Lambda \). In the first term of Eq. 9, the counting of grains \(N_g\) is only related to the void ratio e if the grading is fixed. Therefore, the rate of \(N_g\) is only a function of void ratio (\(\dot{N_g}/N_g = f^{t1}(e){\dot{\varepsilon }}_v\), where \({\dot{\varepsilon }}_v\) is the volumetric strain rate). With the analysis above, Eq. 9 is expressed as the following equation after the all the terms having the same factor regarding \(\dot{\varvec{\varepsilon }}\) are organised together and coefficients are represented by general functions. $$\begin{aligned}&{\dot{C}}(\varvec{n}) = -f^{e}(\Lambda )C(\varvec{n})|\dot{\varvec{\varepsilon }}| + f^{v}(e, \Lambda , C(\varvec{n})){\dot{\varepsilon }}_v \nonumber \\&\qquad \qquad +\,\text {tr}[{{\mathbf {\mathsf{{f}}}}}^s(e, \Lambda , C(\varvec{n}))\dot{\varvec{\epsilon }}] \end{aligned}$$ (10) where \(\dot{\varvec{\epsilon }} = \dot{\varvec{\varepsilon }} -{\dot{\varepsilon }}_v{{\mathbf {\mathsf{{I}}}}} \) is the deviatoric strain rate and \(f^{e}(\Lambda ) = N^\text {disintegration}(\Lambda ) + N^\text {create}(\Lambda ) - N^\text {rotation}(\Lambda )\). Because Eq. 10 is valid for any direction \(\varvec{n}\), it is postulated that the rate of the directional coordination number \({{\mathbf {\mathsf{{C}}}}}\) has a similar format as bellow where the general functions are assumed to be isotropic in this study for the sake of simplicity. $$\begin{aligned}&\dot{{{\mathbf {\mathsf{{C}}}}}}= -\,f^{e}(\Lambda ){{\mathbf {\mathsf{{C}}}}}|\dot{\varvec{\varepsilon }}| \nonumber \\&\quad \qquad +\,f^{v}(e, \Lambda ,{{\mathbf {\mathsf{{C}}}}}){\dot{\varepsilon }}_v {{\mathbf {\mathsf{{I}}}}} + f^s(e, \Lambda , {{\mathbf {\mathsf{{C}}}}})\dot{\varvec{\epsilon }} \end{aligned}$$ (11) The critical state is an important reference state. When reaching the critical state, the material is flowing at constant volume (fixed void ratio e) with fixed contact structure \({{\mathbf {\mathsf{{C}}}}}^c(e)\) and any continuing deviatoric shear deformation will not cause the change of the statistics of the contact structure. Additionally, at this state, the contact force is \(\Lambda ^c(e)\). Therefore, the following equation must be valid. $$\begin{aligned} {{\mathbf {\mathsf{{0}}}}}= -\,f^{e}(\Lambda ^c){{\mathbf {\mathsf{{C}}}}}^c|\dot{\varvec{\varepsilon }}| + f^s(e, \Lambda ^c, {{\mathbf {\mathsf{{C}}}}}^c)\dot{\varvec{\epsilon }} \end{aligned}$$ (12) Open image in new window Fig. 2 Development of contact structure anisotropy. (solid lines are DEM results and dashed lines are model predictions). a, b CV tests. c, d CR tests. e, f CP tests, g CP Cyclic, p = 0.5 MPa, \(e_0\) = 0.77, h CP Cyclic, p = 0.5 MPa, \(e_0\) = 0.60 Subtract Eqs. 12 from 11, the following equation is obtained. $$\begin{aligned}&\dot{{{\mathbf {\mathsf{{C}}}}}}= f^{e}(\Lambda )[f^{d}(\Lambda ){{\mathbf {\mathsf{{C}}}}}^c - {{\mathbf {\mathsf{{C}}}}}]|\dot{\varvec{\varepsilon }}|\nonumber \\&\qquad +\,f^{v}(e, \Lambda ,{{\mathbf {\mathsf{{C}}}}}){\dot{\varepsilon }}_v {{\mathbf {\mathsf{{I}}}}} + f^{s'}(e, \Lambda , {{\mathbf {\mathsf{{C}}}}}) \dot{\varvec{\epsilon }} \end{aligned}$$ (13) where \(f^{d}(\Lambda ) = f^{e}(\Lambda ^c)/f^{e}(\Lambda )\), which is one when \(\Lambda = \Lambda ^c\). Because \({{\mathbf {\mathsf{{A}}}}}\) is the scaled deviatoric part of \({{\mathbf {\mathsf{{C}}}}}\), We postulate the following equation which is based on the first term of Eq. 13. $$\begin{aligned} \dot{{{\mathbf {\mathsf{{A}}}}}}= f^{e,A}|f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}|^{\gamma _A}(f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}) |\dot{\varvec{\varepsilon }}| \end{aligned}$$ (C1) This equation means that the change of \({{\mathbf {\mathsf{{A}}}}}\) is driven by the deformation. Most importantly, the rate depends on a distance (\(f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}\)), which is between the current \({{\mathbf {\mathsf{{A}}}}}\) and its critical state value \({{\mathbf {\mathsf{{A}}}}}^c\) corrected by a coefficient \(f^{d,A}\). The critical state value \({{\mathbf {\mathsf{{A}}}}}^c = A^c(e,\frac{\dot{\varvec{\epsilon }}}{|\dot{\varvec{\epsilon }}|}) \frac{\dot{\varvec{\epsilon }}}{|\dot{\varvec{\epsilon }}|}\) is assumed to coaxial with the deformation direction \(\frac{\dot{\varvec{\epsilon }}}{|\dot{\varvec{\epsilon }}|}\). And the norm \(A^c\) is related to the void ratio e and the loading direction (\(\frac{\dot{\varvec{\epsilon }}}{|\dot{\varvec{\epsilon }}|}\)). However, as mentioned that the range of stress level is limited and only triaxial compression or extension is considered in the presented study, a constant \(A^c\) is used in the present study. This term (\(|f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}|^{\gamma _A}\)) means that the rate depends on the distance non-linearly. \(f^{e,A}\) is a stiffness function. This equation shares some similar ideas to the following evolution equation from [12] : $$\begin{aligned} \dot{{{\mathbf {\mathsf{{F}}}}}}^{ld} = <\lambda > c({{\mathbf {\mathsf{{n}}}}}^{ld} - r {{\mathbf {\mathsf{{F}}}}}^{ld}) \end{aligned}$$ (14) where \({{\mathbf {\mathsf{{F}}}}}^{ld}\) is their normalised fabric tensor and it will evolve towards the tensor-valued loading direction \({{\mathbf {\mathsf{{n}}}}}^{ld}\) if subject to continuous shear. Because the tensor \(A_c\) is along the strain rate direction, which is the same as \(n^{ld}\). \(<\lambda>\) is the plastic multiplier and it coresponds to \(|\dot{\varvec{\varepsilon }}|\) in a hypoplastic formulation that does not distinguish between elastic and plastic strain rates. Also, it is easy to see that the \(f^{d,A}\) is the same as \(\frac{1}{r}\) in Eq. 14. However, they do not specify the correction r in their model and simply use \(r=1\). The development of contact structure anisotropy under various kinds of shear loading is presented in Fig. 2. The first row is some CV tests at different e, the second row is some CR tests at different \(\sigma _r\) and the third row is some CP tests at different p. Some observations are: For every test, the specimen has initially an isotropic contact structure (\({{\mathbf {\mathsf{{A}}}}}|_0 \approx {{\mathbf {\mathsf{{0}}}}}\)). After continuous shear deformation, the contact structure approaches the critical state (\({{\mathbf {\mathsf{{A}}}}} \approx {{\mathbf {\mathsf{{A}}}}}^c\)). The rate is large at the start and decreases after \({{\mathbf {\mathsf{{A}}}}}\) is close to \({{\mathbf {\mathsf{{A}}}}}^c\), which is modelled by the dependence of \(\dot{{{\mathbf {\mathsf{{A}}}}}}\) on the distance (\(f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}\)). Dense specimens show a peak before the contact structure eventually approaches the critical state. This is modelled by the coefficient \(f^{d,A}\), which is a function of \(\Lambda \). A \(f^{d,A}\) greater than one will cause the contact structure to increase to a greater anisotropy than \({{\mathbf {\mathsf{{A}}}}}^c\). For dense specimens, their \(\Lambda \) is smaller than \(\Lambda ^c\) before approaching the critical state and it is also required that \(f^{d,A} =1\) when \(\Lambda = \Lambda ^c\). Therefore, it is modelled as \(f^{d,A} = \text {exp}(\beta _{dA}[\Lambda ^c - \Lambda ])\), which is mathematically similar to the correction to the size of bounding surface in Dafalias and Li [2]’s model. To calibrate \(\beta _{dA}\), the peak anisotropy \({{\mathbf {\mathsf{{A}}}}}^p\) and the corresponding \(\Lambda ^p\) are found for all tests. It is easy to find that they follow an equation \(\text {exp}(\beta _{dA}[\Lambda ^c - \Lambda ^p]){{\mathbf {\mathsf{{A}}}}}^c={{\mathbf {\mathsf{{A}}}}}^p\) from Eq. . The stiffness (\(\partial A_{zz}/\partial \varepsilon _{zz}\) here) depends on the density and \(f^{e,A}\) (a function of \(\Lambda \)) is to account for this. From Eq. , the initial stiffness is \([\partial {{{\mathbf {\mathsf{{A}}}}}}/\partial |\varvec{\varepsilon }|]|_{0} =f^{e,A}[f^{d,A}]^{{\gamma _{A}}+1}|{{\mathbf {\mathsf{{A}}}}}^c|^{\gamma _{A}} {{\mathbf {\mathsf{{A}}}}}^c\). \(f^{d,A}\) is close to one and its influence could be ignored firstly such that \(f^{e,A}\) is proportional to the initial stiffness (\([\partial {{{\mathbf {\mathsf{{A}}}}}}/\partial |\varvec{\varepsilon }|]|_{0} = f^{e,A}|{{\mathbf {\mathsf{{A}}}}}^c|^{\gamma _{A}} {{\mathbf {\mathsf{{A}}}}}^c\)). It can be seen from the tests that dense specimens have a greater initial stiffness and also they have smaller \(\Lambda \) initially. Therefore, the stiffness function is modelled by a non-linear power equation as \(f^{e,A} = c_{eA}/\Lambda ^{\gamma _{eA}}\), where the parameters could be calibrated from the initial stiffness data and the corresponding \(\Lambda |_0\) from the relationship \([\partial {{{\mathbf {\mathsf{{A}}}}}}/\partial |\varvec{\varepsilon }|]|_{0} = c_{eA}/(\Lambda |_0)^{\gamma _{eA}}[{{\mathbf {\mathsf{{A}}}}}^c]^{\gamma _{A}} {{\mathbf {\mathsf{{A}}}}}^c\). \(f^{e,A}\) is a decreasing function of \(\Lambda \), which could be explained by that the initial development of anisotropy of contact structure is dominant by disintegration of contacts in minor principal stress direction. As explained before, a weak contact has a higher possibility to disintegrate than a strong one and therefore, \(P^\text {disintegration}\) and then \(N^\text {disintegration}(\Lambda )\) are decreasing functions of \(\Lambda \). Recapping that \(f^{e}(\Lambda ) = N^\text {disintegration}(\Lambda ) + N^\text {create}(\Lambda ) - N^\text {rotation}(\Lambda )\) and if the effect of disintegration dominates, \(f^e(\Lambda )\) is consequently a decreasing function. This effect could also explain the initial drop of Z later. \(\dot{{{\mathbf {\mathsf{{A}}}}}}\) depends on the distance (\(f^{d,A}{{\mathbf {\mathsf{{A}}}}}^c - {{\mathbf {\mathsf{{A}}}}}\)) non-linearly. The parameter \(\gamma _{A}\) is to model this. For a dense specimen, \(\gamma _{A}\) could change the strain where the peak happens. A larger \(\gamma _{A}\) will cause the peak anisotropy of contact structure happens at a larger shear strain. With Eq.  and parameters calibrated from the data, predictions of the equation are presented in Fig. 2 with dashed lines. It should be noted that before the introduction of any model for \(\Lambda \), interpolated \(\Lambda \) is used for the simulations in Fig. 2. The volumetric strain rate for CR or CP tests is also interpolated from the data between shear strain and void ratio. The results show that the equation could correctly model the development of \({{\mathbf {\mathsf{{A}}}}}\) for various specimens under various loading conditions. The error is generally small. However, because the present model ignores the small variation of \({{\mathbf {\mathsf{{A}}}}}^c\) with e, for some tests where the model parameter deviates from the actual \({{\mathbf {\mathsf{{A}}}}}^c\) slightly (Fig. 2b, e), the error close to the critical state is relatively large. Figure 2g, h also shows some results from cyclic tests. During triaxial compression, the contacts take a triaxial compression structure (\(A_{zz} > 0\)) and they gradually change to a triaxial extension structure (\(A_{zz} < 0\)) when loading direction is switched. The evolution equation for \({{\mathbf {\mathsf{{A}}}}}\) can also reproduce this phenomenon. Open image in new window Fig. 3 Evolution of coordination number (solid lines are DEM results and dashed lines are model predictions). a, b CV tests. c ISOC tests. d ISOD tests. e CR tests, f CP tests The coordination number is modelled by the following equation. $$\begin{aligned} {\dot{Z}}= f^{e,Z}(Z^c - Z)|\dot{\varvec{\epsilon }}| + f^{v,Z}{\dot{\varepsilon }}_v -c^Zf^{e,Z}|\dot{{{\mathbf {\mathsf{{A}}}}}}| \end{aligned}$$ (C2) where the first two terms are based on the first two terms of Eq. 13. The second term is to model the change of Z caused by volumetric deformation and the third term is to model the change of Z associated with the change of anisotropy of contact structure. They both vanish at the critical state and the first term serves to drive Z to its critical state value \(Z^c\) under shear deformations. The evolution of coordination number is presented in Fig. 3. The first row is some CV tests, the second row is ISOC or ISOD tests and the third row is CP or CR tests. Some observations are: In CV tests (Fig. 3a, b), the second term is always zero. Subject to continuous shear, the coordination number is evolving towards \(Z^c\) and this is modelled by the first term in Eq. . Additionally, at early stages when the anisotropy of contact structure changes dramatically, there is an associated large decreasing trend of Z. For example, please see the bottommost line in Fig. 3b. Even if its initial Z is already smaller than \(Z^c\), it will decrease first to a minimal value and then evolves towards \(Z^c\). This decreasing trend is closely related to the rate of \({{\mathbf {\mathsf{{A}}}}}\). The dense specimens have larger \(|\dot{{{\mathbf {\mathsf{{A}}}}}}|\) at early stages, the decreasing trend of Z is also larger (Fig. 3a). These phenomena are modelled by the third term in Eq. . The first and third terms have an identical stiffness function \(f^{e,Z} = c_{eZ}\text {exp}(\beta _{eZ}[\Lambda ^c - \Lambda ])\) which varies with the distance \( \Lambda ^c - \Lambda \) slightly. For calibration, the variation of \(f^{e,Z}\) could be ignored firstly and \(c^Z\) which controls the relative significance of the two terms can be determined from a loose specimen. Subsequently, \(c_{eZ}\) and \(\beta _{eZ}\) can be found by estimating the initial rate \({\dot{Z}}|_0\) from data. As explained before, the initial development of anisotropy of contact structure is dominant by disintegration of contacts at minor principal stress direction. Therefore, this also explains that accompanying with the change of contact structure anisotropy, there is an associated decreasing trend of Z. During ISOC or ISOD tests, the first and the third term are zero and \(f^{v,Z}\) controls the variation of Z with e. Fig. 3c gives some ISOC results in \(e-Z\) plot. The critical state \(e-Z\) line is indicated by a dotted line. The tests suggest that there is a compression \(e-Z\) line \(Z^{com}(e)\) (the black solid line), which is fitted by the following equation. $$\begin{aligned} Z^{com}(e) = 4 + c_{com}(e_{cmax} - e)^{\gamma _{cz}} \end{aligned}$$ (15) During compression, if the state is not on this line, it will move towards this line and if the state is already on this line, it will move along this line. Figure 3d gives the evolution of Z with e for some ISOD tests. During dilation, the rate \({\dot{Z}}\) depends on the contact force \(\Lambda \). For tests at the same initial void ratio, dense specimens have smaller \(\Lambda \) and \(f^{v,Z}\) is larger. In summary, the variation of Z with e is modelled by the following equation. $$\begin{aligned} f^{v,Z} = \left\{ \begin{array}{lll} &{}\text {exp}(\beta _{vZ}[Z^{com} -Z])\frac{dZ^{com}}{de}(1+e) &\quad {}{\dot{\varepsilon }}_v > 0\\ &{}c_{vZ}/\Lambda &\quad{} {\dot{\varepsilon }}_v < 0 \end{array}\right. \end{aligned}$$ (16) where the part \([dZ^{com}/de](1+e)\) measures the variation of Z with e along the compression \(e-Z\) line and the coefficient \(\text {exp}(\beta _{vZ}[Z^{com} -Z])\) is to drive the \(e-Z\) state towards the line if the state is not on the line. Predictions of Eq.  are presented in Fig. 3 with dashed lines. It is shown that this equation captures the phenomena detailed above and correctly predicts the coordination number for various specimens under various loading conditions. To make prediction of CR or CP tests, the volumetric strain rate is again interpolated from the data. This estimated volumetric strain rate can have large error and the variation of Z is very sensitive to the volumetric strain rate. Therefore, for CR or CP tests, only predictions at small strain are shown to illustrate that the model could capture the initial decrease or increase trend. 4 Contact force \(\Lambda \), \({{\mathbf {\mathsf{{F}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}^t\) are all statistics of contact forces. Their response under various loading conditions can be similarly studied and constitutive equations of incremental type are proposed subsequently. The development of the normal contact force anisotropy is shown in Fig. 4. It is similar to the development of contact structure anisotropy. Therefore, the following equation is proposed. $$\begin{aligned} \dot{{{\mathbf {\mathsf{{F}}}}}}= f^{e,F}|f^{d,F}{{\mathbf {\mathsf{{F}}}}}^c - {{\mathbf {\mathsf{{F}}}}}|^{\gamma _F}(f^{d,F}{{\mathbf {\mathsf{{F}}}}}^c - {{\mathbf {\mathsf{{F}}}}}) |\dot{\varvec{\varepsilon }}| \end{aligned}$$ (C3) Open image in new window Fig. 4 Development of contact force anisotropy. (solid lines are DEM results and dashed lines are model predictions). a, b CV tests. c, d CR tests. e, f CP tests From Fig. 1c, the ratio between the peak contact force anisotropy and peak contact structure anisotropy is approximately equal to the ratio of their critical state values (i.e. \(|{{\mathbf {\mathsf{{F}}}}}^p|/|{{\mathbf {\mathsf{{A}}}}}^p| \approx |{{\mathbf {\mathsf{{F}}}}}^c|/|{{\mathbf {\mathsf{{A}}}}}^c|\)), therefore, \(f^{d,F}\) is assumed to be equal to \(f^{d,A}\). Additionally, the stiffness function is found to follow an identical power relationship and \(f^{e,F} = [|{{\mathbf {\mathsf{{F}}}}}^c|/|{{\mathbf {\mathsf{{A}}}}}^c|] f^{e,A}\). In Fig. 1c, draw a line connecting the origin and the point of critical state, all paths are above this line, which means that the development of contact force anisotropy is much faster than that of contact structure. This is because the contact force is much easer to be changed than the structure. In the model, choosing a \(\gamma _F\) smaller than \(\gamma _A\) would be able to model this. With \(\gamma _F\) chosen as 0.2 and the simulated results of Eq.  are presented in Fig. 4 with dashed lines. Again, the equation captures all the features and correctly models the development of contact force anisotropy. The evolution of the average normalised normal contact force \(\Lambda \) in various tests is presented in Fig. 5 and it is modelled by the following equation. $$\begin{aligned} {\dot{\Lambda }}= f^{e,\Lambda }|{{\mathbf {\mathsf{{A}}}}}|(\Lambda ^c - \Lambda )|\dot{\varvec{\epsilon }}| + f^{v,\Lambda }{\dot{\varepsilon }}_v \end{aligned}$$ (C4) where the first term is to recognise that \(\Lambda \) is evolving towards \(\Lambda ^c\) in all three kinds of shear tests (Fig. 1b) and the second term is to model the variation of \(\Lambda \) when the specimen is subject to a volumetric deformation. Open image in new window Fig. 5 Evolution of \(\Lambda \). a, b CV tests (solid lines are DEM results and dashed lines are model predictions). c ISOC tests. d ISOD tests Firstly, focusing on CV tests (Fig. 5a, b) where \({\dot{\Lambda }}\) is only controlled by the first term. It is observed that for all tests, at early stages of shear, \({\dot{\Lambda }}\) is small and \(\Lambda \) converges to \(\Lambda ^c\) more rapidly at later stages. The reason could be that the “comfort” \(\Lambda \) for an assembly is a range whose size is related to the contact structure. When the contact structure is more or less isotropic at early stages, the “comfort” \(\Lambda \) zone is large and \(\Lambda \) does not need to change much under shear deformation. However, when the contact structure is anisotropic, the size of the “comfort” \(\Lambda \) zone shrinks and it reduces to a critical state \(e-\Lambda \) line when \({{\mathbf {\mathsf{{A}}}}}\) reaches \({{\mathbf {\mathsf{{A}}}}}^c\). Therefore, at later stages, the contact structure is anisotropic, the rate of \(\Lambda \) also increases. Mathematically, \(|{{\mathbf {\mathsf{{A}}}}}|\) in the coefficient of the first term is to model this. \({\dot{\Lambda }}\) is also related to the density. Please see the four tests shown in Fig. 5a, the loose specimens (darker lines) converges to \(\Lambda ^c\) faster. Similarly, the density can be indicated by the distance \(\Lambda ^c - \Lambda \). Therefore, \(f^{e,\Lambda }\) is modelled as \(f^{e,\Lambda } = c_{e\Lambda }\Lambda ^{\gamma _{\Lambda }}f^{e,Z}\) in this paper, where the dependence on density is assumed to be related to the stiffness function of Z (i.e. \(f^{e,Z}\)). Figure 5c gives the \(e-\Lambda \) path for ISOC tests. Under isotropic compression, \(\Lambda \) is moving parallel to the critical state \(e-\Lambda \) line for all tests. Figure 5d shows that \(\Lambda \) is moving parallel to another line (\(\Delta \Lambda / \Delta e\) = 0.075 = \(c_{v\Lambda d}\)) in ISOD tests. Therefore, \(f^{v,\Lambda }\) is modelled by the following equation. $$\begin{aligned} f^{v,\Lambda } = \left\{ \begin{array}{lll} &{}c_{c\Lambda }(1+e)&\quad {}{\dot{\varepsilon }}_v > 0\\ &{}c_{v\Lambda d}(1+e)&\quad {}{\dot{\varepsilon }}_v < 0 \end{array}\right. \end{aligned}$$ (17) It is found that the response of shear contact forces is extremely fast and \({{\mathbf {\mathsf{{F}}}}}^t\) reaches a level comparable to its critical state value in a very small shear strain (less than 2\(\%\) in Fig. 6). Even if the response of \({{\mathbf {\mathsf{{F}}}}}^t\) could be complex (Fig. 6), it is simply modelled by the following equation. $$\begin{aligned} \dot{{{\mathbf {\mathsf{{F}}}}}^t}= c_M({{\mathbf {\mathsf{{F}}}}}^{tc} - {{\mathbf {\mathsf{{F}}}}}^t) |\dot{\varvec{\varepsilon }}| \end{aligned}$$ (C5) Open image in new window Fig. 6 Mobilisation of contacts 5 Discussions The proposed constitutive equations (Eqs. –) together with the SFF equation form a constitutive model, which not only predicts the stress, strain and void ratio of the material as traditional constitutive models do, but also predicts the statistical feature of the contact structure and contact forces. Therefore, this model can readily simulate the response of material with inherent contact structure anisotropy as shown in the companion paper. The mathematical format of the constitutive equations (Eqs. –) is largely inspired by that of the simple hypoplastic equation [24]. And the present model is also a hypoplastic model because 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. Actually, as shown in the companion paper, the present model can be seen as an extension to the hypoplastic models by adding more state variables. The fabric tensor \({{\mathbf {\mathsf{{A}}}}}\) and the tensor measuring the anisotropy of normal contact forces \({{\mathbf {\mathsf{{F}}}}}\) are both modelled very well by the corresponding equations. It is also revealed that their stiffness functions are closely related, which implies that there may exist some grain-scale mechanisms that could explain this and it needs further theoretical investigations. Additionally, the model for the tensor measuring the mobilisation of contacts \({{\mathbf {\mathsf{{F}}}}}^t\) is simple in the present study. However, results show that its response is actually complex and similar to that of \({{\mathbf {\mathsf{{A}}}}}\) or \({{\mathbf {\mathsf{{F}}}}}\). It may be better modelled in the future without introducing many more parameters. Z and \(\Lambda \) show totally different response under compression or dilation deformation and different equations are also used (Eqs. 16 and 17) depending on whether \({\dot{\varepsilon }}_v\) is positive or negative. However, it is possible that when there is shear deformation along with volumetric deformation, the response of Z and \(\Lambda \) can be different from that observed in isotropic volumetric deformations and this also needs further investigation. 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 are in the evolution equations of Z or \(\Lambda \) 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 \({{\mathbf {\mathsf{{A}}}}}\) 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 is saved for later research. 6 Conclusions This paper presents 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 SFF equation is proved to predict the stress accurately in various tests in the companion paper. A main contribution of the present paper is the constitutive modelling around the SFF equation. The basic philosophy in Li and Dafalias [12]’s evolution equation of fabric tensor is followed. All the “fabric” or “force” terms are modelled to tend toward their critial state value, which is characterised by a critical state e-Z relationship, a critical state e-\(\Lambda \) relationship and some critical state model parameters concerning the anisotropy of contact structure \({{\mathbf {\mathsf{{A}}}}}^c\), the anisotropy of normalised normal contact forces \({{\mathbf {\mathsf{{F}}}}}^c\) and the mobilisation of contacts \({{\mathbf {\mathsf{{F}}}}}^{tc}\). Regarding the evolution equations of the “fabric” terms, a semi-mechanistic analysis is conducted about the change of the contact number and the obtained results are combined with observations of DEM simulations of virtual experiments to give the constitutive equations. These proposed equations are shown to capture the features observed in DEM simulations and to predict the development of contact structure correctly. Three constitutive equations regarding the contact forces are also summarised from observations. They are also shown to capture the observed features and to predict the evolution of \(\Lambda \), \({{\mathbf {\mathsf{{F}}}}}\) and \({{\mathbf {\mathsf{{F}}}}}^t\) appropriately. The proposed constitutive model not only predicts the stress, strain and void ratio of the material as traditional constitutive models do, but also predicts the statistical feature of the contact structure and contact forces. Notes Acknowledgements Open access funding provided by University of Natural Resources and Life Sciences Vienna (BOKU). The first author acknowledges the financial support from the Otto Pregl Foundation for Geotechnical Fundamental Research, Vienna, Austria and Jiangsu Province natural sciences fund subsidisation Project (BK20170677). The research is also partly funded by the Project “GEORAMP” within the RISE programme of Horizon 2020 under Grant No. 645665. Compliance with ethical standards Conflict of interest The authors declare no conflict of interest. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References 1. Dafalias, Y.F.: Must critical state theory be revisited to include fabric effects? Acta Geotech. 11(3), 479–491 (2016).  https://doi.org/10.1007/s11440-016-0441-0 CrossRefGoogle Scholar 2. Dafalias, Y.F., Li, X.S.: A constitutive framework for anisotropic sand including non-proportional loading. Géotechnique 54(1), 41–55 (2004).  https://doi.org/10.1680/geot.2004.54.1.41 CrossRefGoogle Scholar 3. Fuentes, W., Triantafyllidis, T., Lizcano, A.: Hypoplastic model for sands with loading surface. Acta Geotech. 7(3), 177–192 (2012).  https://doi.org/10.1007/s11440-012-0161-z CrossRefGoogle Scholar 4. Guo, N., Zhao, J.: The signature of shear-induced anisotropy in granular media. Comput. Geotech. 47, 1–15 (2013).  https://doi.org/10.1016/j.compgeo.2012.07.002 ADSMathSciNetCrossRefGoogle Scholar 5. He, X., Cai, G., Zhao, C., Sheng, D.: On the stress–force–fabric equation in triaxial compressions: some insights into the triaxial strength. Comput. Geotech. 85, 71–83 (2017).  https://doi.org/10.1016/j.compgeo.2016.12.011 CrossRefGoogle Scholar 6. He, X., Liang, D., Bolton, M.D.: Run-out of cut-slope landslides: mesh-free simulations. Géotechnique 68(1), 50–63 (2018a).  https://doi.org/10.1680/jgeot.16.P.221 CrossRefGoogle Scholar 7. He, X., Liang, D., Wu, W., Cai, G., Zhao, C., Wang, S.: Study of the interaction between dry granular flows and rigid barriers with an SPH model. Int. J. Numer. Anal. Methods Geomech. 42(11), 1217–1234 (2018).  https://doi.org/10.1002/nag.2782 CrossRefGoogle Scholar 8. He, X., Wu, W., Wang, S.: A constitutive model for granular materials with evolving contact structure and contact forces. Part 1: Framework. Granular. Matter. 21, 16 (2019).  https://doi.org/10.1007/s10035-019-0868-8 CrossRefGoogle Scholar 9. Kruyt, N.P., Rothenburg, L.: A micromechanical study of dilatancy of granular materials. J. Mech. Phys. Solids 95, 411–427 (2016).  https://doi.org/10.1016/j.jmps.2016.01.019 ADSMathSciNetCrossRefGoogle Scholar 10. Kuhn, M.R., Sun, W.C., Wang, Q.: Stress-induced anisotropy in granular materials: fabric, stiffness, and permeability. Acta Geotech. 10(4), 399–419 (2015).  https://doi.org/10.1007/s11440-015-0397-5 CrossRefGoogle Scholar 11. Kumar, N., Luding, S.: Memory of jamming–multiscale models for soft and granular matter. Granul. Matter 18(3), 1–21 (2016).  https://doi.org/10.1007/s10035-016-0624-2 ADSCrossRefGoogle Scholar 12. Li, X.S., Dafalias, Y.F.: Anisotropic critical state theory: role of fabric. J. Eng. Mech. 138(3), 263–275 (2012).  https://doi.org/10.1061/(ASCE)EM.1943-7889.0000324 CrossRefGoogle Scholar 13. Magoariec, H., Danescu, A., Cambou, B.: Nonlocal orientational distribution of contact forces in granular samples containing elongated particles. Acta Geotech. 3(1), 49–60 (2008).  https://doi.org/10.1007/s11440-007-0050-z CrossRefGoogle Scholar 14. Oda, M.: Initial fabrics and their relations to mechanical properties of granular materials. Soils Found. 12(1), 17–36 (1972).  https://doi.org/10.3208/sandf1960.12.17 CrossRefGoogle Scholar 15. Rothenburg, L., Bathurst, R.J.: Analytical study of induced anisotropy in idealized granular material. Géotechnique 40(4), 665–667 (1990).  https://doi.org/10.1680/geot.1990.40.4.665 CrossRefGoogle Scholar 16. Rothenburg, L., Kruyt, N.P.: Critical state and evolution of coordination number in simulated granular materials. Int. J. Solids Struct. 41(21), 5763–5774 (2004).  https://doi.org/10.1016/j.ijsolstr.2004.06.001 CrossRefzbMATHGoogle Scholar 17. Schofield, A., Wroth, P.: Critical state soil mechanics (1968)Google Scholar 18. Shi, J., Guo, P.: Fabric evolution of granular materials along imposed stress paths. Acta Geotech. 13(6), 1341–1354 (2018).  https://doi.org/10.1007/s11440-018-0665-2 CrossRefGoogle Scholar 19. Song, C., Wang, P., Makse, H.A.: A phase diagram for jammed matter. Nature 453(7195), 629–632 (2008).  https://doi.org/10.1038/nature06981 ADSCrossRefGoogle Scholar 20. Sun, J., Sundaresan, S.: A constitutive model with microstructure evolution for flow of rate-independent granular materials. J. Fluid Mech. 682, 590–616 (2011).  https://doi.org/10.1017/jfm.2011.251 ADSMathSciNetCrossRefzbMATHGoogle Scholar 21. Vairaktaris, E., Theocharis, A.I., Dafalias, Y.F.: Correlation of fabric tensors for granular materials using 2D DEM. Acta Geotech. (2018).  https://doi.org/10.1007/s11440-018-0755-1 CrossRefGoogle Scholar 22. Verdugo, R., Ishihara, K.: The steady state of sandy soils. Soils Found. 36(2), 81–91 (1996).  https://doi.org/10.1248/cpb.37.3229 CrossRefGoogle Scholar 23. Wang, R., Fu, P., Zhang, J.M., Dafalias, Y.F.: DEM study of fabric features governing undrained post-liquefaction shear deformation of sand. Acta Geotech. 11(6), 1321–1337 (2016).  https://doi.org/10.1007/s11440-016-0499-8 CrossRefGoogle Scholar 24. Wu, W., Lin, J., Wang, X.: A basic hypoplastic constitutive model for sand. Acta Geotech. (2017).  https://doi.org/10.1007/s11440-017-0550-4 CrossRefGoogle Scholar 25. Yang, L.T., Li, X., Yu, H.S., Wanatowski, D.: A laboratory study of anisotropic geomaterials incorporating recent micromechanical understanding. Acta Geotech. 11(5), 1111–1129 (2016).  https://doi.org/10.1007/s11440-015-0423-7 CrossRefGoogle Scholar 26. Zhao, J., Guo, N.: Unique critical state characteristics in granular media considering fabric anisotropy. Géotechnique 63(8), 695–704 (2013).  https://doi.org/10.1680/geot.12.P.040 CrossRefGoogle Scholar Copyright information © The Author(s) 2019 Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Authors and Affiliations Xuzhen He1View author's OrcID profileWei Wu1Email authorShun Wang11.Institut für GeotechnikUniversität für BodenkulturViennaAustria


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10035-019-0869-7.pdf

Xuzhen He, Wei Wu, Shun Wang. A constitutive model for granular materials with evolving contact structure and contact forces—part II: constitutive equations, Granular Matter, 2019, 20, DOI: 10.1007/s10035-019-0869-7