Coupled solid-fluid FE-analysis of an embankment dam

Abstract

A coupled solid-fluid FE-model for partially saturated soils, characterized by modeling the soil as a three-phase material consisting of a deformable soil skeleton and the fluid phases water and air, is reviewed briefly. As a constitutive model for the soil skeleton, the well-known Barcelona Basic model (BBM) is employed, which is formulated in terms of net stress and matric suction. For the BBM, a computationally efficient return mapping algorithm is proposed, which only requires the solution of a scalar nonlinear equation at the integration point level. The coupled FE-model is applied to the coupled transient numerical simulation of the water flow and the deformations and stresses in an embankment dam.

Keywords

multi-phase model / unsaturated soil model / Barcelona Basic model (BBM) / return mapping algorithm / embankment dam

Introduction

In an embankment dam, the water flow through the dam is coupled with the deformations and stresses of the dam. To consider such solid-fluid coupling effects in a numerical model, a two-phase or three-phase model for the soil, consisting of a deformable soil skeleton and the two interstitial fluid phases water and air, is required. The two-phase and three-phase model differ by the fact that in the former, the pressure in the air phase is assumed to stay at atmospheric conditions, whereas in the latter, variations of the air pressure are taken into account.

Coupled analyses of embankment dams were already conducted, e.g., by participants of the fifth ICOLD Benchmark Workshop on numerical analysis of dams [1]. The main difference between those previous analyses and the present analysis refers to the constitutive model for the soil. Whereas the analyses, documented in Ref. [1], are based on soil models, formulated solely in terms of the effective stress, and it is commonly accepted now that models for partially saturated soils have to be formulated in terms of two independent stress state variables. Probably, the best-known constitutive model for unsaturated soils is the Barcelona Basic model (BBM) [2]. For the latter, an optimized return mapping algorithm will be presented, which is characterized by analytical integration of the hardening law and by solving only a single nonlinear scalar equation at the integration point level. This improves the accuracy and robustness of the stress integration and reduces the required computing time.

The paper is organized as follows. At first, the basic equations of the three-phase formulation and of the Barcelona Basic model are reviewed briefly. Then, the algorithmic treatment of the BBM within the framework of the three-phase formulation is outlined. Finally, the presented coupled FE-model is applied to the numerical simulation of the water flow through an earth dam, showing the impact of the water flow through the dam on the deformations and stresses of the dam.

Three-phase model

In the present three-phase formulation, the soil is considered as a three-phase mixture, consisting of a deformable soil skeleton with incompressible soil grains and the fluid phases water and air. Assuming quasi-static conditions, the equilibrium equations for the soil are given as

(1)-\mathrm{d}\mathrm{i}\mathrm{v}\mathit{\sigma }+\rho \mathit{g}=\mathbf{0}.

In Eq. (1), BoldItalic and ρ denote the total stress tensor and the density of the three-phase mixture, respectively, and BoldItalic is the vector of gravitational acceleration. Compressive stresses are defined as positive quantities.

The mass balance equation for a fluid phase f is derived as [3]

(2)n{\displaystyle \frac{{\mathrm{d}}^s(S_f{\rho }^f)}{\mathrm{d}t}}+S_f{\rho }^f\mathrm{d}\mathrm{i}\mathrm{v}{\mathit{v}}_s+\mathrm{d}\mathrm{i}\mathrm{v}(n{S}_f{\rho }^f{\mathit{v}}_{fs})=\mathrm{0.}

In Eq. (2), n, S_f, ρ^f, and BoldItalic_s denote the porosity of the three-phase mixture, the degree of saturation of fluid phase f, the intrinsic density of fluid phase f, and the velocity of the solid phase s, respectively. Setting f = a and f = w gives the mass balance equation for the air phase and the water phase, respectively. d^s/dt represents the material time derivative with respect to the solid phase, and BoldItalic_fs = BoldItalic_f-BoldItalic_s is the velocity of the fluid phase f relative to the velocity of the solid phase s. It is governed by Darcy’s law

(3)n{S}_f{\mathit{v}}_{fs}={\displaystyle \frac{\mathrm{1}}{{\rho }^{f}g}}{\mathit{k}}_f(-\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{p}^f+{\rho }^f\mathit{g})

