A unified constitutive model for pressure sensitive shear flow transitions in moderate dense granular materials

Granular shear flows exhibit complex transitional regimes that are dramatically affected by the pressure level and shear stress state. New advances in granular shear tests at low pressure have enlightened the understanding of the two granular shear flow transitions: between quasi-static and moderate shear flows, and between steady-state and transient shear flows. However, a unified constitutive model to describe these two transitions is yet to develop. In this work, a simplified and unified model is proposed based on innovative triaxial shear flow tests, using two dimensionless physical variables. Model results validated against experimental data suggest that the shear flow transition between a quasi-static to a moderate Isotach type flow state is highly pressure-dependent. At extremely low pressure, the granular viscosity becomes the primary mechanism, suppressing the quasi-static mechanism even under “quasi-static” shear rates. In transient to steady state granular flow transitions, a mobilized shear stress ratio or mobilized friction coefficient between zero and the critical state ratio for consolidated granular packings is taken into consideration. This is coupled with the mechanism of granular viscosity. These findings have not been discussed before and are of great relevance to granular mechanics as well as space and earthquake engineering.

Constitutive models for granular shear flows with practical engineering relevance have received increasing attention in scientific and engineering fields in recent years [1][2][3][4] . Landslides and subsidence caused by earthquakes or other dynamic impacts have caused tremendous damage. Such granular flow transitions threaten vital marine infrastructure such as submarine pipelines and cables, foundations for various offshore structures, and submerged tunnels. Current engineering practices rely heavily on engineering experience, which is in many cases insufficient. There is also no accepted method to estimate the dynamic forces of granular flows on buried structures during and post-earthquake 5 . Seismic resilient design considering large post-liquefaction deformations has received much attention. It has required increasingly more fundamental research on granular flow transitions or phase transitions under low pressure conditions to accurately describe and predict post-liquefaction behavior.
In space exploration, the rate-dependent behavior of granular materials in microgravity or low confining pressure conditions is important. With accurate constitutive models in microgravity, impact forces from landing gear can be predicted accurately and the design of such landing gear can be optimized. Liquefaction caused by the pressure increase of trapped gas in lunar regolith has also been observed on the Moon 6,7 . These observations have made research on liquefaction properties of lunar regolith or simulants in microgravity important. In particular, the large peak friction angle of sand under microgravity conditions observed in NASA-MGM experiments 8 need further understanding from the perspective of granular physics in low gravity or pressure conditions.
Scientists within the field of geotechnical engineering have proposed some engineering models for the quasistatic response of sand in pre-critical states. The critical state is a material state in which a steady plastic shear flow is assumed. Among these models, elastoplastic theory 9 and empirical models 10 based on extensive experimental data has been found to be representative. However, constitutive models that represent transitional behavior between quasi-static and moderate shear flows, and between steady-state and transient shear flows are yet to develop.
For Isotach clays subjected to an instantaneous increase in strain rates, the stress path will increase instantaneously and follow a new and higher stress path. Such an Isotach behavior of clay and its underlying mechanism of www.nature.com/scientificreports/ pore water viscosity are mostly absent in sand and other granular materials 11 . Recent findings suggest that more complicated viscosity behavior of granular materials such as TESRA and P&N behavior, are present 12 . It is not uniformly agreed upon which fundamental mechanisms such as viscosity, friction, or others, are correct. The occurrence of transitions between TESRA and Isotach behavior are also of question. A unified constitutive model that can describe the granular flow transitional behavior from a quasi-static state to a moderate flow state is urgently needed. The constitutive models representing quasi-static stress-strain relationships in sands are mainly divided into elastoplastic and nonlinear "elastic" models. An elastoplastic model is the critical Camclay model 13 , which uses a specific yield surface, plastic flow criterion, and hardening criterion. Nonlinear "elastic" models are often empirical because of their simple form, which include Duncan-Zhang 14 , K-G 15 , and others. Theoretical Cauchy and Green nonlinear elastic incremental models are also available, though used less frequently 16 .
In general, nonlinear empirical models represent the quasi-static behavior in sand or other granular materials through the relationship between mobilized shear stress ratio or mobilized friction coefficient and strain. For very small shear strains of less than 10 −5 , granular materials are generally in a pure elastic state and the nonlinear shear modulus is constant. In the elastic state, the small-strain shear modulus is affected by volume fraction, effective stress, particle size distribution, and particle shape. When the shear strain increases, the nonlinear shear modulus degrades continuously until it becomes almost zero for very large shear strains. In this state, the granular material has reached the critical state where it can be sheared continuously with a constant shear stress. Oztoprak and Bolton 17 proposed an empirical model for the nonlinear shear modulus by evaluating extensive experimental data of sand stiffness.
The research group MiDi (Groupement de Recherche Milieux Divisés) 1,2 established the µ(I) constitutive rheological model for steady dense granular flows based on extensive experiments and numerical simulations. This model can describe steady flow behavior of granular materials for moderate strain-rates by using the inertial number I , a dimensionless relative strain rate, as a variable. This model was used as the plastic flow rule in order to formulate a nonlinear elastic and perfectly viscoplastic model for dense granular flows 3 . However, it fails to represent transient to steady flow transitions of granular materials.
In space and earthquake engineering problems, sands or other granular materials are often in a quasi-static mobilized stress state far below the critical quasi-static stress state when subject to static loads with a load around one third to one fourth of the critical load. When granular materials are subject to seismic or dynamic loads such as earthquakes and landslides, transient flows are more likely to occur than steady flows. Such transient flows include creep or decreased flow and creep rupture or accelerated flow. Neither the MiDi model nor quasi-static soil mechanics models are suitable in this case.
A series of laboratory tests on sands at low effective isotropic stress were systematically performed to reveal transient flows 8,[18][19][20] . Results from laboratory tests reveal a significant S-shape between strain rates and viscous shear stresses. This behavior is not captured by the rate-dependent Bingham model nor by rate-independent critical state soil mechanics models. Moreover, the underlying mechanism is not yet known, and the pressure sensitive characteristics of such transient flows is yet to be investigated 21,22 .
In this paper, a novel constitutive model is proposed based on advanced triaxial shear flow tests under the conditions of isotropy and moderate density of the sample, using two dimensionless physical variables: the mobilized friction coefficient and the inertial number. The mobilized friction coefficient is used to capture the transient to steady flow transition of granular materials starting from any quasi-static noncritical state. Its development can be formulated in terms of the shear strain variation. The inertial number is adopted to characterize the rate and pressure dependency of granular flows. The model is validated against experimental data and the results suggest that the shear flow transition between a quasi-static to a moderate Isotach type flow state is highly pressure-dependent. At extremely low pressure, the granular viscosity becomes the primary mechanism, suppressing the quasi-static mechanism even under "quasi-static" shear rates. In transient to steady state granular flow transition, the mobilized shear stress ratio between zero and the critical state ratio for any consolidated granular packings should be taken into consideration which interacts with the mechanism of granular viscosity. Predictions and analysis based on the proposed model are done to provide an understanding of the effects of low effective isotropic stresses and shear rates on granular materials with an initial moderate density. Pressure sensitive and rate-dependent transient flow characteristics of granular materials are disclosed for the first time.

