Phenomenological fractional stress–dilatancy model for granular soil and soil-structure interface under monotonic and cyclic loads

Abstract

The stress–dilatancy relation is of critical importance for constitutive modelling of geomaterial. A novel fractional-order stress–dilatancy equation had been developed for granular soil, where a nonlinear stress–dilatancy response was always predicted. However, it was experimentally observed that after a certain extent of shearing, an almost linear response between the stress ratio and the dilatancy ratio, rather than the nonlinear response, usually existed. To capture such stress–dilatancy behaviour, a new fractional stress–dilatancy model is developed in this study, where an apparent linear response of the stress–dilatancy behaviour of soil after sufficient shearing is obtained via analytical solution. As the fractional order varies, the derived stress–dilatancy curve and the associated phase transformation state stress ratio keep changing. But, unlike existing researches, no other specific parameters, except the parameter related to fractional order, concerning such shift are required. Then, the developed stress–dilatancy model is applied to constitutive modelling of granular soil and soil–structure interface, for further validation. A series of test results of different granular soils and soil–structure interfaces under different loading conditions are simulated and compared, where a good model performance is observed.

Appendix

1.1 Constitutive modelling of granular soil

To demonstrate the performance of the developed stress–dilatancy equation, the elegant work by Dafalias and Manzari [2] is employed and modified here, by incorporating the fractional-order stress–dilatancy equation. Accordingly, the general constitutive relation for granular soils under triaxial loading can be provided as:

$$ {\dot{\mathbf{\varepsilon }}} = \left[ {{\mathbf{C}}^{e} { + }\frac{1}{H}{\mathbf{n}}^{{\varvec{T}}} {\mathbf{m}}} \right]{\dot{\mathbf{\sigma }}}{^{\prime}} $$

(27)

where the incremental strain vector ${\dot{\mathbf{\varepsilon }}} = \left[ {\dot{\varepsilon }_{v} ,\dot{\varepsilon }_{s} } \right]^{T}$, in which the superimposed dot (∙) indicates increment; the incremental stress vector ${\dot{\mathbf{\sigma }}}{^{\prime}} = \left[ {\dot{p}{^{\prime}} ,\dot{q}} \right]^{T}$; ${\mathbf{C}}^{e}$ is the elastic compliance vector, which can be expressed as:

$$ {\mathbf{C}}^{e} = \left[ {\begin{array}{*{20}c} {1/K} & 0 \\ 0 & {1/\left( {3G} \right)} \\ \end{array} } \right] $$

(28)

where the bulk modulus and shear modulus can be, respectively, defined as [5]:

$$ K = \frac{{2G\left( {1 + \nu } \right)}}{{3\left( {1 - 2\nu } \right)}} $$

(29)

$$ G = G_{0} \frac{{\left( {2.97 - e} \right)^{2} }}{1 + e}\sqrt {p_{a} p{^{\prime}} } $$

(30)

where $\nu$ is the Poisson’s ratio, defining the lateral deformation ability of the material, which can be determined by fitting the relation between the first and third principal strains; $p_{a} = 100{\text{ kPa}}$ is the atmospheric pressure; $G_{0}$ is the elastic constant, which can be determined by fitting the initial stress–strain curve of sand. $G_{0}$ determines the elastic behaviour of the material.

Moreover, n and m in Eq. (27) are the plastic flow and loading vectors, respectively. n can be defined by using the fractional-order stress–dilatancy equation (21), such that:

$$ {\mathbf{n}} = \left[ {\left( {\mu M - t\eta } \right)\left| \eta \right|^{1 - \alpha } ,{ }t} \right]^{T} $$

(31)

The plastic loading/unloading vector m can be defined by following Dafalias and Manzari [2], such that:

$$ {\mathbf{m}} = \left[ { - \left( {t\alpha_{r} + m} \right),t} \right]^{T} $$

(32)

$\alpha_{r}$ is the back stress ratio defined by Dafalias and Manzari [2], which determines the kinematic hardening of the loading surface of granular soil ($\dot{\alpha }_{r} = \dot{\eta }$); m defines the size of the elastic region of the loading surface, which is set as zero in this study for the sake of simplicity [3]. Therefore, the monotonic loading vector at compression and extension can be derived as:

