Fully coupled global equations for hydro-mechanical analysis of unsaturated soils

Abstract

Fully coupled global equations are proposed for enhancing the performance of Finite Element analysis of unsaturated soils. The governing equation describing mechanical equilibrium is formulated in terms of net stress, and in the mass conservation equation the contribution of this net stress in determining the change of degree of saturation is also included. The novelty of this paper is the development of new global finite element equations that can be used to find an approximate solution to these governing equations. The new equations have a mechanical term appearing in the flow matrix that is additional to the usual hydraulic term. This is in contrast to previous studies in which the coupling matrices ignore this effect. A performance study has been conducted for undrained footing problems, which shows that the additional mechanical term appearing in the flow matrix has a large influence on the accuracy of the numerical results.

Appendices

Appendix 1: supplementary to Eq. (3)

1.1 Fully coupled stress–strain relation

Equation (3) is the general form of the stress–strain equations, and it can be derived from constitutive models of unsaturated soils. In this paper, the fully coupled model by Zhou et al. [37, 38] is adopted and givens as:

$$ \begin{aligned} {\mathbf{D}}_{\text{ep}} & = {\varvec{\Gamma}}^{ - 1} \left( {{\mathbf{D}}_{\text{e}}^{ '} - \frac{{{\mathbf{D}}_{\text{e}}^{ '} {\text{a(a}}^{\text{T}} ){\mathbf{D}}_{\text{e}}^{ '} }}{{A + {\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} {\text{a}}}}} \right) \\ {\mathbf{W}}_{\text{ep}} & = - {\varvec{\Gamma}}^{ - 1} \left[ {\left( {\frac{{{\mathbf{D}}_{\text{e}}^{ '} {\text{a}}C}}{{A + {\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} {\text{a}}}} + \frac{{\partial {\varvec{\upsigma}}}}{{\partial S_{\text{r}} }}} \right)\frac{{\partial S_{\text{r}} }}{{\partial u_{\text{w}} }} + \frac{{\partial {\varvec{\upsigma}}'}}{{\partial u_{\text{w}} }}} \right] \\ {\varvec{\Gamma}} & = \frac{{\partial {\varvec{\upsigma}}'}}{{\partial {\varvec{\upsigma}}}} + \left( {\frac{{{\mathbf{D}}_{\text{e}}^{ '} {\text{a}}C}}{{A + {\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} {\text{a}}}} + \frac{{\partial {\varvec{\upsigma}}'}}{{\partial S_{\text{r}} }}} \right)\frac{{\partial S_{\text{r}} }}{{\partial {\varvec{\upsigma}}}} \\ \end{aligned} $$

(16)

where ${\mathbf{D}}_{\text{e}}^{ '}$ is the elastic stiffness matrix in terms of effective stress, $\partial S_{\text{r}} /\partial u_{\text{w}}$ and $\partial S_{\text{r}} /\partial {\varvec{\upsigma}}$ is determined by Eq. (8), $\partial S_{\text{r}} /\partial u_{\text{w}}$ is given in Eq. (22) and

$$ \begin{aligned} a & = \frac{\partial f}{{\partial {\boldsymbol{\sigma}^{\prime}}}} \\ A & = - \frac{\partial f}{{\partial p_{\text{c}}^{'} }}\frac{{\partial p_{\text{c}}^{'} }}{{\partial \varepsilon_{\text{v}}^{\text{p}} }}\frac{\partial f}{\partial p^{\prime}} \\ C & = \frac{\partial f}{{\partial S_{\text{r}} }} \\ \end{aligned} $$

where $\varepsilon_{\text{v}}^{\text{p}}$ is the plastic volumetric strain, and ${\boldsymbol{\sigma}^{\prime}}$ is Bishop’s effective stress is given as

$$ {\boldsymbol{\upsigma}^{\prime}} = {\varvec{\upsigma}} - {\mathbf{m}}S_{\text{e}} u_{\text{w}} $$

(17)

where ${\varvec{\upsigma}}$ is the net stress vector or total stress vector, ${\boldsymbol{\upsigma}^{\prime}}$ is Bishop’s effective stress vector, m is a column vector ($m = \left\{ {1,1,1,0,0,0} \right\}$), $u_{\text{w}}$ is the pore water pressure, and $S_{\text{e}}$ is the effective degree of saturation given as

$$ S_{\text{e}} = \frac{{S_{\text{r}} - S_{\text{r}}^{\text{res}} }}{{S_{\text{r}}^{ 0} - S_{\text{r}}^{\text{res}} }} $$

(18)

The yield function, $f$ is given as

$$ f = f\left( {{\boldsymbol{\sigma}^{\prime}}, S_{\text{r}} , p_{\text{c}}^{'} } \right) = \left( {\frac{q}{{M^{2} p}^{\prime}} + p^{\prime}} \right) - \left( {p_{\text{c}}^{'} } \right)^{\beta } = 0 $$

(19)

where $q$ is the deviator stress, $p^{\prime}$ is the effective mean stress, M is the stress ratio at the critical state, $p_{\text{c}}^{'}$ is the preconsolidation pressure at $S_{\text{r}} = S_{\text{r}}^{ 0}$, and it is taken as the hardening parameter, and

$$ \beta = \frac{\lambda - \kappa }{{\lambda - \left( {\lambda - \kappa } \right)\left( {1 - S_{\text{e}} } \right)^{{a_{1} }} - \kappa }} $$

(20)

where $\lambda$ is the elastoplastic compression index at the saturated state, $\kappa$ is elastic compression index, and $a_{1}$ is a fitting parameter that defines the variation of compression index with the degree of saturation.

The hardening parameter, $p_{\text{c}}^{'}$ can be determined from

$$ {\text{d}}p_{\text{c}}^{'} = {\mathbf{G}}_{\text{ep}} {\text{d}}{\varvec{\upvarepsilon}} + N_{\text{ep}} {\text{d}}u_{\text{w}} $$

(21)

where

$$ \begin{aligned} {\mathbf{G}}_{\text{ep}} & = \frac{{\partial p_{\text{c}}^{'} }}{{\partial \varepsilon_{\text{v}}^{\text{p}} }}\frac{\partial f}{\partial p^{\prime}}\left( {\frac{{{\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} + C\frac{{\partial S_{\text{r}} }}{{\partial {\varvec{\upsigma}}}}{\mathbf{D}}_{\text{ep}} }}{{A + {\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} {\text{a}}}}} \right) \\ N_{\text{ep}} & = \frac{{\partial p_{\text{c}}^{'} }}{{\partial \varepsilon_{\text{v}}^{\text{p}} }}\frac{\partial f}{\partial p^{\prime}}\frac{{C\left( {\frac{{\partial S_{\text{r}} }}{{\partial u_{\text{w}} }} + \frac{{\partial S_{\text{r}} }}{{\partial {\varvec{\upsigma}}}}{\mathbf{W}}_{\text{ep}} } \right)}}{{A + {\text{a}}^{\text{T}} {\mathbf{D}}_{\text{e}}^{ '} {\text{a}}}} \\ \end{aligned} $$

More details of Eq. (3) can be found in Sect. 2 of Zhang et al. [35].

Appendix 2: supplementary to Eq. (7)

2.1 Enhanced soil–water retention equation

The coupled behaviour of the degree of saturation can be expressed as Eq. (7). The first term $\left( {\frac{{\partial S_{\text{r}} }}{{\partial u_{\text{w}} }}} \right)$ adopted in this paper, is the hydraulic term of soil–water retention soil–water retention curve (SWRC), and stands for the variation of degree of saturation caused by suction change. The first term $\left( {\frac{{\partial S_{\text{r}} }}{{\partial u_{\text{w}} }}} \right)$ can be written as:

$$ \begin{aligned} & \frac{{\partial S_{\text{r}} }}{{\partial u_{\text{w}} }} = \frac{{\partial S_{\text{r}} }}{{\partial S_{\text{e}} }}\frac{{\partial S_{\text{e}} }}{{\partial u_{\text{w}} }} \\ &\quad = \left\{ {\begin{array}{*{20}c} {(S_{\text{r}}^{ 0} - S_{\text{r}}^{\text{res}} )\left( {\frac{{u_{\text{w}}^{\text{d}} }}{{u_{\text{w}} }}} \right)^{ - b} \left[ {\left( {\frac{{\partial S_{\text{e}}^{\text{d}} }}{{\partial u_{\text{w}} }}} \right)_{{u_{\text{w}} = u_{\text{w}}^{\text{d}} }} } \right],} & {{\text{if drying}} ({\text{d}}u_{\text{w}} < 0)} \\ {(S_{\text{r}}^{ 0} - S_{\text{r}}^{\text{res}} )\left( {\frac{{u_{\text{w}}^{\text{w}} }}{{u_{\text{w}} }}} \right)^{b} \left[ {\left( {\frac{{\partial S_{\text{e}}^{\text{w}} }}{{\partial u_{\text{w}} }}} \right)_{{u_{\text{w}} = u_{\text{w}}^{\text{w}} }} } \right],} & { {\text{if wetting}} ({\text{d}}u_{\text{w}} > 0)} \\ \end{array} } \right. \end{aligned} $$

(22)

where the main drying/wetting curve is defined by the VG model as:

$$ \left\{ {\begin{array}{*{20}c} {S_{\text{e}}^{\text{d}} = \left[ {1 + \left( {\frac{{ - u_{\text{w}} }}{{a_{\text{d}} }}} \right)^{{m_{\text{d}} }} } \right]^{{ - n_{\text{d}} }} } & {\left( {\text{for main drying}} \right)} \\ {S_{\text{e}}^{\text{w}} = \left[ {1 + \left( {\frac{{ - u_{\text{w}} }}{{a_{\text{w}} }}} \right)^{{m_{\text{w}} }} } \right]^{{ - n_{\text{w}} }} } & {\left( {\text{for main wetting}} \right)} \\ \end{array} } \right. $$

$$ \left\{ {\begin{array}{*{20}cc} &{\frac{{\partial S_{\text{e}}^{\text{d}} }}{{\partial u_{\text{w}} }} = n_{\text{d}} \left[ {1 + \left( {\frac{{ - u_{\text{w}} }}{{a_{\text{d}} }}} \right)^{{m_{\text{d}} }} } \right]\left( {\frac{{m_{\text{d}} }}{{a_{\text{d}} }}} \right)\left( {\frac{{ - u_{\text{w}} }}{{a_{\text{d}} }}} \right)^{{m_{\text{d}} - 1}} } \\ & {\left( {\text{gredient of main drying branch}} \right)} \\ &{\frac{{\partial S_{\text{e}}^{\text{w}} }}{{\partial u_{\text{w}} }} = n_{\text{w}} \left[ {1 + \left( {\frac{{ - u_{\text{w}} }}{{a_{\text{w}} }}} \right)^{{m_{\text{w}} }} } \right]\left( {\frac{{m_{\text{w}} }}{{a_{\text{w}} }}} \right)\left( {\frac{{ - u_{\text{w}} }}{{a_{\text{w}} }}} \right)^{{m_{\text{w}} - 1}} } \\ & {\left( {\text{gredient of main wetting branch}} \right)} \\ \end{array} } \right. $$

$$ \left\{ {\begin{array}{*{20}c} {u_{\text{w}}^{\text{d}} = - a_{\text{d}} \left[ {\left( {\frac{{S_{\text{r}} - S_{\text{r}}^{\text{res}} }}{{S_{\text{r}}^{ 0} - S_{\text{r}}^{\text{res}} }}} \right)^{{ - \frac{1}{{n_{\text{d}} }}}} - 1} \right]^{{\frac{1}{{m_{\text{d}} }}}} } & {\left( {\text{for drying}} \right)} \\ {u_{\text{w}}^{\text{w}} = - a_{\text{w}} \left[ {\left( {\frac{{S_{\text{r}} - S_{\text{r}}^{\text{res}} }}{{S_{\text{r}}^{ 0} - S_{\text{r}}^{\text{res}} }}} \right)^{{ - \frac{1}{{n_{\text{w}} }}}} - 1} \right]^{{\frac{1}{{m_{\text{w}} }}}} } & {\left( {\text{for wetting}} \right)} \\ \end{array} } \right. $$

where $a_{\text{d}}$, $m_{\text{d}}$ and $m_{\text{d}}$ are fitting parameters for the main drying curve, $a_{\text{w}}$, $m_{\text{w}}$ and $m_{\text{w}}$ are fitting parameters for the main wetting curve, and $b$ is a fitting parameter for non-linear hysteresis.

The second term in Eq. (7), $\frac{{\partial S_{r} }}{{\partial {\varvec{\upsigma}}}}$ is the mechanical term of the SWRC and it can be determined by Eq. (8), where

$$ \begin{aligned} {\mathbf{C}}_{{{\text{v}\upsigma}}} & = \left\{ {\begin{array}{*{20}c} {\frac{{\frac{\kappa }{{p^{\prime}\left( {1 + e} \right)}}\frac{\partial p^{\prime}}{\partial p} + J\frac{\partial f}{\partial p^{\prime}}\frac{\partial p^{\prime}}{\partial p}}}{{1 - \frac{\kappa }{{p^{\prime}\left( {1 + e} \right)}}\frac{{\partial p^{\prime}}}{{\partial S_{\text{r}} }}\frac{{\partial S_{\text{r}} }}{{\partial \varepsilon_{{{\text{v}\upsigma}}} }} - J\left( {\frac{\partial f}{\partial \lambda }\frac{\partial \lambda }{{\partial S_{\text{r}} }} + \frac{\partial f}{\partial p^{\prime}}\frac{\partial p^{\prime}}{{\partial S_{\text{r}} }}} \right)\frac{{\partial S_{\text{r}} }}{{\partial \varepsilon_{{{\text{v}\upsigma}}} }}}}} \\ {\frac{{J\frac{\partial f}{\partial q}}}{{1 - \frac{\kappa }{{p^{\prime}\left( {1 + e} \right)}}\frac{{\partial p^{\prime}}}{{\partial S_{\text{r}} }}\frac{{\partial S_{\text{r}} }}{{\partial \varepsilon_{{{\text{v}\upsigma}}} }} - J\left( {\frac{\partial f}{\partial \lambda }\frac{\partial \lambda }{{\partial S_{\text{r}} }} + \frac{\partial f}{\partial p^{\prime}}\frac{\partial p^{\prime}}{{\partial S_{\text{r}} }}} \right)\frac{{\partial S_{\text{r}} }}{{\partial \varepsilon_{{{\text{v}\upsigma}}} }}}}} \\ \end{array} } \right\}\\ &\quad \left\{ {\left( {\frac{\partial p^{\prime}}{{\partial {\varvec{\upsigma}}}}} \right)^{\text{T}} ,\varvec{ }\left( {\frac{\partial q}{{\partial {\varvec{\upsigma}}}}} \right)^{\text{T}} } \right\} \end{aligned} $$

(23)

where $J = - \left( {p_{\text{c}}^{'} \frac{\partial f}{{\partial p_{\text{c}}^{'} }}\frac{1 + e}{{\lambda_{0} - \kappa }}} \right)^{ - 1}$.

Appendix 3: Newton–Raphson iteration scheme

The backward Euler method with Newton–Raphson iterations can be summarised as following steps:

Enter with the current displacements and pore pressures ${\mathbf{X}}_{t} = \left\{ {{\mathbf{U}}_{t} ,{\mathbf{U}}_{{{\text{w}}t}} } \right\}$, the corresponding derivatives ${\dot{\mathbf{X}}}_{t}$, the current time substep size $h$, the iteration tolerance ITOL.

Compute estimate of new displacements and pore pressures and the corresponding rates using

$$ \begin{aligned} {\dot{\mathbf{X}}}_{t + h}^{0} & = {\dot{\mathbf{X}}}_{t} \\ {\tilde{\mathbf{X}}}_{t + h}^{0} & = {\mathbf{X}}_{n - 1} + h{\dot{\mathbf{X}}}_{t} \end{aligned} $$

Set the initial the relative change in ${\tilde{\mathbf{X}}}$ as:

$$ \alpha^{0} = h{\dot{\mathbf{X}}}_{t} /{\tilde{\mathbf{X}}}_{t + h}^{0} $$

Repeat steps 5 to 7 until $\alpha^{i} \le ITOL$, where ITOL is a given tolerance.

Compute the residual vector

$$ {\mathbf{R}}^{i} = \left[ {\begin{array}{*{20}c} {\frac{{{\mathbf{F}}_{t + h}^{\text{ext}} - {\mathbf{F}}_{t + h}^{\text{int}} }}{h}} \\ {{\mathbf{Q}}_{t + h}^{\text{ext}} - {\mathbf{Q}}_{t + h}^{\text{int}} } \\ \end{array} } \right] $$

and solve for $\delta {\dot{\mathbf{X}}}^{i}$ using

$$ \delta {\dot{\mathbf{X}}}^{i} = \left[ {{\mathbf{C}}_{\text{ep}} + h{\mathbf{K}}} \right]^{ - 1} {\mathbf{R}}^{i} $$

where ${\mathbf{F}}_{t + h}^{\text{int}}$, ${\mathbf{Q}}_{t + h}^{\text{int}}$, ${\mathbf{C}}_{\text{ep}} = \left[ {\begin{array}{*{20}c} {{\mathbf{K}}_{\text{ep}} } & {\mathbf{L}} \\ {{\mathbf{L}^{\prime}}} & {\mathbf{S}} \\ \end{array} } \right]$ and ${\mathbf{K}} = \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & {\mathbf{H}} \\ \end{array} } \right]$ are evaluated at ${\tilde{\mathbf{X}}}_{t + h}^{i - 1}$.

Update the displacements and pore pressures and the corresponding rates to

$$ \begin{aligned} {\dot{\mathbf{X}}}_{t + h}^{i} & = {\dot{\mathbf{X}}}_{t + h}^{i - 1} + \delta {\dot{\mathbf{X}}}^{i} \\ {\tilde{\mathbf{X}}}_{t + h}^{i} & = {\mathbf{X}}_{t} + h{\dot{\mathbf{X}}}_{t + h}^{i} \end{aligned} $$

Compute convergence criterion

$$ \alpha^{i} = h{\dot{\mathbf{X}}}^{i} /{\tilde{\mathbf{X}}}_{t + h}^{i} $$

If $\alpha^{i} \le {\text{ITOL}}$ then go to step 8.

Exit with displacements and pore pressures, ${\tilde{\mathbf{X}}}_{t + h} = {\tilde{\mathbf{X}}}_{t + h}^{i}$, and there rates, ${\dot{\mathbf{X}}}_{t + h} = {\dot{\mathbf{X}}}_{t + h}^{i}$, at time $h + t$.