with BoldItalic_f as the permeability tensor, g=\Vert \mathit{g}\Vert, and p^f as the pressure of fluid phase f.

The capillary pressure or matric suction is defined as the difference between the pressure of the air phase and the pressure of the water phase, i.e.,

(4)p^c=p^a-p^w.

The constitutive relationship between the capillary pressure and the effective degree of water saturation S_w,eff can be approximated by [4]

(5)S_{w,\mathrm{e}\mathrm{f}\mathrm{f}}={\displaystyle \frac{S_w-S_{w,r}}{S_{w,s}-S_{w,r}}}={\left[\mathrm{1}+{\left({\displaystyle \frac{p^c}{p_b^c}}\right)}^n\right]}^{-m}.

In Eq. (5), p_b^c denotes the air entry value and n and m are parameters for fitting Eq. (5) to experimental data. The effective degree of water saturation 0≤S_w,eff≤1 is computed in terms of the maximum and residual degree of water saturation, S_w,s and S_w,r, respectively.

For partially saturated soils also, the coefficients of the permeability tensor BoldItalic_f depend on the degree of saturation of the respective fluid phase. Various empirical approximations for the relative permeability in terms of the degree of water saturation are reported in Ref. [4,5]. In the present study, isotropic permeability is assumed, i.e., {\mathit{k}}_f=k^{of}{k}_{r,f}BoldItalic, where k^of denotes the permeability for saturated conditions, and BoldItalic denotes the identity matrix. The relative permeabilities k_{r, f} of the fluid phases air and water account for partial saturation and are assumed as [6]

(6)\begin{array}{l}{k}_{r,a}=\sqrt{\mathrm{1}-S_{w,\mathrm{e}\mathrm{f}\mathrm{f}}}{[\mathrm{1}-S_{w,\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{1}/m}]}^{\mathrm{2}m},\hfill \\ k_{r,w}=\sqrt{S_{w,\mathrm{e}\mathrm{f}\mathrm{f}}}{[\mathrm{1}-{(\mathrm{1}-S_{w,\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{1}/m})}^m]}^{\mathrm{2}}.\hfill \end{array}

The constitutive equations for the fluid phases water and air relate the intrinsic density to the respective fluid pressure, i.e.,

(7)\begin{array}{ll}{\dot{\rho }}^w={\rho }_{\mathrm{0}}^w{C}^w{\dot{\rho }}^w,& {\dot{\rho }}^a={\displaystyle \frac{\mathrm{1}}{R_{a}T}}{\dot{\rho }}^a\end{array}

with C^w denoting the compressibility of the water phase, {\rho }_{\mathrm{0}}^w the density of the water phase at p^w = 0, and R_a and T representing the gas constant and the absolute temperature, respectively.

The constitutive relations for the soil skeleton will be described in the subsequent section. For describing the basic features of partially saturated soils, the constitutive relations have to be formulated in terms of two independent stress state variables. Among several proposed definitions of stress state variables (e.g., Ref. [7-9]), the use of net stress