Model development
Conceptual framework. Figure 1 shows the relationship between the shear strain γ and the effective friction coefficient µ and its components. This shear stress-strain plane with shear flow rate effects is commonly used to characterize the shear flow behavior of sand and other granular materials. The upper flat line µ 1 signifies the threshold in the granular material between a steady dense moderate flow for µ < µ 1 and an accelerated unsteady dense flow for µ > µ 1 . In the latter case, which is beyond the scope of this paper, the long-term contact between particles is no longer maintained and binary collision becomes the main mechanism of momentum exchange. The lower curved line OC is the quasi-static mobilized friction coefficient µ m . The mobilized friction coefficient µ m represents the stress-strain relation in a quasi-static state and its upper bound µ s is the critical state of Drucker-Prager yield surface. It signifies the threshold in the granular material between a quasi-static state for µ < µ m and a moderate shear flow for µ > µ m .
In a steady state when µ > µ s and γ > γ s , the viscous stress in a steady dense moderate flow is believed to be a function of flow rate γ and effective isotropic stress p ′ . For most of the cases in seismic and dynamic engineering, the soils are statically loaded in the transient domain at point A and not in the steady domain at point C in Fig. 1. Thus, the granular materials cannot be characterized by a steady shear flow and will mainly inherit properties of a transient shear flow. www.nature.com/scientificreports/ The concept for granular materials is analogous with the concept of phase transitions for matters such as water, but the scale in which the particles interact are different. In water, viscosity is determined by intermolecular forces such as van der Waals forces and hydrogen bonds. In granular matters, it is thought that the main intergranular mechanism is friction. The inertial number I is used to characterize the granular friction-based viscosity, which is related to the change of intrinsic time scale or frequency of the granular flow due to different effective isotropic stress levels.
In this paper, three friction coefficients have been introduced, namely the critical friction coefficient µ s , the upper limit friction coefficient µ 1 and the mobilized friction coefficient µ m . Of the three mentioned friction coefficients, the upper limit friction coefficient µ 1 is not well known in granular mechanics. In fluid mechanics and granular physics, the mobilized friction coefficient µ m is not well known. The mobilized friction coefficient µ m is the quasi-static friction coefficient of a granular material in any pre-consolidated state, or for any effective isotropic stress p ′ . When a granular material with zero deviatoric stress τ is quasi-statically sheared from a pre-consolidated isotropic state to a critical state, the mobilized friction coefficient µ m = µ m p ′ , γ gradually develops from zero to the critical friction coefficient µ s .

