COMSOL 6.4 - Slope Stability in an Embankment Dam

Slope stability analysis is an essential technique for predicting the settlement, deformation, and slippage of soil due to various loading and environmental conditions. In an embankment dam, slope stability analysis is important to determine the safety of the dam.

The present model is inspired by an example included in Ref. 1. The pore pressure in the soil is modeled by Darcy’s law, while the Mohr–Coulomb criterion is used for the elastoplastic analysis.

The technique used for studying the slope stability is called the strength reduction method, where the Mohr–Coulomb material parameters are functions of the factor of safety (FOS). In this method, the FOS is gradually increased, which reduces the material parameters and shear strength of the soil, which eventually becomes unstable. This phenomenon can produce a collapse of the embankment at a certain FOS level for given loading conditions. More details of this technique are given in Ref. 1.

Model Definition

Figure 1 shows the cross section of the embankment dam. The dimensions L1, L2, and L3 are 24 m, 5 m, and 24 m, respectively, and the height of the embankment L4 is 12 m. The water level is 10 m and the possible seepage height is 4 m. The total width of the embankment is L1 + L2 + L3. To avoid boundary effects, a soil domain is added below the embankment (not shown in Figure 1).

Figure 1: Illustration of the embankment dam.

In this example, a plane strain approximation is used to model the embankment dam in 2D. The effects of gravity and hydrostatic pressure are also included. The material properties for the Mohr–Coulomb model are parameterized with respect to a factor of safety parameter, FOS. A parametric study increases the FOS parameter, thereby reducing the strength of the soil with every parameter step. The model does not converge for FOS values above 1.92, which indicates a collapse of the slope.

The Mohr–Coulomb yield function and associated plastic potential is

(1)

with

(2)

The parameterized cohesion C and angle of internal friction Φ are given as

(3)

where c is the material cohesion,

and

are angles of internal friction for unsaturated and saturated soils, and p is the pore pressure given by Darcy’s law.

Results and Discussion

The pressure head in the embankment dam is shown in Figure 2; it varies from 0 m to 10 m on the submerged wall, while it is 0 m on the seepage face. Positive pressure head means positive pore pressure, which indicates a saturated soil, while unsaturated soil is represented with zero pressure head. The zero pressure head line in the figure is the location of a phreatic surface that divides saturated from the unsaturated soil.

The elastoplastic analysis does not converge for FOS values greater than 1.92; hence, the simulation is performed until its value becomes 1.92, which is the value at which the slope collapses due to the increase in plastic strains and subsequent reduction in shear strength.

Figure 2: Pressure head in the embankment dam. The zero pressure head line shows the location of phreatic surface.

The equivalent plastic strain just before the collapse shows a different pattern, which gives an indication of the failure mechanism (see Figure 3).

The slip surface is shown in Figure 4. The arrows show the direction of displacement for the soil particles. This figure illustrates the phenomenon of soil slippage. The soil near the lower-right corner does not slip because of the fixed constraint on the lower boundary. The slip surface figure matches well qualitatively with the results given in Ref. 1.

A 3D visualization of the displacement field is shown in Figure 5 with the help of an extrusion dataset. A 2D analysis of an embankment dam is an efficient way to predict soil instability for wide dams where a plane strain approximation is reasonable; a 3D visualization is still possible with the help of the postprocessing tools in COMSOL Multiphysics. This simplification avoids solving a larger numerical problem in 3D.

The plot of the maximum displacement versus the FOS is shown in Figure 6. The maximum displacement increases significantly at around FOS = 1.9, which indicates the onset of the collapse of the slope.

Figure 3: Equivalent plastic strain just before the collapse of the slope.

Figure 4: Slip circle just before the collapse of the slope.

Figure 5: Displacement magnitude in the embankment dam just before the collapse of the slope.

Figure 6: Maximum displacement versus the FOS.

Notes About the COMSOL Implementation