(8){\mathit{\sigma }}^{"}=\mathit{\sigma }-p^a\mathit{I},

i.e., the total stress in excess of the pore air pressure, and matric suction p^c, and the use of the generalized effective stress tensor

(9){\mathit{\sigma }}^{\prime }=\mathit{\sigma }-p^a\mathit{I}+S^w(p^a-p^w)\mathit{I},

and matric suction p^c are popular choices for constitutive models for partially saturated soils. A claimed advantage of the stress states variables {\mathit{\sigma }}^{"} and p^c is given by the fact that they degenerate to the total stress and to the negative pore water pressure in the case of atmospheric pore air pressure. However, using these two stress variables, the transition from a partially saturated to a saturated state requires a switch of the stress state variables, because in the case of water saturated conditions, {\mathit{\sigma }}^{"} does not contain the effective stress tensor, proposed by Terzaghi [10], as a special case. The choice of {\mathit{\sigma }}^{\prime }and p^c as stress state variables is advantageous, because for the limiting case of a water saturated soil, characterized by S^w = 1, Eq. (9) degenerates to the well-known effective stress tensor for water saturated soils. Thus, material models for partially saturated soils, formulated in terms of {\mathit{\sigma }}^{\prime } and p^c, allow a straightforward transition from partially saturated to saturated conditions.

Unsaturated soil model

For the present three-phase formulation, the Barcelona Basic model (BBM) is employed as constitutive model for describing unsaturated soil behavior [2]. Although since the publication of the pioneering paper [2], this model was extended in many subsequent papers, and it is employed here in its original version, because the latter was agreed as the basis for extensive benchmark activities within the framework of the MUSE network (http:muse.dur.ac.uk).

The BBM is formulated in terms of the net stress tensor (8) and matric suction p^c. For stress states located within the elastic domain, enclosed by the yield surface, the elastic volumetric and deviatoric strain rates are given as

(10a){\dot{ϵ}}_v^e={\displaystyle \frac{\kappa }{\mathrm{1}+e}}{\displaystyle \frac{{\dot{p}}^{"}}{p^{"}}}+{\displaystyle \frac{{\kappa }_s}{\mathrm{1}+e}}{\displaystyle \frac{{\dot{p}}^c}{p^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}},

(10b){\dot{e}}_{ij}^e={\displaystyle \frac{{\dot{s}}_{ij}}{\mathrm{2}G}}

with the material parameters κ and κ_s, representing the elastic stiffness for changes of the mean net pressure p^{"}=({\mathit{\sigma }}^{"}:\mathit{I})/\mathrm{3} and for changes of matric suction p^c, respectively. e, p_atm, {\dot{s}}_{ij}, and G denote the void ratio, the atmospheric air pressure, the deviatoric stress rate, and the shear modulus, respectively. It follows from Eq. (10a) that the elastic volumetric strain rate {\dot{ϵ}}_v^e depends on both the mean net pressure p^{"} and matric suction p^c. Figure 1 shows the volumetric behavior of the BBM in an \ln p^{"}-e diagram. Within the elastic domain, the stress point (p^{"},e) lies on the unloading-reloading line (URL) with the slope κ and for isotropic plastic conditions on the isotropic compression line (ICL) with the suction-dependent slope λ(p^c). The intersection point of the URL and the ICL is the preconsolidation pressure p_{\mathrm{0}}^{"}. The ICL is defined by the slope λ(p^c) and the void ratio N(p^c)-1 at p^{"}=\mathrm{1} with N(p^c) denoting the respective specific volume. From geometric considerations in Fig. 1, we have

(11)e=N(p_c)-\mathrm{1}-\lambda (p^c)\ln p_{\mathrm{0}}^{"}+\kappa \ln \left({\displaystyle \frac{p_{\mathrm{0}}^{"}}{p^{"}}}\right)

from which

(12)p_{\mathrm{0}}^{"}=\exp \left({\displaystyle \frac{-N(p^c)+(\mathrm{1}+e)}{\kappa -\lambda (p^c)}}\right)p^{"\frac{\kappa }{\kappa -\lambda (p^c)}}

is obtained.

The value of N(p^c) follows from the requirement that the void ratio at point C in Fig. 1 is the same for the paths D → C and A → B → C, i.e.,

(13)N(p^c)-\mathrm{1}-\lambda (p^c)\ln p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}=N(\mathrm{0})-\mathrm{1}-\lambda (\mathrm{0})\ln p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}-{\kappa }_s\ln \left({\displaystyle \frac{p^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}{p_{\mathrm{a}\mathrm{t}\mathrm{m}}}}\right)

from which

(14)N(p^c)=N(\mathrm{0})-(\lambda (\mathrm{0})-\lambda (p^c))\ln p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}-{\kappa }_s\ln \left({\displaystyle \frac{p^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}{p_{\mathrm{a}\mathrm{t}\mathrm{m}}}}\right)

follows, where N(0) and p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"} denote material parameters. Since the LC yield curve is a straight line for {(p_{\mathrm{0}}^{"})}^{\ast }=p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}, the wetting path B→C is elastic. It follows from Eq. (10) that the change of the void ratio for the wetting path B→C is obtained as

(15)\mathrm{\Delta }{e}_{BC}=-{\displaystyle {\int }_B^C(\mathrm{1}+e){\dot{ϵ}}_v}=-{\kappa }_s\ln \left({\displaystyle \frac{p^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}{p_{\mathrm{a}\mathrm{t}\mathrm{m}}}}\right).