Mathematical formulation. The original MiDi model in tensor form assumes granular flows incompress-
ible and is expressed as Eqs. (1) and (2): with µ(I) and I defined in Eqs. (3) and (4): www.nature.com/scientificreports/ In the above equations, µ(I) is the rate-dependent and pressure dependent friction coefficient, p ′ is the effective isotropic stress or the effective hydrostatic stress, γ ij is the deviatoric strain rate tensor and γ = γ ijγij /2 is the second invariant of γ ij . In Eqs. (3) and (4), µ s = tan θ s is the critical friction coefficient in a quasi-static state, µ 1 = tan θ 1 is the upper limit friction coefficient, I 0 is a material parameter reflecting the rheological properties of the granular material, I is the inertial number and d is the mean particle diameter. The inertial number can be expressed as the ratio between the microscopic time scale T p = d/ p ′ /ρ s and the macroscopic time scale T γ =γ −1 . The microscopic time scale T p represents the time required for a particle to undergo a displacement of the mean particle diameter under an effective isotropic stress p ′ . The macroscopic time scale T γ represents the time required to generate a unit of strain.
In this paper, an extension and modification of the MiDi model is formulated by separating the deviatoric stress τ ij into a rate-independent quasi-static component τ s ij and a rate-dependent viscous component τ v ij as shown in Eq. (5): In the quasi-static component τ s ij , the critical friction coefficient µ s is replaced by the mobilized friction coefficient µ m with an upper bound limit of µ m ≤ µ s as shown in Eq. (6): In the original MiDi model, the direction of the critical shear stress is associated with the direction of the deviatoric strain rate γ ij . This is corrected in the present model as the static shear stress τ s ij is a function of the nonlinear shear modulus G and deviatoric strain tensor γ ij as shown in Eq. (7): where G is the shear modulus. τ v ij is extended to include transient flows as shown in Eq. (8).
We propose to calculate the quasi-static deviatoric stress by using the nonlinear expression for the shear modulus G proposed by Oztoprak and Bolton 17 . The nonlinear shear modulus can be expressed as Eq. (9): where p a is the standard atmospheric pressure, e is the void ratio, p ′ is the effective isotropic stress and A(γ ) and m(γ ) are fitting parameters that are dependent on the shear strain γ.
We also propose that A(γ ) and m(γ ) can be expressed as Eqs. (10) and (11): where m 0 , m 1 and A 0 are parameters directly determined from experimental results, γ r1 and γ r2 are reference strain levels and M = µ s / √ 3 is the critical stress ratio in a triaxial compression test. γ r1 can be obtained from the curve fitting of m(γ ) , and γ r2 can be determined from Eq. (12). From the extreme cases of G | γ →0 = G 0 and G | γ →γ s = 0 , Eq. (9) can be calibrated to successfully predict the nonlinear behavior of the shear modulus G of the sand sample with an initial moderate density and subjected to unrotated principal stress conditions. Since the experimental database established by Oztoprak and Bolton does not include data for effective isotropic stresses of p ′ < 10 kPa , the applicability of the model for such effective isotropic stresses must be verified. In this paper, it is assumed that the model is still valid for low effective isotropic stresses according to 23 .
In Fig. 2, observe that Oztoprak and Bolton's expression for the nonlinear shear modulus correlates well with experimental data. Thus, it can be deemed to represent the behavior accurately.
Model characteristics. The proposed model described in "Mathematical formulation" section, is a simple unified constitutive model for incompressible granular flows starting from any quasi-static stress state, includ-  (14) and (15): where µ v is the viscous friction coefficient. Together, the quasi-static component q s and the viscous component q v can describe transient granular flows under the conditions of isotropy and moderate density of the sample starting from any quasi-static stress state but the critical state.

