FLAC3D Plastic Hardening Model

How can I get the PH Model for FLAC3D 5.0?
   Triaxial Compression: PH Model versus Mohr-Coulomb Model
   Rough Strip Footing on a Cohesive Frictionless Material
Plastic Hardening Model - Theory and Examples (PDF manual)


Soils being excavated usually exhibit a decrease in stiffness accompanied by irreversible deformation or failure. This plasticity effect accounts for grain crushing and rearrangement, resulting in a nonlinear stress-strain relationship. Numerical modeling of plastic soil hardening is well known inside the civil geotechnical community worldwide. It has become a common constitutive model for civil engineering design in many areas, such as excavations, foundations, tunneling, and soil‐structure interaction in general. Several regulatory agencies, especially in Europe, require (or at least strongly endorse) this type of constitutive model for civil applications.

Itasca has formulated a new Plastic Hardening (PH) constituent model for FLAC and FLAC3D. The PH model is a shear and volumetric hardening constitutive model for the simulation of plastic soil behavior (i.e., when subject to deviatoric loading, they usually exhibit a decrease in stiffness, accompanied by irreversible deformation). The PH model has been developed based on the work by Schanz et al. (1999), which extends the hyperbolic soil linear elastic model (Duncan and Chang, 1970) to an elasto-plastic material model. This avoids the main drawbacks of the hyperbolic formulation:

  • difficulty in detecting and characterizing unloading/reloading, and
  • in specific cases, producing a nonphysical bulk-modulus value that can lead to an erroneous energy generation in the model.

The main features of the new PH model are:

  • hyperbolic stress-strain relationship in axial drained compression;
  • plastic strain in mobilizing friction (shear hardening);
  • plastic strain in primary compression (volumetric hardening);
  • stress-dependent stiffness according to a power law;
  • elastic unloading/reloading compared to virgin loading;
  • memory of pre-consolidation stress;
  • Mohr-Coulomb failure criterion;
  • easy to calibrate using either lab tests or in-situ tests; and
  • it uses familiar properties, names, and conventions in civil engineering.







How can I get the PH Model for FLAC3D 5.0?

The new Plastic Hardening model will be a standard constitutive model in FLAC3D 6.0 when it is released next year. However, Itasca is making the PH model available immediately for free to FLAC3D 5.0 owners who pre-purchase FLAC3D 6.0.

To pre-purchase FLAC3D 6.0, please contact your local Itasca software agent. You must already own a FLAC3D 5.0 license; upgrade options from earlier versions to FLAC3D 5.0 are available.

  • If you have already pre-purchased FLAC3D 6.0 and would like to request the PH Model for FLAC3D 5.0, please contact your local Itasca software agent. The PH model for FLAC3D 5.0 is a C++ User-Defined Model (UDM) and requires the C++ UDM option. If you do not currently own this option, Itasca will provide a free temporary lease license until FLAC3D 6.0 is officially released.


The PH model is very easy and straightforward to calibrate from either lab tests or in-situ tests. Guidelines are provided in the Plastic Hardening Model Manual to streamline the calibration process.The following plots show the very good calibration fits that were achieved using triaxial compression test data (with loading as well as unloading/reloading) of fine Monterey Sand with confining pressures of 0.3, 0.6, and 1.2 MPa and initial void ratios of 0.783, 0.786, and 0.781 (Lade, 1972).


Triaxial Compression: PH Model versus Mohr-Coulomb Model

A comparison between the PH and MC models was done using a single-zone triaxial compression test with 100 kPa confinement stress in FLAC3D.

Stress-strain curves for the two models can be seen in the following plot and it can be verified that:

  • the ultimate failure deviatoric stresses are the same for both models;
  • for the pre-failure curve, the PH and MC models are crossing at the half of the failure stress, which is consistent to the concept of stiffness; and
  • the unloading stiffness is the same as the loading stiffness for the MC model, while these are different for  the PH model

Rough Strip Footing on a Cohesive Frictionless Material

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 (q) is obtained when steady plastic flow has developed underneath the footing, providing a measure of the ability of the software to model this condition.

For this verification problem, half-symmetry and plane-strain conditions are assumed as shown below. The model is 10 m high, 20 m wide, and 1 m thick. The apparent width of the footing is taken to be 3 m, plus half the zone width adjacent to the footing edge (since forces are exerted on the footing by this zone, we assume that the forces are divided equally between left and right grid points). A velocity of magnitude of 0.5 × 10−5 m/step is applied at the contact nodes for a total of 25,000 calculation steps.

The strip footing is located on soil with the following properties:

Table 1. The soil bearing capacity can be calculated as part of the solution to the “Prandtl’s wedge” problem as given by Terzaghi and Peck (1967).

Property Value Material Model
shear modulus (G, GPa) 0.1


Plastic Hardening

bulk modulus (K, GPa) 0.2
cohesion (c, MPa) 0.1
friction angle (φ, °)* 0.01 (~0)
dilation angle (ψ, °) 0
Poisson's Ratio (-) 0.2857
Young's Modulus (GPa) 0.257
E50_ref (GPa) 0.257 Plastic Hardening
Eoed_ref (GPa) 0.2056
pref (MPa) 0.1
m (-) 0.8
Rf (-) 0.9
Knc (-) 0.6
OCR (-) 1, 2, 100