The yield surface is defined as (Fig. 2)

(16)J_{\mathrm{2}}-{\displaystyle \frac{M^{\mathrm{2}}}{\mathrm{3}}}(p^{"}+p_s^{"})(p_{\mathrm{0}}^{"}-p^{"})=\mathrm{0}

with the second invariant of the deviatoric stress tensor J₂ = s_ijs_ij/2 and

(17a)p_s^{"}(p^c)=k_s{p}^c,

(17b)p_{\mathrm{0}}^{"}=p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}{\left({\displaystyle \frac{{(p_{\mathrm{0}}^{"})}^{\ast }}{p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}}}\right)}^{\frac{\lambda (\mathrm{0})-\kappa }{\lambda (p^c)-\kappa }},

(18)\lambda (p^c)=\lambda (\mathrm{0})[(\mathrm{1}-r)\exp (-\beta p^c)+r].

In Eq. (16), M defines the slope of the critical state line. The preconsolidation pressure p_{\mathrm{0}}^{"} and the one for saturated conditions {(p_{\mathrm{0}}^{"})}^{\ast } are located on the so-called loading collapse yield curve (LC curve) according to Eq. (17b). This curve is the intersection of the yield surface with the plane J₂ = 0 for positive values of p″ (see Fig. 2). For negative values of p″, the intersection with the plane J₂ = 0 is given by p_s^{"} according to (17a) with the material parameter k_s describing the increase in cohesion due to matric suction (Fig. 2). p_s^{"} and p_{\mathrm{0}}^{"} both depend on matric suction according to Eq. (17). p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"} serves as a reference pressure such that for {(p_{\mathrm{0}}^{"})}^{\ast }=p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}, (17b) degenerates to p_{\mathrm{0}}^{"}=p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}. λ(p^c) describes the soil stiffness during plastic loading in a hydrostatic test for a given matric suction in terms of the respective stiffness λ(0) at saturated conditions and the material parameters r and β.

The plastic strain rate is determined from the non-associated flow rule

(19){\dot{ϵ}}_{ij}^p=\dot{\gamma }{\displaystyle \frac{\partial g}{\partial {\sigma }_{ij}^{"}}}

(20)g=\alpha J_{\mathrm{2}}-{\displaystyle \frac{M^{\mathrm{2}}}{\mathrm{3}}}(p^{"}+p_s^{"})(p_{\mathrm{0}}^{"}-p^{"}),

where

(21)\alpha ={\displaystyle \frac{M(M-\mathrm{9})(M-\mathrm{3})}{\mathrm{9}(\mathrm{6}-M)}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\displaystyle \frac{\kappa }{\lambda (\mathrm{0})}}}}

is a constant. The hardening law relates the rate of the preconsolidation pressure at saturated conditions {(p_{\mathrm{0}}^{"})}^{\ast }, which serves as the hardening parameter, to the volumetric plastic strain rate {\dot{ϵ}}_v^p by

(22){\left({\dot{p}}_{\mathrm{0}}^{"}\right)}^{\ast }={(p_{\mathrm{0}}^{"})}^{\ast }{\displaystyle \frac{\mathrm{1}+e}{\lambda (\mathrm{0})-\kappa }}{\dot{ϵ}}_v^p.

Equation (22) describes the evolution of the yield surface. The latter is shown for two different values of {(p_{\mathrm{0}}^{"})}^{\ast } in Fig. 2.

The suction increase yield surface used in Ref. [2] is not included in the present implementation.

Finite element formulation

The finite element formulation of the three-phase model relies on 1) the weak formulations of the governing Eqs. (1) and (2), the latter for both fluid phases, water, and air and on the 2) subdivision of the (spatial) domain under consideration into a number of finite elements. Within a single finite element, the primary variables, consisting of the displacements of the soil skeleton, matric suction, and the pressure of the air phase, are interpolated from the respective nodal values by suitable shape functions. This results in the coupled system of equations for the Newton-type iteration step i within the time step Δt_n+1:

(23){\mathit{C}}_{n+\mathrm{1}}^{(i)}\mathrm{\Delta }{\dot{\mathit{q}}}_{i+\mathrm{1},n+\mathrm{1}}+{\mathit{K}}_{n+\mathrm{1}}^{(i)}\mathrm{\Delta }{\mathit{q}}_{i+\mathrm{1},n+\mathrm{1}}={\mathit{f}}_{n+\mathrm{1}}^{(i)}

with BoldItalic and BoldItalic denoting tangent matrices and BoldItalic denoting a residual vector resulting from the spatial discretization and BoldItalic representing the global vector of nodal variables. Numerical integration in the time domain by means of the implicit Euler method yields

(24)\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_{n+\mathrm{1}}}}{\mathit{C}}_{n+\mathrm{1}}^{(i)}+{\mathit{K}}_{n+\mathrm{1}}^{(i)}\right]\mathrm{\Delta }{\mathit{q}}_{i+\mathrm{1},n+\mathrm{1}}={\mathit{f}}_{n+\mathrm{1}}^{(i)},

where ΔBoldItalic_i+1,n+1 represents the contribution of iteration step i + 1 to the global vector of the nodal variables at t_n+1

(25){\mathit{q}}_{n+\mathrm{1}}^{(i+\mathrm{1})}={\mathit{q}}_{n+\mathrm{1}}^{(i)}+\mathrm{\Delta }{\mathit{q}}_{i+\mathrm{1},n+\mathrm{1}}.

Stress update algorithm

Within the framework of an FE-analysis, the BBM, described in Section 3, is employed for determining the stresses at a given time instant t_n+1, provided that the equilibrium state at time instant t_n is known, and, in addition, for the time instant t_n+1, an estimate of the nodal variables (25) is known.

A popular stress update algorithm is the return mapping algorithm. Basically, it consists of two steps. In the first step from the known estimate of the primary variables at t_n+1, a trial stress, also denoted as elastic predictor, is computed by assuming elastic material behavior for the current time step. If the yield condition is violated by the trial stress, then in the second step, the so-called corrector step, the trial stress is returned to the yield surface.

Since the BBM is characterized by 1) a single yield surface (16), 2) a single plastic potential (20) only depending on the first and second invariant of the net stress tensor, I_{\mathrm{1}}^{"} and J₂, 3) a single hardening parameter that can be expressed for the corrector step in terms of only the mean net pressure using (12) and (17b), and 4) a relation for the deviatoric stress rate (10b) depending on the deviatoric elastic strain rate through the shear modulus G, the return mapping algorithm can be formulated in a very efficient manner in terms of a single scalar nonlinear equation for the mean net pressure p^{"}=I_{\mathrm{1}}^{"}/\mathrm{3} at t_n+1. A similar algorithm was proposed by Ref. [11] for the modified cam clay model for water saturated soils.

The flow rule (19) can be split into a volumetric and deviatoric part

(26)\begin{array}{ll}{\dot{ϵ}}_v^p=\dot{\gamma }\mathrm{3}{\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}^{"}}},& {\dot{e}}_{ij}^p=\dot{\gamma }{\displaystyle \frac{\partial g}{\partial J_{\mathrm{2}}}}{s}_{ij}.\end{array}

Backward Euler integration of (26) yields

(27a)\mathrm{\Delta }{ϵ}_v^p=\gamma \mathrm{3}{\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}^{"}}}=\mathrm{\Delta }{ϵ}_v-\mathrm{\Delta }{ϵ}_v^e,

(27b)e_{ij}^p=\gamma {\displaystyle \frac{\partial g}{\partial J_{\mathrm{2}}}}{s}_{ij}=\mathrm{\Delta }{e}_{ij}-\mathrm{\Delta }{e}_{ij}^e,

where \gamma =\dot{\gamma }\mathrm{\Delta }t with t = t_n+1-t_n. From Eq. (10b), we have

(28)s_{ij}=\mathrm{2}G{e}_{ij}^e=\mathrm{2}G[(e_{ij}-e_{ij,n}^p)-(e_{ij}^p-e_{ij,n}^p)]=s_{ij}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}-\mathrm{2}G\mathrm{\Delta }{e}_{ij}^p.

Note that quantities with the subscript n refer to the converged values at t_n, whereas all other quantities refer to the current values at t_n+1.

Using Eqs. (27) in Eq. (28) gives

(29)\left(\mathrm{1}+\gamma {\displaystyle \frac{\partial g}{\partial J_{\mathrm{2}}}}\mathrm{2}G\right)s_{ij}=s_{ij}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}.

Since the term enclosed by the brackets is a scalar quantity, s_ij and s_{ij}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}} differ only by a scalar factor. Hence, from Eq. (29), we have

