Examples • Verification Problems

Rough Strip Footing on a Cohesive Frictionless Material (FLAC2D)

Problem Statement


The project file for this example is available to be viewed/run in FLAC2D. The data file used for structured grid with quadrialteral zones is shown at the end of this example. The remaining data files can be found in the project.

The prediction of collapse loads under steady plastic-flow conditions can be difficult for a numerical model to simulate accurately (Sloan and Randolph, 1982). As a two-dimensional example of a steady-flow problem, we consider the determination of the bearing capacity of a strip footing on a cohesive frictionless material (Tresca model). The value of the bearing capacity is obtained when steady plastic flow has developed underneath the footing, providing a measure of the ability of the code to model this condition. The model uses plane strain condition, meaning that the the domain and footing are infinitie in out-of-plane direction.

The strip footing is located on an elasto-plastic material with the following properties:

shear modulus (\(G\))

0.1 GPa

bulk modulus (\(K\))

0.2 GPa

cohesion (\(c\))

0.1 MPa

friction angle (\(\phi\))

dilation angle (\(\psi\))

Closed-Form Solution

The bearing capacity obtained as part of the solution to the “Prandtl’s wedge” problem is given by Terzaghi and Peck (1967) as

\[q = (2 + \pi)c\]

in which \(q\) is the average footing pressure at failure, and \(c\) is the cohesion of the material. The corresponding failure mechanism is illustrated in the figure below:


Figure 1: Prandtl mechanism for a strip footing.

FLAC2D Model

For this problem, half-symmetry and plane-strain conditions are assumed in the numerical simulation. The bottom left corner of the domain is located at zero coordinate, the far \(x\)-boundary is at a distance of 20 m from the axis of symmetry, the height of the domain is 10 m, and the area representing the strip footing has a half-width \(a\) = 3 m.

The boundary conditions applied to this domain are sketched in the next figure. A velocity is applied to the model in the negative \(y\)-direction to simulate the footing load. In the data file “PrandtlsWedge.dat”, the rightmost gridpoint of the footing is free in the \(x\)-direction. This condition can be justified because the physical constraint exactly at the edge is ambiguous and can be chosen arbitrarily. Releasing the constraint leads to a more uniform stress distribution under the footing, but does not affect the limit load.


Figure 2: Boundary conditions for FLAC2D analysis — half symmetry.

For the initial model with structured grid, the domain is discretized into 520 zones organized in a graded pattern, as represented in Figure 3, with the grading arranged to increase zone density in the areas of high strain gradient. The area representing the footing overlaps six zones, and a velocity of magnitude 0.5 × 10-5 m/step is applied at the contact nodes for a total of 25,000 calculation steps.

When a velocity is applied to gridpoints to simulate a footing load, the average footing pressure is found by assuming that the footing width is represented by a velocity that varies from the value at the last gridpoint to zero at the next gridpoint. The half-width, \(a\), is then

\[a = A (x_l +x_{l+1})\]

where \(x_l\) is the \(x\)-location of the last applied gridpoint velocity, \(x_{l+1}\) is the \(x\)-location of the gridpoint adjacent to \(x_l\), and \(A\) accounts for the variation. If a linear variation is assumed, then \(A\) = 0.5, and \(a\) = 3.5 m for this problem. The effect of this assumption is discussed in Results And Discussion.

The FISH function load computes the numerical value of the normalized average footing pressure, \(p / c\).


Figure 3: FLAC2D grid view.

Results and Discussion

The load-displacement curve corresponding to the numerical simulation is presented in Figure 4, in which load is the normalized average footing pressure, \(p/c\), and disp is the magnitude of the normalized vertical displacement, \(u_y/a\), at the center of the footing. The numerical value of the bearing capacity, \(q\), is 513.6 kPa, and the relative error is 0.2% when compared to the analytical value of 514.2 kPa.

The apparent width of the footing is taken to be 3 m, plus half the zone width adjacent to the footing edge (because forces are exerted on the footing by this zone, it is assumed that the forces are divided equally between left and right gridpoints).

Note that the error in the bearing capacity is related to the indeterminacy in the apparent width of the footing. The mechanism illustrated in Figure 1 implies a velocity singularity at the ends of the footing. In a numerical simulation, this singularity is spread over the width of one zone. The apparent position of the velocity jump within that zone depends on the exact geometry of the velocity field that develops. In deriving \(q\), it is assumed that the jump occurs half a zone width from the end of the controlled boundary segment (\(A\) = 0.5 in Figure Figure #quad-domain-2d). For finer grids, the indeterminacy in footing width decreases, and the match to the exact solution improves.


Figure 4: Load-displacement curve.

Velocity contours and velocity vectors at the end of the run are presented in Figure 5, showing good agreement with the mechanism in Figure 1.


Figure 5: Velocity field at collapse load.

The same problem was run again using FLAC2D’s nodal mixed discretization (NMD) feature (file “nmd.dat”) for the case of an unstructured grid containing traingular and quadrilateral zones (file “pran-mix-mesh.f2grid”). The NMD calculations are automatically executed for all triangular zones if NMD option is not explicitly specified (however, it can be turned on or off via zone nodal-mixed-discretization command). With NMD, the numerical value of the bearing capacity, \(q\), is 514.5 kPa which is very close to the analytical value of 514.2 kPa. Figure 6 and Figure 7 show the NMD results, and correspond directly with Figure 4 and Figure 5 (results with an all-hex grid).


Figure 6: Load-displacement curve (NMD results).


Figure 7: Velocity field at collapse load (NMD results).


Sloan, S. W., and M. F. Randolph. “Numerical Prediction of Collapse Loads Using Finite Element Methods,” Int. J. Num. & Analy. Methods in Geomech., 6, 47-76 (1982).

Terzaghi, K., and R. B. Peck. Soil Mechanics in Engineering Practice, 2nd Ed. New York: John Wiley and Sons (1967).


Data File


;   2D rough strip footing on Tresca material (Prandtl's wedge problem)
;    -associated plastic flow-
model new
model large-strain off
fish automatic-create off
; Create zones
zone create2d quadrilateral size 6 20 point 1=3.0,0.0 ratio 0.9 0.97
zone create2d quadrilateral size 20 20 point 0=3.0,0.0 ...
                            point 1=20.0,0.0 ratio 1.08 0.97
zone gridpoint initialize position-y 0.5 multiply
; Assign constitutive model and properties
zone cmodel assign mohr-coulomb
zone property bulk 2.e8 shear 1.e8 cohesion 1.e5
zone property friction 0. dilation 0. tension 1.e10
; Assign group name to footing surface
zone face group 'Footing' range position-x 0 3 position-y 10
; Boundary Conditions
zone face apply velocity-normal 0.0 range position-x 0
zone face apply velocity (0.0,0.0) range union position-x 20 position-y 0
zone gridpoint fix velocity (0,-0.5e-5) range position-x 0 3 position-y 10
zone gridpoint free velocity-x range position-x 3 position-y 10

; program call FISH function that monitors footing load
program call 'footing-load'
; Take some histories
fish history load
fish history solution
fish history disp
model history mechanical ratio-local
; Cycle 
model cycle 25000
; Save the model
model save 'pran'