Simulations of transient granular flows
Multi-stage triaxial creep tests in low effective isotropic stresses. To investigate post-liquefaction flow mechanisms of sand, Towhata et al. 19 performed triaxial experiments on sand for low effective isotropic stresses p ′ . A specific multi-stage creep stress path was used to obtain the correlation between the transient strain rate γ and the driving viscous stress q v . In every loading stage, the discrepancy between the loading state and the quasi-static or reference curve is the viscos stress component q v . As creep develops at each load stage, the quasistatic component q s increases with the development of shear strain γ , and the viscous component q v decreases as the shear strain rate γ decreases.
Although the triaxial setting has the advantage of investigating granular transient flow studies under controlled effective isotropic pressure, it cannot introduce the stress and strain-rate tensor rotation. However, a triaxial setting allows us to achieve an insight into the pressure sensitive granular flows. In this paper, numerical simulations of the two kinds of triaxial experiments are performed. General material parameters from Towhata et al. 19,20 are used and model specific parameters are derived from experimental data. The critical stress ratio M is obtained in the critical state from the data set obtained by Been et al. 24 , A 0 is determined by measuring the smallstrain shear modulus G 0 in Eq. (9) for γ → 0 , and γ r1 and γ r2 are calibrated to fit the quasi-static stress-strain relationship. The parameters used in the numerical simulations are shown in Table 1, where M and µ 1 are material constants, and I 0 is a state related parameter. Figure 3 show experimental and simulation results of monotonic multi-stage triaxial creep flow on Yurakucho sand. In Fig. 3a, the red dotted line is the measured deviatoric stress, the red solid line is the quasi-static reference curve of experiment data, the black solid line is the simulated deviatoric stress, and the green solid line is the simulated reference curve. In Fig. 3b, the red dotted line is the experimental line and the black solid line is www.nature.com/scientificreports/ the simulated time-deviatoric strain curve. In the experiment, the reference curve represents the quasi-static stress-strain relationship defined by the stress-strain points at the end of each loading step. In Fig. 3a, large deviatoric step loads q are applied. We assume that the experiment sample reaches a quasistatic state when the deviatoric strain rate γ is less than 10 −8 s −1 . For constant deviatoric strains γ , the sudden application of each deviatoric step load q will result in a sudden increase of q due to the sudden increase of the effective isotropic stress p ′ .These effects cause the jump and following derivation in the simulated reference curve according to Eq. (14) compared to experimental reference curve that is simply smooth connection of ending points of each load step. The experiments did not accurately maintain a constant stress which is the source of the discrepancy between the experimental and simulated deviatoric stress curve.
In Fig. 3b, the deviatoric strain γ develops rapidly after the application of every step load q . The deviatoric strain γ appears to reach a steady state at each loading step after sufficient time has passed. The development of the deviatoric strain γ may be explained through the interplay between the quasi-static component q s and the viscous component q v . When a deviatoric stress step q is applied to any quasi-static state q s . A transient viscous component q v = q − q s is generated. Due to the nonzero viscous component q v = 0 , deviatoric strain γ develops. Because of the increase in deviatoric strain γ , the quasi-static component q s in the sample increases gradually towards its critical value q s c and the viscous component q v decreases in the same manner. This leads to the increase of the deviatoric strain γ with a decreasing deviatoric strain rate γ . Eventually, the sample reaches a quasi-static state.
The creep time and strain in each loading step increases stepwise. This can largely be explained by the shear modulus G . When the deviatoric stress q is small, the shear modulus G is large and the strain increment �γ is small. When the deviatoric stress τ is higher, the deviatoric strain increment �γ becomes larger for similar  www.nature.com/scientificreports/ deviatoric stress increments q . Increasing deviatoric strain �γ for increasing values of the deviatoric stress q means that more creep time is needed for the sample to reach a new quasi-static state. The simulation and experimental results shown in this section reflect the correlation between the viscous component q v and the deviatoric strain rate γ in the sample, and that rate-dependency is affected by the effective isotropic stress p ′ . In the following sections, We will present predictions and analyze the rate-dependent properties of sandy materials at different effective isotropic stress states p ′ and deviatoric strain rates γ . These properties are represented through the ratio q v /q s c and the deviatoric strain rate γ. The main factors influencing the viscous component q v are the effective isotropic stress p ′ , deviatoric strain rate γ and volume fraction . In this section, the effects of these factors are studied. Based on Eqs. (13)-(15), the rate-dependent behavior of Yurakucho sand for different effective isotropic stresses p ′ are simulated with the parameters calibrated in Table 1. q v /q s c for different effective isotropic stresses p ′ and deviatoric strain rates γ are then obtained.
In general soil mechanics experiments, the effective isotropic stress p ′ is in the range between 50 kPa and several thousand kPa , and the deviatoric strain rate γ is less than 10 −5 s −1 . As shown in Fig. 4a, when the deviatoric strain rate γ is less than 10 −5 s −1 of a quasi-static strain rate commonly adopted in soil mechanics lab, q v /q s c approaches zero. It implies that the viscous component q v is negligible compared to the critical quasi-static component q s c , and is consistent with the perception that sand is rate-independent in general soil mechanics experiments. The corresponding inertial number in general soil mechanics experiments can be calculated as around 10 -8 based on the Eq. (4) and Table 1. When the deviatoric strain rate γ is higher than 10 −4 s −1 or the effective stress level is decreased in orders, the viscous component q v is no longer negligible as the inertial number increases in a few orders. The effects of both effective isotropic stress and strain rate can be alternately characterized by the dimensionless Inertial number in the range between 10 −8 and 10 −4 shown in Fig. 4b. If the effective isotropic stress p ′ is reduced to less than 0.5 kPa, the range in which the viscous component q v become significant is shifted towards smaller deviatoric strain rates γ . The increasing significance of the viscous component q v for low effective isotropic stresses p ′ implies that the sand is no longer rate-independent in such conditions. Earthquake induced liquefaction typically cause deviatoric strain rates of about γ ≈ 10 −3 s −1 . Thus, the viscous component becomes non-negligible and may severely affect underground structures such as tunnels, pipelines, pile foundations and more. Experimental data from triaxial creep tests at low effective isotropic stresses p ′ obtained by 19 are reproduced in Fig. 4b. The experimental data fit the model predictions in general tendency fairly well.