(30){\left(\mathrm{1}+\gamma {\displaystyle \frac{\partial g}{\partial J_{\mathrm{2}}}}\mathrm{2}G\right)}^{\mathrm{2}}{J}_{\mathrm{2}}=J_{\mathrm{2}}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}.

Making use of \gamma ={\displaystyle \frac{\mathrm{\Delta }{ϵ}_v-\mathrm{\Delta }{ϵ}_v^e}{{\displaystyle \frac{\mathrm{3}\partial g}{\partial I_{\mathrm{1}}^{"}}}}}, resulting from Eq. (27a), yields

(31){\left(\mathrm{3}{\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}^{"}}}-{\displaystyle \frac{\partial g}{\partial J_{\mathrm{2}}}}\mathrm{2}G(\mathrm{\Delta }{ϵ}_v^e-\mathrm{\Delta }{ϵ}_v)\right)}^{\mathrm{2}}{J}_{\mathrm{2}}-{\left(\mathrm{3}{\displaystyle \frac{\partial g}{\partial I_{\mathrm{1}}^{"}}}\right)}^{\mathrm{2}}{J}_{\mathrm{2}}^{\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}=\mathrm{0.}

In Eq. (31), the incremental volumetric strain \mathrm{\Delta }{ϵ}_v is known from the current estimate of the displacements at t_n+1 and the converged values of the displacements at t_n. J₂, and \mathrm{\Delta }{ϵ}_v^e in Eq. (31) can be replaced by

(32)J_{\mathrm{2}}={\displaystyle \frac{M^{\mathrm{2}}}{\mathrm{3}}}(p^{"}+p_s^{"})(p_{\mathrm{0}}^{"}-p^{"}),

and

(33)\mathrm{\Delta }{ϵ}_v^e={\displaystyle \frac{\kappa }{\mathrm{1}+e}}\ln \left({\displaystyle \frac{p^{"}}{p_n^{"}}}\right)+{\displaystyle \frac{{\kappa }_s}{\mathrm{1}+e}}\ln \left({\displaystyle \frac{p^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}{p_n^c+p_{\mathrm{a}\mathrm{t}\mathrm{m}}}}\right)

following from Eqs. (16) and (10a). The rate of the void ratio is given by

(34)\dot{e}=-(\mathrm{1}+e){\dot{ϵ}}_v.

Integration of Eq. (34) yields the value of the void ratio at t_n+1

(35)e=(\mathrm{1}+e_n)\exp (-\mathrm{\Delta }{ϵ}_v)-\mathrm{1.}

Equation (31) together with Eqs. (32), (33), and (12) represents a nonlinear scalar equation for the unknown I_{\mathrm{1}}^{"} (or p^{"}=I_{\mathrm{1}}^{"}/\mathrm{3}), which can be solved, e.g., by the Newton method.

Once p″ has been determined from this equation, it is inserted into Eq. (12), yielding p_{\mathrm{0}}^{"}, and the latter into the recast form of Eq. (17b), yielding the hardening parameter:

(36){(p_{\mathrm{0}}^{"})}^{\ast }=p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}{\left({\displaystyle \frac{p_{\mathrm{0}}^{"}}{p_{\mathrm{r}\mathrm{e}\mathrm{f}}^{"}}}\right)}^{\frac{\lambda (p^c)-\kappa }{\lambda (\mathrm{0})-\kappa }}.

The Newton iteration for solving Eq. (31) together with Eq. (32), (33), and (12) can be started by the trial values p^{"\mathrm{T}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}} and BoldItalic^Trial, which are computed by assuming elastic material behavior for the current time step. However, an improved trial value for the mean net stress p^{"\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}} will be obtained by projecting p^{"\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}} on the isotropic compression line by

(37)p^{"\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}=\exp \left({\displaystyle \frac{M(p^c)-(\mathrm{1}+e)}{\lambda (p^c)}}\right).

p^{"\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}} is only used as an improved trial stress, when holds. A graphical interpretation of Eq. (37) is contained in Fig. 3.

Application

It is known from reported cases that even relatively small homogeneous dams may be affected by damage [1,12]. In this context, an important role plays the magnitude and distribution of the pore water pressure after construction of the dam, which depends on the permeability and the degree of water saturation of the dam material [13]. For a realistic analysis, the material parameters have to be determined for the employed dam material, and the construction process must be taken into account in the numerical analysis. Ideally, for the construction process, numerical results should be compared with measurement data and, if necessary, the material parameters should be calibrated accordingly.

In this section, the application of the three-phase model for unsaturated soils with the BBM as the constitutive model for the soil skeleton is demonstrated by a simplified coupled solid-fluid FE-analysis of water flow through a homogeneous earth dam. The geometric properties, shown in Fig. 4, are taken from an example problem outlined in Ref. [5]. The employed hydraulic parameters, taken from Ref. [14], are summarized in Table 1. The material parameters for the BBM, taken from Ref. [2], are provided in Table 2.

Figure 5 shows the FE-mesh and the applied boundary conditions. The foundation of the dam is considered by the respective boundary conditions at the bottom of the discretized domain. The displacements are fully constrained along the foundation, and for the undrained part of the base, an impermeable boundary is assumed. For the drained part, a permeable boundary is considered by applying a mass flux with a pressure dependent velocity

(38)\begin{array}{lll}{v}_n=k_{sc}{p}^w\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\hfill & p^w>\mathrm{0},\hfill \\ v_n=\mathrm{0}\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\hfill & p^w\le \mathrm{0},\hfill \end{array}

where p^w is the water pressure at the boundary, and k_sc is a sufficiently large seepage coefficient to approximately enforce the requirement of a zero water pressure for a freely draining surface. Similar boundary conditions for the fluid phase are applied for the free surfaces, i.e., for the upstream slope above the water level, the crest and the downstream slope.

The numerical analysis is subdivided into several parts. In the first part, as the initial condition, a constant value of matric suction of 100 kPa, corresponding to an initial degree of water saturation of S^w = 0.729, is assumed for the dam body and the primary stresses due to dead load are computed presuming elastic response. In the second part of the analysis, the net stresses, the void ratio, and the matric suction, computed in the first part, are taken as initial values; whereas the displacements are set equal to zero, and the matric suction at the upstream boundary is reduced to zero by specifying the respective boundary conditions. In the next part of the analysis, the transient seepage flow due to a water table of H = 10 m is computed until steady-state conditions are attained.

Figure 6 shows the distribution of matric suction for selected time instants. The isoline p^c = 0 represents the saturation front. Note that negative values of p^c indicate the pore water pressure in water saturated regions (this follows from setting p^a = 0 in Eq. (4)). Ten days after the impoundment, only a small region close to the upstream slope is affected by the water flow. With advancing time, the water saturated domain is propagating until steady-state conditions are attained. The transient evolution of the inflow and outflow of water is shown in Fig. 7. Whereas the inflow of water through the upstream face is decreasing with increasing time, the outflow through the drainage is increasing until steady-state conditions are attained with nearly identical values of the inflow and outflow.

Soils are characterized by complex material behavior during wetting. On one hand, the increase of the degree of water saturation yields an increase of the self weight of the soil, and hence, the total stresses are increasing. On the other hand, matric suction is decreasing, resulting in decreasing generalized effective stresses. For water saturated conditions, the soil skeleton is unloaded elastically due to an increase of the water pressure, resulting in volumetric expansion. For wetting during partial saturation, the soil strength depends on matric suction, and in this case, either elastic or plastic material behavior may occur. In case of elastic material behavior, the BBM predicts volumetric expansion, whereas in case of plastic material behavior, the volumetric behavior depends on the location of the stress state with respect to the critical state line. For a sufficiently large hydrostatic pressure of the soil skeleton, hardening with plastic volumetric compaction is predicted, which is also denoted as plastic collapse upon wetting, whereas for small values of the hydrostatic pressure of the soil skeleton softening with plastic, volumetric expansion is predicted. The distribution of the volumetric strain with advancing saturation front is shown in Fig. 8. Volumetric compaction, indicated by positive values of the volumetric strain, occurs in the lower central region and in the downstream regions that are characterized by partial saturation. The computed distribution of the vertical displacements of the soil skeleton (negative values denote settlements) is shown in Fig. 9. Shortly after the impoundment settlements in the vicinity of the upstream face of the dam occur, whereas at later stages, uplifting is predicted in this region. By contrast, in the upper central region of the dam, the settlements increase with advancing saturation front.

Conclusions

In embankment dams, both regions of water saturated soil and regions with partial saturation are encountered. During the impoundment of a dam, these regions are changing, and consequently, deformations and changes of the stress state in the dam body occur. For a realistic description of the structural behavior of a dam, a coupled solid-fluid FE-model for partially saturated soils is required, which is characterized by modeling the soil as a three-phase material consisting of a deformable soil skeleton and the fluid phases water and air. An essential part of a coupled model is the constitutive model for partially saturated soil, which allows describing the essential features of partially saturated soil behavior. Here, the well-known Barcelona Basic model, which is formulated in terms of net stress and matric suction, is applied. A computationally efficient return mapping algorithm for this model is proposed, which only requires the solution of a scalar nonlinear equation at the integration point level. The application of the developed coupled FE-model is demonstrated by the coupled transient numerical simulation of the behavior of an embankment dam due to impoundment.

References

[1] Carrere A J, Brunet C, Fry J J, La Barbera A. Bani, Mazza G, Griffths D V, Torres R L, Lane P A. Embankment dams – stability analysis of a homogeneous embankment dam. In: Proceedings of ICOLD Fifth Benchmark Workshop. Denver: ICOLD International Commission of Large Dams, 1999

[2] Alonso E E, Gens A, Josa A. A constitutive model for partially saturated soils. Geotechnique, 1990, 40(3): 405–430

[3] Lewis R W, Schreer B A. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. 2nd ed. New York: John Wiley & Sons, 1998

[4] van Genuchten M Th, Nielsen D R. On describing and predicting the hydraulic properties of unsaturated soils. Annales Geophysicae, 1985, 3(5): 615–628

[5] Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils. New York: John Wiley & Sons, 1993

[6] Hochgürtel T. Numerische Untersuchung zur Beurteilung der Standsicherheit der Ortsbrust beim Einsatz von Druckluft zur Wasserhaltung im schildvorgetrriebenen Tun-nelbau. Dissertation for the Doctoral Degree. Aachen: Rheinisch-Westf alische Technische Hochschule Aachen, 1998

[7] Houlsby G T. The work input to an unsaturated granular material. Geotechnique, 1997, 47(1): 193–196

[8] Borja R I. Cam-clay plasticity. Part V: A mathematical framework for three-phase deformation and strain localization analyses of partially saturated porous media. Computer Methods in Applied Mechanics and Engineering, 2004, 193(48-51): 5301–5338

[9] Schrefler B A. Mechanics and thermodynamics of saturated/unsaturated porous materials and quantitative solutions. Applied Mechanics Reviews, 2002, 55(4): 351–388

[10] Terzaghi K. The shear resistance of saturated soils. In: The 1st International Conference on Soil Mechanics and Foundation Engineering. Cambridge: Harvard University Press, 1936, 1: 54–56

[11] Hickman R J, Gutierrez M. An internally consistent integration method for critical state models. International Journal for Numerical and Analytical Methods in Geomechanics, 2005, 29(3): 227–248

[12] Laigle F, Poulain D, Magnin P. Numerical simulation of La Ganne dam behaviour by a three phases approach. In: Alonso E E, Delage P, eds. Unsaturated Soils/Sol Non Saturés. Proceeding of the First International Conference on Unsaturated Soils, UNSAT 95. Rotterdam: Balkema A. A., 1995, 1101–1108

[13] Gens A, Alonso E E, Delage P. Unsaturated Soil Engineering Practice. In: Geotechnical Special Publication, chapter Computer Modeling and Application to Unsaturated Soils. Washington: ASCE American Society of Civil Engineers, 1997, 68: 299–330

[14] Kohler R. Numerical modelling of partially saturated soils in the context of a three-phase FE-formulation. Dissertation for the Doctoral Degree. Innsbruck: Leopold-Franzens-Universität Innsbruck, 2006

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg