Second-Order Work Criterion in the Stability Analysis of an Earth Dam Subjected to Seepage

One of the main problems of geomechanics is the definition of the failure criteria, experimental analysis and simulation of the failure occurrence. The definition of failure is based on the experimentally found fact that a certain group of load paths leads to a limit stress state that cannot be exceeded along these load paths. Therefore, the limit states are asymptotic states.

If we limit ourselves to the class of rate independent elastic–plastic constitutive relationships, they take the general form: (1) dσ=D:dε d{\boldsymbol {\sigma }} = {\bf {D}}:d{\boldsymbol {\varepsilon }} where the elastic–plastic tangent operator D depends on the history of stress and strain expressed by state variables and memory parameters, as well as on the direction of strain increment dɛ (determined by the matrix v=dɛ/‖dɛ‖).

The limit stress state is given by the following relationship: (2) dσ=0 when dε≠0 d{\boldsymbol {\sigma }} = 0\;{\rm{when}}\;\left\| {d{\boldsymbol {\varepsilon }}} \right\| \ne 0 Equation (2) implies: (3) det D=0 and D:dε=0 {\rm{det\;}}{\bf {D}} = 0\;{\rm{and}}\;{\bf {D}}:d{\boldsymbol {\varepsilon }} = 0 If the constitutive equation is incrementally linear, D is independent of v. The first of equations (3) represents the limit state criterion (it is the equation of a hypersurface in six-dimensional stress space). The second equation represents the plastic flow rule, as it determines the direction of dɛ after reaching the limit stress state.

In the case of incrementally linear constitutive relationships (e.g. such as in classical elastic–plastic models, where the tangential operator has a different form under loading and unloading), where D discontinuously depends on v, det D=0 is also an equation of hypersurface in stress space involving several plastic mechanisms. In this case, it is a peculiar generalised plastic flow rule with vertices represented locally by pyramids [8].

Such a perception of failure in an elastic–plastic material resulted in a widely understood analysis of the limit stress state, in which it is assumed that the limit condition is reached at every point of the area under consideration. In reality, however, there are various forms of failure: localised in shear bands ([2], [15], [19]), localised in compression or dilatancy bands, a diffuse form associated with geometric instability (buckling, torsional buckling) or a field of chaotic displacements. In the case of materials with the non-associated flow rule (with an asymmetric elastic–plastic matrix D), the theory of plasticity shows that by applying the appropriate bifurcation criteria, localised or diffuse forms of failure may occur before the plastic limit condition is reached [1].

According to the different forms of failure, different bifurcation criteria are given in the literature. Concerning the formation of the shear band as a result of the localisation of plastic deformations, the Rice criterion [15], based on the description of such an emerging shear band with normal n, corresponds to the vanishing value of the determinant of the acoustic tensor: (4) det nTDn=0 {\rm{det\;}}\left( {{{\bf {n}}^T}{\bf {Dn}}} \right) = 0 For a material with the non-associated flow rule, condition (4) may occur before the limit state criterion (3) is met. This was verified experimentally for dense sand [4]. The strain localisation corresponds to the bifurcation of the strain field from a diffuse to strictly discontinuous form. This type of bifurcation can be called ‘discontinuous bifurcation’.

Continuous bifurcations [4] are also a form of failure without strain localisation. This form is called ‘diffuse failure’, a response path that bifurcates with the loss of constitutive uniqueness at the bifurcation point. Depending on the type of load path control (stress, strain or mixed), different response paths from the bifurcation point are possible. For certain types of control, the load ceases to be controllable in the sense of Nova [12], [13].

Continuous bifurcations and the associated diffuse forms of failure can be recognised by Hill's stability criterion [8], which is related to the nil second-order work for the unstable states, that is to the zero value of the determinant of the symmetric part of the matrix D in incrementally linear constitutive relationships.

This article presents the results of stability calculations of an earth dam subjected to seepage through its body due to changing upstream water levels. The Hardening Soil model was used as the constitutive model in the analysis. The safety factor was calculated using the shear strength reduction method for different upstream water levels using the finite element software ZSoil. Calculated increments of stresses and strains during changes in the water levels were used to calculate the values of the normalised second-order work at Gauss points, that is local Hill criterion. Maps of these values in the body of the earth dam and subsoil were made using a script in Python. The decreasing safety factor corresponds to the decreasing second-order work during water rising until the appearance of a zone with negative values of the second-order work. Negative second-order work values appear at the foot of the downstream slope for the highest upstream water level and the lowest factor of safety.

Finite element stability analyses using the second-order work criterion are very rare in the literature ([7], [9], [10], [14]). The results of the calculations presented in this article complement the achievements in this area.

The basic definition of stability was given by Lyapunov [11]. In continuum mechanics, it is as follows: a stress–strain state in a material with a defined load history is said to be stable when, for any positive scalar ɛ, there is a positive number η(ɛ) such that for all load increments limited by η, the resulting responses remain bounded by ɛ.