Three stationary studies are created in order to account for the effects of pore pressure and gravity loads on the stability of the slope. In the first study, only Darcy’s Law is computed to solve for the pore pressure profile. In a second study, the initial stresses in the embankment are computed given the hydrostatic, pore pressure, and gravity loads. In the third study step, these initial stresses are taken into account by adding an node. Moreover, the Mohr–Coulomb soil plasticity model is used to analyze the elastoplastic failure of the soil under combined hydraulic and gravity loading.

The angle of internal friction is different for saturated and unsaturated soils, hence an expression that depends on the pore pressure is used. No external stresses are applied in regions of unsaturated soil, since pores are considered interconnected and at constant atmospheric pressure.

The studied problem involves localization of deformation into a band, the slip surface. To maintain mesh objectivity of the solution, a length scale is introduced to the material model to limit strain localization to a predefined width. Here the nonlocal plasticity model is used with a length scale lint = 0.1 m. The value lint is here chosen so that the band of plastic strains is distributed over several mesh elements.

Including plasticity in the material model requires local computations at each Gauss point during the assembly process, which is expensive. By using , the number of Gauss points are reduced by a factor two for the given displacement shape function order and mesh element type. This will speed up the computation significantly.

An additional extrusion dataset is created to generate a 3D plot from the 2D dataset, and the 3D view is adjusted in order to properly visualize it. As stated in Ref. 1, the nonconvergence of the simulation is considered as an indicator of slope failure.

References

1. D.V. Griffiths and P.A. Lane, “Slope Stability Analysis by Finite Elements,” Géotechnique, vol. 49, no. 3, pp. 387–403, 1999.

2. E. Ediger, “Analyzing Slope Stability Through the Shear Strength Reduction Method,” COMSOL Blog, 28 May 2019, www.comsol.com/blogs/analyzing-slope-stability-through-the-shear-strength-reduction-method/.

Geomechanics_Module/Soil/slope_stability

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Geometry 1

Model parameters and interpolation function data are available in the appended text files.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file slope_stability_parameters.txt.

Interpolation 1 (int1)

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file slope_stability_interpolation.txt.

In the text field, type cond.

Locate the section. In the table, enter the following settings:


Argument	Unit
t	m

In the table, enter the following settings:


Function	Unit
cond	m/s

Geometry 1

Construct the 2D geometry using a polygon.

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


x (m)	y (m)
0	0
L1	L4
L1+L2	L4
L1+L2+L3	0
L1+L2+L3*2	0
L1+L2+L3*2	-L4*2
-L1	-L4*2
-L1	0

Add a fillet to remove the geometric discontinuity.

Fillet 1 (fil1)

In the toolbar, click

On the object , select Point 4 only.

In the window for , locate the section.

In the text field, type 5.

Click

Add points at the reservoir level and the possible seepage level to partition the sides of the dam.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Hw*L1/L4.

In the text field, type Hw.

Point 2 (pt2)

Right-click and choose .

In the window for , locate the section.

In the text field, type L1+L2+L3-Hs*L1/L4.

In the text field, type Hs.

Click

Definitions

Variables 1

In the window, expand the > node.

Right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Saturated	dl.Hp>=0		Boolean variable for saturated region
Unsaturated	dl.Hp<0		Boolean variable for unsaturated region
K	cond(dl.Hp)	m/s	Hydraulic conductivity
C	c/FOS	Pa	Parameterized cohesion
PHI	atan(tan(phi_un)/FOS)Unsaturated+atan(tan(phi_sat)/FOS)Saturated	rad	Parameterized friction angle

Maximum 1 (maxop1)

In the toolbar, click

and choose .

Select Domain 1 only.

Add two blank materials, one for the soil and one for the water, then rename them accordingly. For the water material, keep the domain selection empty.

Global Definitions

Soil

In the window, under right-click and choose .

In the window for , type Soil in the text field.

Water

Right-click and choose .

In the window for , type Water in the text field.

Materials

Porous Material 1 (pmat1)

In the window, under right-click and choose > .

In the window for , locate the section.

From the list, choose .

Fluid 1 (pmat1.fluid1)

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat1.solid1)

In the window, right-click and choose .

In the window for , locate the section.

In the θs text field, type 1-psi.

Darcy’s Law (dl)

