Soil Reinforcement in Finite Element Analysis
The integration of soil reinforcement elements such as geotextiles, soil nails, and ground anchors into finite element slope stability analysis represents a significant advancement in modeling stabilized slopes. These reinforcement systems fundamentally alter the stress distribution and failure mechanisms within slopes, requiring sophisticated modeling approaches to capture their beneficial effects accurately (Duncan & Wright, 2005).
The modeling of reinforced slopes presents unique challenges because the reinforcement elements typically have dramatically different mechanical properties compared to the surrounding soil. Reinforcement elements are usually much stiffer in tension and often have negligible compressive strength, creating a highly anisotropic composite material that requires specialized finite element formulations.
Truss Element Approach
While there are numerous ways to simulate soil reinforcement in the finite element method including the equivalent force method and interface element modeling, the most straightforward method involves representing reinforcement elements as one-dimensional truss elements embedded within the two-dimensional soil continuum. These truss elements are characterized by their axial stiffness \(EA/L\), where \(E\) is the elastic modulus of the reinforcement material, \(A\) is the cross-sectional area, and \(L\) is the element length.
This approach is particularly effective for modeling geosynthetic reinforcement, soil nails, and tie-back anchors. The truss elements can only carry tension loads up to a specified tensile strength limit \(T_{max}\), beyond which they either yield plastically to a residual strength \(T_{res}\), or fail completely. The inability to carry compression loads accurately reflects the behavior of flexible reinforcement materials like geotextiles and ensures that the reinforcement cannot resist compressive buckling. The truss elements are oriented along the centerline of the physical reinforcement and connected to the surrounding soil elements through shared nodes. This connection ensures that the reinforcement participates in the overall deformation pattern of the slope while contributing its tensile resistance to improve stability.
Truss elements are incorporated into the XSLOPE finite element mesh by passing the geometry of the reinforcement lines from the input template to the mesh generation process. The reinforcement lines are discretized into multiple truss elements based on the specified mesh density (target_size). The 1D elements are fully integrated with the 2D elements - each 1D element corresponds to the edge of two adjacent 2D elements and both the 1D and 2D elements share the same nodes. The 1D elements have their own set of material properties corresponding to the properties of the corresponding reinforcement lines input by the user and include \(T_{max}\), \(T_{res}\), \(E\), and cross-sectional area \(A\).
Only the two end nodes of each 1D element are used for computing truss stiffness and axial forces. When the mesh uses quadratic (8-node) 2D elements, the 1D elements may have a mid-side node shared with the adjacent 2D elements. This mid-node participates in the 2D element shape functions for displacement interpolation but is ignored for the truss element formulation. A 2-node linear truss element is exact for a prismatic bar with constant axial stiffness, so the quadratic interpolation adds no physical fidelity for the 1D element.
The meshing algorithms used in XSLOPE, including the integration of 1D and 2D elements for problems involving soil reinforcement are documented in the Automated Mesh Generation page.
Mathematical Formulation
Truss Element Stiffness Matrix: Each 1D truss element contributes to the global stiffness matrix through its element stiffness matrix. For a truss element with nodes \(i\) and \(j\), the element stiffness matrix in local coordinates is:
\([K_e]_{local} = \dfrac{AE}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\)
where \(A\) is the cross-sectional area, \(E\) is the elastic modulus, and \(L\) is the element length.
Coordinate Transformation: The local stiffness matrix must be transformed to global coordinates using the transformation matrix \([T]\):
\([K_e]_{global} = [T]^T [K_e]_{local} [T]\)
For a truss element oriented at angle \(\theta\) to the horizontal:
\([T] = \begin{bmatrix} \cos\theta & \sin\theta & 0 & 0 \\ 0 & 0 & \cos\theta & \sin\theta \end{bmatrix}\)
Assembly Process: The global stiffness matrix combines contributions from both 2D soil elements and 1D truss elements:
\([K]_{global} = \sum_{soil} [K_e]_{soil} + \sum_{truss} [K_e]_{truss}\)
Force Vector Assembly: The global force vector includes both soil body forces and any applied forces on reinforcement:
\(\{F\}_{global} = \{F\}_{soil} + \{F\}_{reinforcement}\)
Force Behavior and Failure Modes
The forces and failure modes in 1D truss elements are analyzed in an iterative fashion. In general, each truss element has a maximum allowable tensile capacity \(T_{allow}\) and a residual tensile capacity \(T_{res}\) derived from the user-specified reinforcement parameter inputs. The axial force in each element is calculated as \(T = (AE/L) \cdot \delta\), where \(\delta\) is the element elongation (the component of relative nodal displacement along the element axis). When the tensile force meets or exceeds \(T_{allow}\), the element's effective capacity drops to \(T_{res}\). These two parameters can be adapted to simulate a variety of behaviors and failure modes, depending on the reinforcement material and loading conditions:
Complete Failure Model:
A complete failure model is appropriate for brittle materials (some steel cables, fiber reinforcement).
- When \(T > T_{allow}\), the element's effective force is capped at zero via body-force corrections
- This is simulated by setting \(T_{res} = 0\) in the input for the reinforcement line
Peak-Residual Model:
This model is appropriate for ductile materials (geotextiles, some geosynthetics) and is the most common scenario.
- When \(T > T_{allow}\), the element's effective force is capped at \(T_{res}\) via body-force corrections
- The element remains in the stiffness matrix at full elastic stiffness, but corrective forces prevent it from carrying more than \(T_{res}\)
- Typical residual strength ratios for geosynthetics: \(T_{residual}/T_{allow} = 0.3-0.7\)
Pullout Failure Model:
For each reinforcement line, it is assume that the tension force in the reinforcement is zero at the two ends and increases linearly with distance along the line as frictional resistance between the reinforcement and the surrounding soil developed. Full tension force develops over a pullout distance, \(L_p\).
- Pullout failure may occur in elements where the embedment length is less than the pullout length \(L_p\)
- The available strength is limited by pullout resistance rather than material strength
- For elements at distance \(d\) from the reinforcement end where \(d < L_p\):
\(T_{available} = T_{allow} \times \frac{d}{L_p}\)
- Pullout failure is typically sudden and complete (no residual capacity)
- Progressive pullout can occur as elements near the ends fail sequentially
Tension-Only Behavior:
Truss elements are restricted to carry only tension forces. This is implemented through body-force corrections within the viscoplastic iteration loop:
- After each iteration, the axial force in each element is computed from the current displacement field
- If compression develops (\(T < 0\)), a corrective body force is applied that cancels the compressive force
- The element remains in the stiffness matrix at full elastic stiffness; the correction enters through the load vector
- This approach is consistent with the viscoplastic initial-stiffness method used for soil elements, where the stiffness matrix is factored once and all nonlinearity is driven through load corrections
Integration with Viscoplastic Iteration
The 1D truss element nonlinearity (tension-only behavior, capacity limits, and failure) is handled through body-force corrections within the Griffiths & Lane (1999) viscoplastic iteration loop, using the same initial-stiffness approach that governs the 2D soil elements. The key principle is that the global stiffness matrix \([K]\) is assembled once with the full elastic stiffness of all elements (both 2D soil and 1D truss) and pre-factored for efficient repeated solves. All nonlinear behavior is then driven entirely through corrections to the right-hand-side load vector.
Assembly: The global stiffness matrix includes contributions from both 2D soil elements and 1D truss elements:
\([K]_{global} = \sum_{soil} [K_e]_{soil} + \sum_{truss} [K_e]_{truss}\)
This matrix is factored once (via sparse LU decomposition) and reused for all viscoplastic iterations.
Iteration procedure: At each viscoplastic iteration, after solving for the updated displacement field \(\{u\}\):
-
For each 1D truss element, compute the axial force from the current displacements:
\(\delta = (u_{x,j} - u_{x,i})\cos\theta + (u_{y,j} - u_{y,i})\sin\theta\) \(T = \dfrac{AE}{L} \cdot \delta\)
-
Determine whether a correction is needed:
- If \(T < 0\) (compression): set \(\Delta T = -T\) (cancel the compressive force entirely)
- If \(T > T_{allow}\) and the element has not previously failed: mark the element as failed and set \(\Delta T = T_{res} - T\)
- If \(T > T_{res}\) and the element has previously failed: set \(\Delta T = T_{res} - T\)
- Otherwise: no correction (\(\Delta T = 0\))
-
Convert the axial force correction to equivalent nodal forces and add to the load vector:
\(\{F\}_{correction} = \Delta T \begin{Bmatrix} -\cos\theta \\ -\sin\theta \\ \cos\theta \\ \sin\theta \end{Bmatrix}\)
These corrections are added to the same load vector that receives the soil viscoplastic strain corrections (\([B]^T [D] \{\varepsilon_{vp}\}\)). The factored stiffness matrix then solves the corrected system, and the process repeats until convergence.
Failure irreversibility: Once an element is marked as failed (having exceeded \(T_{allow}\)), it remains failed for all subsequent iterations within that analysis. Its effective capacity permanently drops from \(T_{allow}\) to \(T_{res}\). This models the irreversible nature of material yielding or pullout failure.
Strength Reduction and Reinforcement
In the Shear Strength Reduction Method (SSRM), only the soil shear strength parameters \(c\) and \(\tan\phi\) are reduced by the strength reduction factor \(F\). The reinforcement properties (\(T_{max}\), \(T_{res}\), \(E\), \(A\)) are held constant throughout the SSRM bisection. The resulting factor of safety represents the margin of safety in the soil strength, given the structural reinforcement as-designed. This follows standard practice (Duncan & Wright, 2005) where reinforcement capacity is treated as a structural property independent of the soil strength reduction.
This also means that the truss element contributions to the global stiffness matrix do not change between SSRM bisection steps — only the soil yield parameters change — so the same pre-factored stiffness matrix approach remains efficient.
Reinforcement Line Input Parameters and Element Properties
In the Excel input template used by XSLOPE, the user can define up to 20 reinforcement lines by entering the reinforcement line geometry and properties into the lines of a table. Each row of the table includes the following:
| Item | Description |
|---|---|
| x1, y1 | The x and y coordinates of the left end of the line |
| x2, y2 | The x and y coordinates of the right end of the line |
| Tmax | Maximum allowable tensile force |
| Tres | Residual tensile force |
| Lp1 | The pullout length on the left side |
| Lp2 | The pullout length on the right side |
| E | The modulus of elasticity of reinforcement material |
| Area | The cross-sectional area of the reinforcement material |
The units for E and Area need to be compatible with each other and with the other weight and length units used. For metric units, E should be in \(kPa\) and Area should be in \(m^2\). For English units, E should be in \(psf\) and Area should be in \(ft^2\). Alternately, E could be in \(psi\) as long as Area is in \(in^2\).
Element Discretization and Capacity Assignment
A separate pullout length (Lp) is used for each end since each end may be embedded in a separate soil with different shear resistance values.
During mesh generation, each reinforcement line is discretized into multiple truss elements based on the specified mesh density. The discretization process follows these steps:
- Material Property Assignment: Each truss element along the line receives the same material properties:
Cross-sectional area: \(Area\)
Elastic modulus: \(E\)
Element stiffness: \(K_e = AE/L\) (where L varies based on element length)
- Tensile Capacity Assignment: Each truss element is assigned an allowable \(T_{allow}\) and residual \(T_{res}\) tensile capacity based on the following logic:
For an element whose center is at distance \(d\) from the nearest reinforcement end, where \(L_p\) is the pullout length corresponding to that end (\(L_{p1}\) if nearest to end 1, \(L_{p2}\) if nearest to end 2):
\(T_{allow} = \begin{cases} T_{max} \cdot \dfrac{d}{L_p} & \text{if } d < L_p \\ T_{max} & \text{if } d \geq L_p \end{cases}\)
\(T_{res} = \begin{cases} 0 & \text{if } d < L_p \\ T_{residual} & \text{if } d \geq L_p \end{cases}\)
This approach ensures that elements near the reinforcement ends have reduced capacity (starting from zero at the ends), while elements beyond the pullout length carry the full design strength. The linear variation within the pullout length reflects the gradual development of pullout resistance through interface friction. Since each end of a reinforcement line may be embedded in a different soil with different shear resistance, the appropriate pullout length (\(L_{p1}\) or \(L_{p2}\)) is selected based on which end is nearest to the element centroid. For the residual capacity, if \(d < L_p\), the residual capacity is set to zero because it is assumed that a pullout failure is sudden and complete and there is no residual capacity. If \(d \geq L_p\), the residual strength specified by the user for the element is used.
Axial Stiffness (EA)
The analysis depends only on the product \(EA\) (the axial stiffness, sometimes called the tensile stiffness or \(J\) in geosynthetic specifications), not on \(E\) and \(A\) independently. Any combination of \(E\) and \(A\) that produces the same \(EA\) will give identical results. The axial stiffness controls how much the reinforcement must elongate before mobilizing its tensile capacity:
\(EA = \dfrac{T_{max}}{\varepsilon_{rupture}}\)
where \(\varepsilon_{rupture}\) is the strain at which the reinforcement reaches its ultimate tensile strength.
Determining Reinforcement Line Pullout Lengths
The pullout length \(L_p\) represents the distance from each end of the reinforcement over which the full tensile strength is mobilized. This variation captures the physical reality that pullout resistance must develop over a finite distance from the reinforcement ends through interface friction between the reinforcement and surrounding soil. This friction cannot be mobilized instantaneously but requires relative displacement to develop, creating the gradual strength buildup characteristic of all reinforcement systems. Pullout length can be estimated as follows:
For Soil Nails:
\(L_p = \dfrac{T_{max}}{\alpha \pi D \sigma_n' \tan \phi_{interface}}\)
where:
\(T_{max}\) = design tensile capacity of the nail
\(\alpha\) = surface roughness factor (0.5-1.0 for grouted nails)
\(D\) = effective nail diameter
\(\sigma_n'\) = average effective normal stress along the nail
\(\phi_{interface}\) = interface friction angle (typically 0.8-1.0 times soil friction angle)
For Geotextiles:
\(L_p = \dfrac{T_{max}}{2 \alpha \sigma_n' \tan \phi_{interface}}\)
where the factor of 2 accounts for friction on both sides of the geotextile.
These equations are a general guide that can be used to come up with reasonable estimates of Lp. Typical values are as follows:
| Reinforcement Type | Pullout Length \(L_p\) (m) | Notes |
|---|---|---|
| Soil Nails | 1.5 - 3.0 | Depends on soil conditions and nail diameter |
| Geotextiles | 0.5 - 1.5 | Depends on normal stress and surface texture |
| Geogrid | 1.0 - 2.0 | Depends on aperture size and bearing resistance |
Wished-in-Place Analysis and EA Selection
XSLOPE currently uses a wished-in-place approach: the entire slope (all soil layers and all reinforcement) is assumed to exist in its final geometry, and gravity is applied in a single step. The reinforcement starts at zero strain and zero force. Tension develops only through deformation that occurs during the gravity application and subsequent SSRM strength reduction. This differs from reality, where reinforcement accumulates tension progressively as each soil lift is placed during construction.
The wished-in-place approach is conservative — it underestimates the reinforcement contribution because it misses the construction-induced pre-tension. However, for the reinforcement to provide meaningful stabilization in a wished-in-place SSRM analysis, it must be stiff enough to develop significant force from the relatively small displacements that occur near the incipient failure state. Low-stiffness reinforcement may undergo insufficient strain during SSRM to mobilize its capacity, leading to an unrealistically low factor of safety.
Parametric studies show that the computed factor of safety increases with \(EA\) and then plateaus above a threshold stiffness, beyond which further increases in \(EA\) have negligible effect. For typical reinforced slope geometries, this plateau is reached at approximately \(EA \approx 100\)–\(200 \times T_{max}\). Below approximately \(EA \approx 50 \times T_{max}\), the reinforcement may not mobilize enough force to significantly improve the factor of safety in a wished-in-place analysis.
The following table provides recommended \(EA\) values for wished-in-place SSRM analysis. These values are deliberately at the stiffer end of the physical range for each material type, because the wished-in-place approach requires sufficient stiffness to compensate for the absence of construction-induced pre-tension:
| Reinforcement Type | Recommended \(EA/T_{max}\) | \(\varepsilon_{rupture}\) | Notes |
|---|---|---|---|
| Woven geotextiles | \(50\)–\(100\) | 1–2% | Use stiffer end for SSRM |
| HDPE geogrids | \(50\)–\(100\) | 1–2% | Uniaxial, reinforcement grade |
| PET geogrids | \(100\)–\(200\) | 0.5–1% | Higher stiffness than HDPE at same strength |
| Steel strips | \(500\)–\(2{,}000\) | 0.05–0.2% | Very stiff, minimal elongation |
| Soil nails (grouted) | \(1{,}000\)–\(5{,}000\) | 0.02–0.1% | Based on steel bar + grout composite |
Typical E and Area Values
The following tables provide representative \(E\) and \(Area\) values for common reinforcement materials. These are intended as starting points when manufacturer-specific data is not available. Any combination of \(E\) and \(Area\) producing the target \(EA\) will give identical results.
English Units:
| Material | \(E\) (psf) | \(Area\) (ft\(^2\)/ft) | \(EA\) (lb/ft) | \(T_{max}\) (lb/ft) |
|---|---|---|---|---|
| Woven geotextile (light) | 500,000 | 0.01 | 5,000 | 100 |
| Woven geotextile (heavy) | 2,000,000 | 0.02 | 40,000 | 500 |
| HDPE geogrid | 2,000,000 | 0.03 | 60,000 | 800 |
| PET geogrid | 5,000,000 | 0.02 | 100,000 | 1,000 |
| Steel strip (galvanized) | 400,000,000 | 0.0003 | 120,000 | 5,000 |
| Soil nail (grouted, #8 bar) | 600,000,000 | 0.0006 | 360,000 | 15,000 |
Metric Units:
| Material | \(E\) (kPa) | \(Area\) (m\(^2\)/m) | \(EA\) (kN/m) | \(T_{max}\) (kN/m) |
|---|---|---|---|---|
| Woven geotextile (light) | 25,000 | 0.003 | 75 | 1.5 |
| Woven geotextile (heavy) | 100,000 | 0.006 | 600 | 7 |
| HDPE geogrid | 100,000 | 0.009 | 900 | 12 |
| PET geogrid | 250,000 | 0.006 | 1,500 | 15 |
| Steel strip (galvanized) | 20,000,000 | 0.00009 | 1,800 | 75 |
| Soil nail (grouted, #8 bar) | 30,000,000 | 0.0002 | 6,000 | 220 |
Staged Construction Alternative
In practice, reinforced slopes and walls are built in lifts. Each soil layer is placed and compacted on top of previously installed reinforcement, which develops tension in the reinforcement before the next lift is added. By the time the slope reaches its final geometry, the lower reinforcement layers may have accumulated significant construction-induced pre-tension.
A staged construction analysis models this process by activating soil layers and reinforcement elements sequentially in the FEM, solving for equilibrium at each stage. The stress state — including locked-in reinforcement tension — carries forward from each stage to the next. The SSRM analysis then begins from the end-of-construction stress state rather than from a zero-strain condition.
Staged construction analysis is supported by all major commercial geotechnical FEM packages (PLAXIS, FLAC, RS2/Phase2, SIGMA/W) and is the recommended approach when:
- The reinforcement has low axial stiffness (extensible geosynthetics)
- The slope is tall with many reinforcement layers
- Accurate prediction of reinforcement forces is important (not just the factor of safety)
- The analysis needs to match instrumented field measurements
For the wished-in-place approach currently used by XSLOPE, selecting \(EA\) values at the stiffer end of the recommended range (see table above) partially compensates for the absence of construction-induced pre-tension and produces conservative but reasonable factors of safety. Staged construction analysis may be implemented in a future version of XSLOPE.
References
Duncan, J.M., & Wright, S.G. (2005). Soil Strength and Slope Stability. John Wiley & Sons.
Griffiths, D.V., & Lane, P.A. (1999). Slope stability analysis by finite elements. Geotechnique, 49(3), 387-403.
Smith, I.M., & Griffiths, D.V. (2004). Programming the Finite Element Method (4th ed.). John Wiley & Sons.