According to this definition, all limit stress states are unstable. If we consider a ‘small’ increase in stress directed from the limit state, it will cause a ‘small’ response in the form of an increase in strain, often treated as elastic. On the other hand, a ‘small’ increase in stress directed beyond the limit state causes a ‘large’ deformation response. Lyapunov's definition clearly indicates that material instability can occur in an elastic–plastic medium. Since some limit states occur inside the domain bounded by the Mohr–Coulomb limit condition (and this is related to the non-associated flow rule), some instabilities are expected to occur before the Mohr–Coulomb limit state is reached.

Nevertheless, the Lyapunov definition is unsuitable for studying instability in geomaterials. For this purpose, a sufficient criterion of Hill's stability is used [8]. According to Hill, a stress–strain state is unstable if one load direction can be continued in the sense of an infinitesimal increment without the supply of energy from an external source. In practice, this means that the deformation process will continue spontaneously. And in fact, it is; usually, failure requires some energy from an external source (e.g. foundation load), and in some cases, no energy is required (landslides, rock falls).

Hill's stability condition states that a stress–strain state is stable if for all increments (dσ, dɛ) related by the constitutive law the value of the second-order work is positive: (5) d2W=dσ:dε>0 ∀dε≠0 {d^2}W = d{\boldsymbol {\sigma }}:d{\boldsymbol {\varepsilon }} > 0\;\;\;\forall \left\| {d{\boldsymbol {\varepsilon }}} \right\| \ne 0

For Drucker's hyperelastic–plastic constitutive hypotheses (dσ:dɛ^p>0) implies the Hill condition since in hyperelasticity dσ:dɛ^e>0 and dσ:dɛ=dσ:dɛ^e+dσ:dɛ^p, due to the decomposition of the strain increment into the elastic, dɛ^e, and plastic part, dɛ^p, made in the classical theory of plasticity.

The second-order work (5) for a material with non-associated flow rule depends only on the symmetric part of the constitutive operator D: (6) d2W=dσ:dε=dε:D:dε=dε:Dsym:dεDsym=1/2Dijkl+Dklij \matrix{{{d^2}W = d{\boldsymbol {\sigma }}:d{\boldsymbol {\varepsilon }} = d{\boldsymbol {\varepsilon }}:{\bf {D}}:d{\boldsymbol {\varepsilon }} = d{\boldsymbol {\varepsilon }}:{{\bf {D}}^{sym}}:d{\boldsymbol {\varepsilon }}} \cr {{{\bf {D}}^{sym}} = 1/2\left( {{D_{ijkl}} + {D_{klij}}} \right)} \cr } Thus, for incrementally linear constitutive relations, condition (5) is equivalent to (7) det Dsym>0 {\rm{det\;}}{{\bf {D}}^{sym}} > 0 In the axial symmetry and plane strain conditions, a unique relationship was found between the second-order work sign d²W and the sign of the determinant of the symmetric part of the constitutive matrix D: the area where d²W≤0 is contained in the area where detD^sym≤0, for the general case of incrementally non-linear constitutive relationships [7].

If conditions (5) and (7) are met, stability is guaranteed. In the associated plasticity, stable behaviour implies a positive hardening modulus. The loss of stability occurs at yielding when the stress meets the plastic limit condition (the stress remains on the limit surface).

The loss of stability in materials with the non-associated flow rule is related to material properties that cause certain reactions to occur within the failure condition, such as static liquefaction [3].

Simulations of material instability using an incrementally linear and incrementally non-linear model are carried out by, among others, Darve with co-workers ([4,5,6,7], [9]). They showed that there is an area of instability located strictly inside the boundary surface. Its existence was confirmed by simulation of element tests of loose sand [6] and dense sand [4] in the axi-symmetric conditions, as well as in a plane strain condition and for radial stress paths in the deviatoric plane for loose and dense sand [9].

The area of instability in the stress space contains points of potentially unstable states, that is such states of stress and strain related by the constitutive law, in which there is at least one direction of stress or strain increment leading to a zero or negative second-order work value. After the loss of stability, the unstable state will be maintained if the load direction coincides with the unstable direction and if the load is constant. In other words, if the load change occurs under stress (or force) control, the deformation of the material is not kinematically constrained. In such a case, the deformation response will be arbitrarily large and, therefore, will not satisfy the Lyapunov stability criterion.

Determination of the domain of instability is done by simulating the element test (triaxial compression or tension, radial stress paths in the deviatoric plane, load in a plane strain condition and others) to the assumed stress values. Along the stress path, a check is made for the existence of a load direction resulting in zero second-order work. When such a direction exists, a given point on the stress path is part of the boundary of the region of instability. Fig. 1 shows a scheme for determining the domain of instability in the principal stress space.

Figure 1:

Domains of instability were found by the results of simulations of element tests of both loose and dense sands, in both drained and undrained conditions [7]. Dense sand exhibits a smaller domain of instability than loose sand. The domains of material instability identified in the simulations of triaxial compression and tension of dense sand described by the incrementally non-linear Darve model are shown in Fig. 2 as shaded areas before reaching the Mohr–Coulomb limit condition [7].

Figure 2:

The domain of instability was found in the analysis of changes in the homogeneous stress and strain fields, which occur in the element tests. In the case of solving the boundary problem, failure (material instability) will occur only when the direction of the stress path is appropriate and the boundary conditions allow the formation of the failure mechanism. Predicting Hill instability in engineering problems is possible, which was previously demonstrated by, among others, Darve and Laouafa [7], Sternik [17], [18], Prunier et al. [14].

Hill proposed a criterion for the stability of an elastic–plastic body, which, assuming small deformations, takes the form: (8) ∫Vdσ:dεdV>0 ∀dε≠0 \int_V {d{\boldsymbol {\sigma }}:d{\boldsymbol {\varepsilon }}dV > 0\;\;\;\forall \left\| {d{\boldsymbol {\varepsilon }}} \right\| \ne 0} where dσ is the Cauchy stress increment, and dɛ is the strain increment associated with the stress increment by the constitutive law. This criterion is sufficient but not necessary to maintain stability. This means that as long as condition (8) is met, we are sure that the body's stability is maintained, but when condition (8) does not hold, there is a possibility of unstable behaviour. Condition (8) at the material point has the form (5).

In finite element calculations, the value of second-order work is calculated at all points of Gaussian integration (pi) of the discretised area as the product of the increments of the stress and strain increments between two successive equilibrium states of the analysed problem. A convenient form of presentation that increases the readability of the distribution of the second-order work values is to normalise the d₂W value by the product of stress and strain norms: (9) d2Wpinorm=d2Wpidσpidεpi=dσpi:dεpidσpidεpi {d^2}W_{pi}^{norm} = {{{d^2}{W_{pi}}} \over {\left\| {d{{\boldsymbol {\sigma }}_{pi}}} \right\|\left\| {d{{\boldsymbol {\varepsilon }}_{pi}}} \right\|}} = {{d{{\boldsymbol {\sigma }}_{pi}}:d{{\boldsymbol {\varepsilon }}_{pi}}} \over {\left\| {d{{\boldsymbol {\sigma }}_{pi}}} \right\|\left\| {d{{\boldsymbol {\varepsilon }}_{pi}}} \right\|}} In the appropriate conjugate spaces, the value (9) corresponds to the cosine of the angle between the stress and strain vector and is in the range <−1;1> [10].

Hill's criterion in the domain (the global second-order work) given by equation (8) in a finite element model is the sum of the weighted d²W_pi values at all integration points N_pi of the domain: (10) D2W=Σpi=1Npidσpi:dεpiωpiJpi=Σpi=1NpiD2Wpi {D^2}W = \mathop \sum \nolimits_{pi = 1}^{{N_{pi}}} d{{\boldsymbol {\sigma }}_{pi}}:d{{\boldsymbol {\varepsilon }}_{pi}}{\omega _{pi}}\left| {{J_{pi}}} \right| = \mathop \sum \nolimits_{pi = 1}^{{N_{pi}}} {D^2}{W_{pi}} where the product of the Jacobian at the point of integration |J_pi| and the weight ω_pi expresses the part of the surface area of the element assigned to the given point of integration.

The calculations were carried out in several stages, in which the water level in the subgrade was first raised and then lowered, separated by checking the stability of the earth dam slopes using the shear strength reduction method. Changes in boundary conditions are shown in Fig. 6. The coupled hydro-mechanical analysis was carried out in a way that ensured free flow of water in the soil mass. The phreatic surface moved up steadily with the increase in the value of hydraulic heads on the vertical edges of the model (Fig. 6). The rate of changes in the hydraulic boundary conditions enabled the stabilisation of the phreatic surface at every stage of the calculations. Changes in the phreatic surface location at selected moments of time t = 0, 10, 18 days are shown in Fig. 7.

Figure 6:

Figure 7:

As previous work by Darve and co-workers has shown material instability may occur in both loose and dense sands in undrained and drained conditions. In order to be able to predict it, constitutive relationships should take into account hardening and non-associated plasticity.

Instability can be detected using Hill's criterion, which can be applied to a point (local criterion) or a domain (global criterion). This is the first criterion of bifurcation understood as a sudden or discontinuous change in the material response while maintaining constant or slightly disturbed load parameters.

The article presents the FEM calculations of the ‘earth dam – subgrade’ system, in which changes in stress and strain appear due to upstream water pressure on the earth dam slope and seepage through the earth dam body. The stress–strain relationship of the Hardening Soil model was used in the calculations. The increasing pore water pressure due to the upstream water level rise leads to a zone of negative second-order work values at the downstream slope's toe. This is the area where slope instability usually begins. Thus, the local Hill criterion indicates instability that corresponds to reality. Negative second-order work also means that in the analysed case, the material instability occurs before reaching the limit stress condition for which d²W = 0. It is to note that in the case considered, points of instability were detected for dense sand in drained conditions.

The decreasing values of the second-order work in the entire area of the modelled problem of rising upstream water level correspond to the decreasing value of the earth dam safety factor determined by a more classical shear strength reduction method.