*In order to avoid possible numerical instability, a zero-degree friction angle is replaced by a small value of 0.01 degree, which will have an insignificant impact to the final magnitude of the result.

Displacement contours and velocity vectors at the end of the run are shown below and are in good agreement with the Prandtl mechanism for a strip footing.

The load-displacement curve corresponding to the numerical simulation is shown below, in which p load is the normalized average footing pressure, p/c, and c disp is the magnitude of the normalized vertical displacement, uz/a, at the center of the footing. The bearing capacity predicted by the MC model was close to the Prandtl’s Wedge closed-form solution. However, the pre-failure part of the load-displacement curve by the MC model is based on the assumption that the soil has a constant stiffness before failure, which may not be realistic.  This problem is re-modeled with the PH model. The cohesion, friction, and dilation angles are the same. Three cases (OCR = 1, 2, and 100) are considered for the PH model. A high value of OCR equal to 100 approximately represents no cap hardening effect.

The MC and PH models all predict approximately the same ultimate bearing capacity (with a relative error of 1.7% from the analytical solution). However, from the figure below, the pre-failure parts of the load-displacement curves are quite different. First, all pre-failure parts predicted by the PH model are smooth curves, while the MC model load-displacement curve is initially linear with a sudden change to the smooth curve later on. When the MC curve hits the closed-form line, there is another non-smooth corner. In addition, the lower the OCR in the PH model, the lower the initial tangent stiffness of the load-displacement curve occurs. This case illustrates that although the MC model is sufficient to predict the ultimate bearing capacity, when pre-failure behavior is significant, use of the PH model is recommended.


The benchmark exercise described by Schweiger (2002), consisting of a deep excavation problem in Berlin sand, is the basis for this example application. The problem studied by Schweiger (2002) is also documented in the Plaxis Material Models Manual (2015). Both references are used here. The geometry, basic assumptions, and computational steps adopted for this example are taken from the benchmark exercise. The model parameters are adapted from the Plaxis example description.

A sketch of the problem conditions in the benchmark exercise is shown below. The soil profile consists of two horizontal sand layers. The thickness of the top layer (Layer 1) is 20 m, and the bottom layer (Layer 2) extends to a depth of 100 m below the surface. The excavation is 60 m wide, and the final depth is 16.8 m. Plane-strain conditions and half symmetry are used for the problem. The soil is modeled with the PH model. The wall is modeled by the linear elastic liner structural elements. The anchors are modeled by cable structural elements. The connection between the liner and the cables is through breast-beams with rigid links. Several specifications and the construction sequence, taken from Schweiger (2002) are adopted for this example:
  • The influence of the diaphragm wall construction is neglected.
  • The diaphragm wall is modeled using liner structural elements (Young’s modulus = 30 GPa, Poisson’s ratio = 0.15, thickness = 0.8 m).
  • The horizontal hydraulic cutoff at a 30-m depth is not considered as structural support; the mechanical properties are assumed to be the same as for the surrounding soil.
  • Hydrostatic water pressures, corresponding to water levels, hold inside and outside the excavation. (Full groundwater lowering inside the excavation is performed before the excavation starts.)
  • Anchors are modeled as cables, which are prestressed. The grouted part allows load transfer to the soil.

Stage 1 – Initialize stress state, including groundwater table, 3 m below soil surface.

Stage 2 – Activate diaphragm wall and lower water level to −17.90 m in pit.

Stage 3 – Excavation step 1 (to level −4.80 m).

Stage 4 – Activate anchor row 1 at level −4.30 m and prestress anchors.

Stage 5 – Excavation step 2 (to level −9.30 m).

Stage 6 – Activate anchor row 2 at level −8.80 m and prestress anchors.

Stage 7 – Excavation step 3 (to level −14.35 m).

Stage 8 – Activate anchor row 3 at level −13.85 m and prestress anchors.

Stage 9 – Excavation step 4 (to level −16.80 m).

The simulated wall deflection for both the MC and PH models are shown below as compared with the measured data (Schweiger, 2002). While the MC model predicts several unrealistic observations for wall deflection (ground lifting behind the wall and over-lifting for the excavation base at the final excavation stage), the Plastic Hardening model predicts realistic ground behaviors and compares very well to the measured site data.


PH Model - Theory and Examples Manual




Duncan, J.M. and C.Y. Chang. (1970) “Nonlinear analysis of stress and strain in soil,” J. Soil Mech. Found. Div. 96(5), 1629-1653, October.

Lade, P.V. (1972) The stress-strain and strength characteristics of cohesionless soils. Ph.D. Thesis, University of California, Berkeley, Berkeley, CA.

Plaxis bv, (2015) Material Models Manual, 216 pages.

Rowe, P.W. (1962) “The stress-dilatancy relation for static equilibrium of an assembly of particles in contact,” Proc. Roy. Soc. A. 269(1339), 500-527, October.

Schanz, T., P.A. Vermeer, P.G. Bonnier. (1999) “The hardening soil model: formulation and verification,” in Beyond 2000 in Computational Geotechnics - 10 Years of Plaxis, R.B.J. Brinkgreve, Ed. Rotterdam: Balkema.

Schweiger, H.F. (2002) “Results from numerical benchmark exercises in geotechnics,” in 5th European Conference Numerical Methods in Geotechnical Engineering, Vol. 1, pp. 305–314. Presses de l’ENPC/LCPC, Paris.

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).