$$ {\mathbf{m}} = \left\{ {\begin{array}{*{20}c} {\left[ { - \left( {\alpha_{r} + m} \right),1} \right]^{T} ,{\text{ compressive loading}}} \\ {\left[ {\left( {\alpha_{r} - m} \right), - 1} \right]^{T} ,{\text{ extensive loading }}} \\ \end{array} } \right. $$

(33)

while for cyclic unloading from compression and extension, Eq. (32) can be derived as:

$$ {\mathbf{m}} = \left\{ {\begin{array}{*{20}c} {\left[ {\left( {\alpha_{r} - m} \right), - 1} \right]^{T} ,{\text{ unloading from compression}}} \\ {\left[ { - \left( {\alpha_{r} + m} \right),1} \right]^{T} ,{\text{ unloading from extension }}} \\ \end{array} } \right. $$

(34)

Furthermore, in Eq. (27), H is the hardening modulus. The hardening modulus describes a criterion to determine the incremental plastic strains induced by the particle rearrangement and even particle degradation, i.e. the change of phase position between particles and the amount of particle breakage, at a given stress level. The hardening function in the constitutive model represents the increment of yield surface under certain loading conditions, and is generally a function of the plastic strain, stress level, material state and even fabric [39]. For the sake of simplicity, the following formula for H is defined:

$$ H = h_{0} \left( {1 - e} \right)c_{l} G{\text{exp}}\left( {k\Psi } \right)\frac{{\left( {M_{p} - m - t\alpha_{r} } \right)}}{{\left| {\alpha_{r} - \alpha_{in} } \right|}} $$

(35)

where $h_{0}$,${ }c_{l}$ and k are the material constants, in which $h_{0}$ captures the influence of void ratio on the plastic hardening modulus, while $c_{l}$ reflects the impact of reloading on the hardening behaviour of granular soil. $h_{0}$ and $c_{l}$ can be determined by fitting the $q - \varepsilon_{1}$ relationship during virgin loading and first reloading, respectively. The role of $c_{l}$ is only activated for reloading, in other words, $c_{l} = 1$ when unloading takes place. Furthermore, $M_{p} = M\exp \left( { - k\Psi } \right)$ is the peak stress ratio, in which k can be determined by fitting the peak stress ratio points, as comprehensively discussed in Sun et al. [36]. Note that k defines the dependence of peak failure state on material state. $\alpha_{in}$ is the stress ratio at the last stress reversal; therefore, $\alpha_{in} = 0$ for the initial loading case. Therefore, in Eq. (35), the effects of plastic strain level (e), stress level ($\alpha_{r}$) and material state ($\Psi$) are all considered. However, the material fabric is not taken into account in the current study and needs further comprehensive investigation.

1.2 Constitutive modelling of soil–structure interface

To model the stress–displacement behaviour of soil–structure interface, the constitutive relation should be reformulated in terms of the incremental normal/tangential stress and the incremental normal/tangential displacement. Accordingly, the constitutive relation for soil–structure interface under monotonic and cyclic loads can be reformulated as [13, 23]:

$$ {\dot{\mathbf{\sigma }}} = \frac{1}{{t_{n} }}\left( {{\mathbf{D}}^{e} - \frac{{{\mathbf{D}}^{e} {{\tilde{\mathbf{n}}\tilde{\mathbf{m}}}}^{T} {\mathbf{D}}^{e} }}{{K_{p} + {\tilde{\mathbf{m}}}^{T} {\mathbf{D}}^{e} {\tilde{\mathbf{n}}}}}} \right){\dot{\mathbf{U}}} $$

(36)

where ${\dot{\mathbf{\sigma }}} = {[}\dot{\sigma }_{n} ,\dot{\tau }{]}^{T}$ is the incremental stress vector, in which $\sigma_{n}$ and $\tau$ are the normal and tangential stresses, respectively; ${\dot{\mathbf{U}}} = \left[ {u,v} \right]^{T}$ is the displacement vector, in which $u$ and $v$ are the normal and tangential displacements, respectively; $t_{n}$ is the thickness of the interface, which is 5–10 times of the median particle size of the interface material; following Saberi et al. [23], $t_{n} = 5d_{50}$ is used in this study. The elastic stiffness tensor ${\mathbf{D}}^{e}$ can be expressed as:

$$ {\mathbf{D}}^{e} = \left[ {\begin{array}{*{20}c} {D_{n} } & 0 \\ 0 & {D_{t} } \\ \end{array} } \right] $$

(37)

where the elastic moduli in the normal and tangential directions can be, respectively, defined as:

$$ D_{n} = D_{n0} \sqrt {\sigma_{n} /p_{a} } $$

(38)

$$ D_{t} = D_{t0} \sqrt {\sigma_{n} /p_{a} } $$

(39)

where $D_{n0}$ and $D_{t0}$ are material constants, defining the elastic behaviour of the material, which can be determined by fitting the stress–displacement curve at the initial loading stage.

Moreover, ${\tilde{\mathbf{n}}}$ and ${\tilde{\mathbf{m}}}$ in Eq. (36) describe the plastic flow and loading directions of the interface, respectively. The plastic potential function for soil–structure interface can be provided by reformulating Eq. (13) in terms of the stresses, $\tau$ and $\sigma_{n}$, where $\tilde{f} = \tau + t\tilde{\eta }_{c} \sigma_{n} {\text{ln}}\sigma_{n} - t\tilde{\eta }_{c} \sigma_{n} {\text{ln}}\sigma_{n0} = 0$. $\tilde{\eta }_{c}$ is the stress ratio obtained at the ultimate state of the interface, where the ultimate state is similar to the critical state of soil; $\sigma_{n0}$ is the intercept of $\tilde{f}$ with the $\sigma_{n} -$ axis. Due to the similar expressions of $\tilde{f}$ and $f$, ${\tilde{\mathbf{n}}}$ can be simply defined through the same algebraic operations shown in Eqs. (15)–(21). Therefore, analogous to Eq. (31), one has the following plastic flow direction for soil–structure interface:

$$ {\tilde{\mathbf{n}}} = \left[ {\left( {\mu \tilde{\eta }_{c} - t\tilde{\eta }} \right)\left| {\tilde{\eta }} \right|^{1 - \alpha } ,{ }t} \right]^{T} $$

(40)

where the stress ratio $\tilde{\eta }$ within the interface is defined as $\tilde{\eta } = \tau /\sigma_{n}$. In addition, cyclic experimental results [46, 48, 49] indicated that interface contraction took place, followed by dilation after each shear stress reversal, where the dilation rate was lower than the contraction rate, resulting in accumulated interface contraction [25, 26]. Therefore, to reasonably capture such dilatancy behaviour of the interface, the following modified expression for $\alpha$ is suggested:

$$ \alpha = {\text{exp}}\left( {\frac{{\tilde{\Delta }\Psi }}{{1 + W_{p}^{\chi } }}} \right) $$

(41)

where $\tilde{\Delta }$ and $\chi$ are the material constants, defining the dependence of plastic flow direction of the interface on material state; $\tilde{\Delta }$ can be determined from the phase transformation state of the interface at virgin loading, while $\chi$ can be obtained by fitting the unloading stress–displacement curve of the interface during cyclic loading; $W_{p}$ is the total plastic work at the last shear stress reversal; therefore, $W_{p} = 0$ for virgin loading case. For cyclic loading case, $W_{p}$ can be computed as:

$$ W_{p} = \frac{1}{{t_{n} }}\smallint \sigma_{n} du^{p} + \tau dv^{p} $$

(42)

where the superscript (p) indicates plastic displacement; $\langle \rangle$ is the McCauley bracket. For soil–structure interface, the state parameter $\Psi$ is defined as:

$$ \Psi = e - \tilde{e}_{{\Gamma }} + \tilde{\lambda }{\text{ln}}\left( {\frac{{\sigma_{n} }}{{p_{a} }}} \right) $$

(43)

where $\tilde{e}_{{\Gamma }}$ and $\tilde{\lambda }$ are the material constants, describing the ultimate state line of the interface, which can be determined by fitting the ultimate state data points in the $e - {\text{ln}}\sigma_{n}$ plane. This procedure is similar to the one used with soil; however, the void ratio (e) for the interface is defined as:

$$ de = \frac{du}{{t_{n} }}\left( {1 + e_{0} } \right) $$

(44)

where $e_{0}$ is the initial void ratio, prior to virgin shearing.

The plastic loading/unloading vector ${\tilde{\mathbf{m}}}$ can be defined by following Saberi et al. [24], such that:

$$ {\tilde{\mathbf{m}}} = \left[ { - \left( {t\tilde{\alpha }_{r} + \tilde{m}} \right),t} \right]^{T} $$

(45)

where $\tilde{\alpha }_{r}$ is the back stress ratio, defining the kinematic hardening (location) of the yielding surface ($\mathop {\tilde{\alpha }_{r} }\limits^{.} = \mathop {\tilde{\eta }}\limits^{.}$); $\tilde{m} = 0.01\tilde{\eta }_{c}$ is the size of the yielding surface, defining the elastic region of the interface. During numerical implementation, the sign function (t) is determined by the following logic:

$$ t = \left\{ {\begin{array}{*{20}c} { + 1,{ }\tilde{\eta } - \tilde{\alpha }_{r} \ge 0} \\ { - 1,{ }\tilde{\eta } - \tilde{\alpha }_{r} < 0} \\ \end{array} } \right. $$

(46)

Finally, the plastic modulus $K_{p}$ in Eq. (36) can be defined by measuring the distance between the current yielding surface and the bounding surface of the interface, such that:

$$ K_{p} = K_{0} D_{n} \frac{{\left( {\tilde{\eta }_{p} - \tilde{m} - t\tilde{\alpha }_{r} } \right)}}{{\left| {\tilde{\alpha }_{r} } \right|}} $$

(47)

where the peak stress ratio $\tilde{\eta }_{p} =$ $\tilde{\eta }_{c} {\text{exp}}\left( { - \tilde{k}\Psi } \right)$, in which $\tilde{k}$ is a material constant, defining the dependence of peak state of interface on material state. $\tilde{k}$ can be determined by fitting the peak stress ratio data points in the $\tilde{\eta }_{p} - \Psi$ plane. $K_{0}$ is a material constant, controlling the hardening behaviour of interface, which can be determined by fitting the tangential stress–displacement relation.

As the aim of this study is to develop a simple stress–dilatancy equation, details on how to identify model parameters will not be provided in a separate section, but instead described directly after each constitutive equation, as shown before. For more details on parameter identification, one can refer to Dafalias and Manzari [2] and Saberi et al. [24].