Rate-controlled triaxial shear flow tests.
In the previous sections, creep properties of sands were observed in stress-controlled experiments and simulations. In this section, rate-dependent properties of sands are studied by the proposed model using experimental data by Enomoto et al. 25 .
Sand exhibits different viscous behavior when subject to different deviatoric strain rates γ . In experiments with varying deviatoric strain rates γ , different types of sands will have unique viscous reactions which can be seen in stress-strain curves. The rate-dependent behavior of sand can be classified into viscosity types such as Isotach, TESRA (Temporary Effects of Strain Rate and Acceleration) and P&N (Positive and Negative). The three mentioned types of granular viscosity can be seen in Fig. 5.
The most common granular viscosity type is Isotach viscosity. If Isotach sands are subjected to an instantaneous increase in strain rates, the stress path will also increase instantaneously and follow a new and higher stress path. For TESRA viscosity, the stress path does not exhibit the same uniquely defined stress path as Isotach viscosity. If TESRA sands are subjected to an instantaneous increase in strain rate, the stress path will only temporarily follow a higher stress path before it reverts to the original stress path independent of the new strain rate. For P&N viscosity, the stress path may fall to a lower stress path after an instantaneous increase in strain rate. However, P&N viscosity is rare and has not been explicitly confirmed by experiments.  Table 1. Figure 6 shows both the simulated relationship between the deviatoric stress q and the axial strain ε a , and experimental results from Enomoto et al. for a strain-rate controlled triaxial shear experiment. In both figures  Figure 6a shows the experimental results and simulations for Miho sand with an effective confining pressure of 50 kPa. In the experiment, the axial strain rate ε a is varied in the interval 0.1ε 0 ≤ε a ≤ 10ε 0 where the reference strain rate ε 0 is defined as ε 0 = 0.0625 % min −1 = 1.04 × 10 −5 s −1 . When the axial strain rate ε a is 0.1ε 0 , the experimental results and the quasi-static reference curve coincides, indicating that the rate-dependent viscous component q v is small. When the axial strain rate ε a changes suddenly with a magnitude �ε a , the deviatoric stress q changes accordingly with a magnitude q . The change in deviatoric stress q shows a strong positive correlation with the change in axial strain rate �ε a which is typical Isotach behavior as described in Fig. 5 as any positive change of the axial strain rate ε a will significantly change the viscous component q v defined in Eq. (15).
When the deviatoric stress q reaches 174.0 kPa, the deviatoric stress q is kept constant and the sand is creeped for one hour. The axial strain ε a gradually approaches the reference curve during the one-hour creep. After the one-hour creep, loading is resumed with an axial strain rate ε a =ε 0 . The stress path makes an immediate jump and follows a new and higher stress path.
It is shown that the model simulates Isotach behavior of Miho sand well. The rate-dependency of Miho sand in the simulation is consistent with experimental results at different strains. The change in deviatoric stress q is strongly correlated with the change of axial strain rate �ε a , and the magnitude of the change in deviatoric stress q does not decay with continued axial strain ε a . Figure 6b shows the experimental results and simulations for Toyoura sand with an effective confining pressure of 400 kPa. In the experiment, the axial strain rate ε a is varied in the interval 0.1ε 0 ≤ε a ≤ 20ε 0 where the reference strain rate ε 0 is defined as ε 0 = 0.005% min −1 = 8.33 × 10 −7 s −1 . But from Fig. 6b, the experimental and simulated results show that Toyoura sand exhibit weak Isotach behavior. The discrepancy between the experimental and simulated curves are small, and there is no distinct rate-dependent viscous component q v .The absence of a distinct rate-dependent viscous component q v suggests that Toyoura sand with an effective confining pressure of 400 kPa is TESRA viscous. However, the sudden increase in the deviatoric stress q and corresponding increase in axial strain rate �ε a is not well simulated. This is to be studied further. The limitation of this model does not affect its overall rate-independency for TESRA viscous materials.