Porous Matrix 1

In the window, under > > click .

In the window for , locate the section.

From the list, choose .

In the K text field, type K.

Global Definitions

Soil (mat1)

In the window, under > click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho_soil	kg/m³	Basic
Young’s modulus	E	E_soil	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	nu_soil	1	Young’s modulus and Poisson’s ratio

Water (mat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho_wat	kg/m³	Basic

Darcy’s Law (dl)

Add a pressure head to the submerged parts of the downstream and upstream sides of the dam.

Pressure Head 1

In the toolbar, click

and choose .

Select Boundaries 1, 3, and 4 only.

In the window for , locate the section.

In the Hp0 text field, type Hw-y.

Pressure Head 2

In the toolbar, click

and choose .

Select Boundaries 8, 9, and 11 only.

Pressure Head 3

In the toolbar, click

and choose .

Select Boundary 10 only.

In the window for , locate the section.

In the Hp0 text field, type -y.

In the window, click .

In the window for , locate the section.

Select the checkbox.

Gravity 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Use a finer mesh in the region where we expect the slip surface to form.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Size Expression 1

In the window, right-click and choose .

Drag and drop below .

In the window for , locate the section.

In the text field, type if(Y>-L4&&X>L1/2&&X<(L1+L2+L3*1.5),0.75,10).

In the text field, type 50.

Click

Study 1: Darcy’s Law

Disable the default plots for this study.

In the window, click .

In the window for , type Study 1: Darcy's Law in the text field.

Locate the section. Clear the checkbox.

Disable the physics from the study to solve only for the pore pressure.

Step 1: Stationary

In the window, under click .

In the window for , locate the section.

In the column of the table, under , clear the checkbox for .

In the toolbar, click

Add a contour plot of the pressure head.

Results

Pressure Head

In the toolbar, click

In the window for , type Pressure Head in the text field.

Contour 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

Locate the section. From the list, choose .

In the text field, type range(0,3,30).

Locate the section. From the list, choose .

Select the checkbox. In the associated text field, type 0.1.

In the toolbar, click

Solid Mechanics (solid)

Add the water pressure as a boundary load on the downstream side of the dam.

Boundary Load 1

In the toolbar, click

and choose .

Select Boundaries 3 and 4 only.

In the window for , locate the section.

From the list, choose .

In the p text field, type p.

Fixed Constraint 1

In the toolbar, click

and choose .

Select Boundary 2 only.

Roller 1

In the toolbar, click

and choose .

Select Boundaries 1 and 10 only.

Gravity 1

In the toolbar, click

and choose .

Linear Elastic Material 1

Use reduced integration to speed up the simulation.

In the window, click .

In the window for , locate the section.

Select the checkbox.

External Stress 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the pA text field, type p*Saturated.

From the αB list, choose . In the associated text field, type 1.

In the pref text field, type 0.

Add another study to compute the in situ stresses generated by gravity and pore pressure. Disable the default plots for this study.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2: Solid Mechanics (In Situ Stress Initialization)

In the window for , type Study 2: Solid Mechanics (In Situ Stress Initialization) in the text field.

Locate the section. Clear the checkbox.

In the toolbar, click

Solid Mechanics (solid)

Linear Elastic Material 1

In the window, under > click .

Soil Plasticity 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the Ff list, choose .

Use a nonlocal plasticity model to improve the solution during strain localization.

Click to expand the section. From the list, choose .

In the lint text field, type 0.1.

Linear Elastic Material 1

In the window, click .

Initial Stress and Strain 1

In the toolbar, click

and choose .

Add the stresses computed in the second study step as initial stresses. You can access these stresses by using the withsol operator as follows:

In the window for , locate the section.

Specify the S0 matrix as


withsol('sol2',solid.sx)	withsol('sol2',solid.sxy)	withsol('sol2',solid.sxz)
withsol(’sol2’,solid.sxy)	withsol('sol2',solid.sy)	withsol('sol2',solid.syz)
withsol(’sol2’,solid.sxz)	withsol(’sol2’,solid.syz)	withsol('sol2',solid.sz)