Pressure dependent viscosity of granular materials. In the two previously presented experiments,
Toyoura sand showed different transient flow behavior. In the triaxial creep experiment with low effective isotropic stress by Towhata et al., significant granular viscosity was demonstrated. In Enomoto et al. 's experiment, Toyoura sand exhibited much less viscosity behavior with TESRA viscosity when subject to an effective isotropic stress of 400 kPa. It is generally believed that Isotach and TESRA behavior of sands are determined mainly by intrinsic material properties. However, no clear physical mechanism and quantitative description are explicitly agreed upon. In the model presented in this paper, different viscous behavior of granular materials is highly sensitive to the effective stress level shown in Fig. 4. The model simulations uncover the pressure dependent physical mechanism of viscosity of granular materials. Based on proposed model in this paper, simulations in Fig. 7 show that when the effective confining pressure of Miho sand is raised from 50 to 500 kPa, the flow behavior will make transition from strong Isotach to weak Isotach or TESRA for both varying triaxial strain rates and creep tests. In Fig. 7a, the strain rate is suddenly increased to 100 times faster, in Fig. 7b, the deviatoric stress is suddenly increased and held on to measure the creeping. Similarly, when the effective stress of Toyoura sand is reduced from 400 to 40 kPa, the rate-dependent behavior will make transition from weak Isotach or TESRA to strong Isotach.  www.nature.com/scientificreports/ Observe that the flow behavior of sand cannot be singly attributed to material properties, but that it is also the result of granular interaction properties which are dependent on external parameters such as the effective isotropic stress. In the description of this model, the inertial number I reflects the influence of the effective isotropic stress on the viscosity behavior of granular materials. Sands under very low effective isotropic stresses may respond in a pressure sensitive viscous and not quasi-static manner, even if the loading rate is very low. The described response is due to the low natural frequency and increased viscosity in such conditions. This behavior is believed to be responsible for the large peak friction angles observed by Sture et al. in the NASA-MGM program.

Conclusions
In this paper, the proposed model has demonstrated a good ability to accurately simulate pressure sensitive shear flow behavior of granular materials under the conditions of isotropy and moderate density of the sample. This has been achieved by separating the total stress τ ij into a viscous part τ v ij described by the viscous friction coefficient µ v , and a quasi-static component τ s ij described by the mobilized friction coefficient µ m . Both the viscous friction coefficient µ v and the mobilized friction coefficient µ m are quantified by extensive experimental data.
We propose that the pressure sensitivity of the granular viscous component τ v exists in all transient and steady flow states. Based on the quasi-static component τ s , the viscous component τ v must always be considered to determine whether the viscous mechanism is of significance or not. In our proposed model, the viscous component τ v is generally neglected in quasi-static loading when the deviatoric strain rate γ is equal to or smaller than 10 −4 s −1 for the effective isotropic stresses used in laboratories.
In traditional studies, the flow behavior of granular materials is mainly attributed to material properties. However, lower effective isotropic stresses p ′ may increase the effective friction coefficient significantly which can result in a non-negligible viscous component τ v . The granular viscosity of Isotache behavior in low effective isotropic stresses can occur for granular materials with the presence of a non-negligible viscous stress component τ v . This should be accounted for when designing buried infrastructure against earthquake induced liquefaction. The pressure sensitive granular viscosity found in this paper has implications on the understanding of large peak friction angles found in NASA-MGM experiments for extremely low effective isotropic stresses.
In the field of granular physics, the MiDi rheological model describes a steady dense granular flow. The MiDi rheological model cannot be directly applied to problems concerning transient flows such as creep and strain rate-controlled shear flow and earthquake liquefaction. Simulations and experiments on sands in this paper for low effective isotropic stresses and varying strain rates, have revealed that the sands and other granular materials also exhibit Isotach behavior pre-critical state. The sensitivity of the experiments to low effective isotropic stresses are in general neglected in the experiments. Understanding and modelling transient flows of sand starting from pre-critical states will have a profound influence in engineering.
For problems concerning transient behavior of sands such as earthquakes, landslides, dynamic actions on subsea structures, granular phase-transition problems and problems with high deviatoric strain rates γ or low effective isotropic stresses p ′ , the proposed model is more accurate than current accepted constitutive models, most of which are either rate independent or in the steady state. For the experimental data in Fig. 4, the total stress will be underestimated up to 20% using rate independent model. The original MIDI model is not capable of modeling the Fig. 4 data. Meanwhile the proposed model highlights the role of effective stress in the granular shear flow transitions under the conditions of isotropy and moderate density of the sample. We believe the uncovered new physical mechanism in our proposed model should be valid for sand as well as other granular systems. www.nature.com/scientificreports/ However, quantitative verification of the proposed model for different scenarios and different granular systems must be carried out to clarify its limitations and relevance to different problems in constitutive or representative element volume and boundary-valued configurations.