Jump to navigation Jump to search



In this section, the basic concepts of the mechanics of solids are briefly summarized. For later reference, we consider a solid with volume and boundary subjected to body forces . On part of the boundary, , the displacements are prescribed while on the other part, , the tractions, , are prescribed.

Figure 1: Solid of volume with boundary subjected to tractions on and supported on .

Standard continuum mechanics sign conventions are adopted throughout, i.e. tensile stresses and strains are positive.

1.1 Stress and equilibrium

In the general three-dimensional case, the state of stress at a point is defined in terms of the six-dimensional stress vector:

In plane strain, two of the shear stresses are zero. Let be the out-of-plane direction. We then have:

For a two-dimensional plane strain body, the equilibrium equations are given by

where are body forces. These equations may be expressed in matrix form as


The static boundary conditions may be written as

where and are the components of the traction vector and is the outward normal to the boundary. These equations may be expressed in matrix form as


1.2 Displacement, strain and compatibility

In the general three-dimensional case, the state of strain at a point is defined in terms of the six-dimensional strain vector:

In plane strain, the out-of-plane strain as well as two of the shear strains are zero. Let be the out-of-plane direction. We then have:

Assuming small displacements, the strains are related to the displacements by

where are the displacements. This may also be written as

where is given by Eq. 2. The fact that serves a dual purpose: as the equilibrium operator and as the strain-displacement operator is fundamental to much of the discussion that follows and indeed to much of the theory behind OPTUM G2.
The kinematic boundary conditions may be expressed as

where are the boundary displacements.

1.3 Principle of virtual work

In what follows, extensive use will be made of the principle of virtual work. This principle may be stated as follows. Consider a stress field, , that satisfies the static equilibrium and boundary conditions:

Furthermore, consider a displacement field, , and a strain field related by:

The following identity then holds:

This is the principle of virtual work. It is emphasized that the stress and displacement/strain fields are not necessarily related. Indeed, the principle of virtual work makes no assumptions on constitutive behaviour at all.
A different statement of the principle of virtual work is as follows. Let be a stress field satisfying

for all displacement and strain fields with


This statement of the principle of virtual work is at the heart of the standard finite element method.  


The strength and deformation characteristics of geomaterials are usually accounted for by a combination of elasticity and plasticity as summarized in the following.

2.1 Stress and yield

The fundamental premise of plasticity theory is that there exists a threshold beyond which stress states cannot exist. This limit is defined by a yield function such that defines the domain of permissible stress states while defines the yield surface (see Figure 2). The yield surface may be open or closed but is always convex.

Figure 2: Yield surface.

2.2 Strain and flow

Classic plasticity theory operates with the postulate of an additive decomposition of the total strain into elastic and plastic parts:

where are the total strains, are the elastic strains and are the plastic strains.
The elastic strains are related to the stresses via a relation of the type

where is an elastic compliance modulus.
The rate of change in plastic strain is usually specified via a flow rule of the type

where is a plastic multiplier and is a flow potential. While the value of is unknown a priori, it must be such that plastic strains occur only for , i.e. when the stress state is at yield. This is guaranteed by the condition:

Regarding the flow potential, there are in principle no limits to the form of and whatever fits the experiments the best may be chosen. However, from a mathematical point of view, there are significant benefits from choosing . The flow rule is then referred to as being associated as opposed to the general case where . Also, in practice, even if and are not identical, is often of the same functional form as but involves a different set of material parameters.

2.3 Hardening

While the most basic material models make use of a steady yield surface, more advanced models incorporate one or more mechanisms of hardening. From the point of the view of the yield function, hardening can be modeled by assuming the dependence of set of additional stress-like hardening variables . The elastic domain is then given by

while again defines the yield surface.
The evolution of the hardening variables is specified by a hardening rule which may be written in the following general form:

where is an array of hardening functions. This rule implies that hardening occurs as a result of plastic straining only ().

2.4 Summary

The above constitutive relations can be expressed in the following concise incremental format:

where the notation has been used to denote the gradient of a function with respect to .  


A distinguishing feature of geomechanics problems is that pore pressures often play an important role and must be accounted for in detail in many analyses. In a fluid-saturated medium, the total pore fluid pressure, , can be divided into two parts:

where is the pore pressure due to seepage or a static water table and is an excess pore pressure generated in response to the deformation of the material. In the following, we will assume that pore pressures are negative, corresponding to compression and consistent with the sign convention adopted for stresses.
In principle, excess pore pressures are generated for any material in response to loading. However, for coarse grained materials, the dissipation of excess pressures is usually assumed to be sufficiently rapid for the excess pressures to be neglected at all times. For fine grained materials, on the other hand, excess pressures will initially be generated in response to loading. These will then gradually dissipate so that all excess pore pressures again are zero after a sufficiently long time. These two extreme states are usually referred to as drained and undrained respectively.

3.1 Effective stress principles

Besides the equilibrium, boundary and strain-displacement relations, the complete solution of a problem of solid mechanics requires a constitutive equation relating the stresses to the strains. For porous solids, the establishment of such relations are complicated by the fact that the pores may contain fluids that have a non-negligible effect on the deformations. The most common approach to establishing constitutive laws for fluid saturated porous media is Terzaghi’s effective stress principle. This principle states that the appropriate constitutive relations should be established between the strains and the effective stress, , defined by

where are the total stresses and . For example, for an linear elastic material, the relevant constitutive law is given by rather than . It should be borne in mind that while Terzaghi’s effective stress principle has proved highly successful as a basis for modeling many materials of practical interest, there is nothing ‘fundamental’ about it. Indeed, it may be viewed as part of the constitutive law and for some materials, notably rock, concrete and similar materials, the effective stress principle proposed by Biot is often more appropriate. An in-depth discussion of this issue is provided by .

3.1.1 Undrained compressibility

In general, problems of fluid-saturated porous media involve the determination of the total stresses, , as well as the excess pore pressures, . The latter may in many ways be viewed as an additional stress variable. As such, an additional constitutive relation is required (in addition to the equations relating the six stress components to the six strain components). Since only volumetric deformations can arise as a result of changes in pore pressure, the additional constitutive relation must necessarily be of the type

where is the volumetric strain and is the effective bulk modulus of the solid/fluid mixture. This material constant may be estimated in various ways, but can be shown to fall within the bounds provided by the harmonic and arithmetic between the bulk moduli of the fluid and of the solid phase material (the grains):

where and are the bulk moduli of the fluid and solid phases respectively and is the volume fraction of the fluid phase (the porosity). Considering that is usually much larger than , we have

where the lower bound usually is the more accurate one. The point, though, is that for most materials is sufficiently large for the medium to be considered incompressible, i.e. . As an example, for water at atmospheric conditions, GPa while for quartz, is of order GPa. With a typical porosity of , we have a lower bound of GPa, which is several orders of magnitude larger than the effective bulk modulus, , of clays and sands measured under drained conditions. Consequently, is often assumed.

3.1.2 Unsaturated conditions

Under partially saturated (unsaturated) conditions, the pore water pressure is in tension corresponding to . It is common practice to ignore such pressures when calculating the effective stress. In other words, the effective stress is defined as:

In OPTUM G2, this definition is denoted Terzaghi effective stress. Alternatively, it is possible to use the definition originally due to Bishop:

where is a parameter. The basic idea is to specify this parameter such that Terzaghi’s effective stress is recovered in the limits of a fully saturated and a completely dry material. The choice satisfies this requirement and is used in conjunction with the Bishop effective stress. The switch between Terzaghi (default) and Bishop effective stress can be made under Project/Physical Parameters.

3.2 Equilibrium equations

The native total stress equilibrium and boundary conditions are given by:

where the total tractions have been split into two parts: one, , due to externally applied mechanical loading and another, , due to steady-state pore pressures (for example those acting on fluid submerged boundaries).
Alternatively, Eq. 8, the equilibrium and boundary conditions may be expressed in terms of effective stresses and pore pressures as

3.3 Strain-displacement relations

The strain-displacement relations are unaffected by the presence of a pore fluid and are given by

as usual.

3.4 Principle of virtual work

The principle of virtual work is given by

or, alternatively, by

where in both cases satisfies the kinematic boundary conditions.  

3.5 Summary of governing equations

The full set of governing equations for drained and undrained conditions are summarized in the tables below. We have here used and assumed perfect compressibility, .

Table 1: Governing equations for quasi-static mechanical problem under drained and undrained conditions.
 Drained  Undrained
Eff. stress princip.
Static BC
Kinematic BC
Incompressibility  –
Table 2: Governing equations for quasi-static mechanical problem under drained and undrained conditions (effective stresses and excess pressures as primary variables).
 Drained  Undrained
Static BC
Kinematic BC
Incompressibility  –



Most problems in geomechanics involve pore pressures in the form of either seepage pressures or excess pore pressure, or both. While the governing equations summarized in the previous section are well established, the actual implementation of these may be undertaken in a number of different ways. Sometimes a particular choice is motivated solely by convenience while, at other times, it is motivated more by the actual physics. In the following, some of the more common modeling approaches are discussed along with some of the consequences associated with the generation of excess pore pressures.

4.1 Drainage Conditions and Time Scope

In principle, excess pore pressures are generated for any material in response to loading. With time, these excess pressures will dissipate, i.e. the imbalance created in the total pore pressure distribution resulting from excess pressure generation will tend to a new equilibrium state following Darcy’s law. In practice, this equilibrium state is given by the state that existed before the application of the load, i.e. a state with zero excess pore pressures.
For coarse grained materials, the dissipation of excess pressures is usually assumed to be sufficiently rapid for the excess pressures to be neglected at all times. For fine grained materials, on the other hand, excess pressures will initially be generated in response to loading. These will then gradually dissipate so that all excess pore pressures again are zero after a sufficiently long time. The two different situations are illustrated in Figure 3. The coarse grained material never experiences a significant build-up of excess pore pressure. As such, the conditions are said to be drained at all times. The fine grained material, on the other hand, experiences significant excess pressures immediately after the application of the load. These excess pressure remain at close to their initial value for an appreciable amount of time under which the conditions are said to be undrained. After a sufficiently long time, the excess pore pressure will have dissipated and the conditions are again drained.

Figure 3: Response of fluid saturated coarse and fine grained materials to external loading.

In OPTUM G2 , the term ‘sufficiently long time’ is quantified by two settings:

  • Time Scope: Short Term or Long Term. This is a stage setting and indicates whether the loading and/or changes to the geometry happen instantaneously or over an extended (in principle infinitely long) period of time. Generally speaking, the structural integrity of geotechnical structures have to be verified with respect to both scenarios. For example, a foundation must be stable both when the load is first placed (possibly leading to generation of excess pressures) and after a long time (when the excess pressure have dissipated).
  • Drainage Conditions: Always Drained, Drained/Undrained, Non-Porous. This is a material setting that indicates whether excess pressures are generated or not. For Always Drained materials, excess pressures are never generated, regardless of the Time Scope. Gravel, sands and similar coarse grained materials are examples of such materials. For Drained/Undrained materials, excess pore pressures will be generated for Time Scope = Short Term, but not for Time Scope = Long Term. Clays and similar fine grained materials follow this pattern of behavior. Finally, for Non-Porous materials, neither seepage nor excess pressures exist (this setting may be relevant to some types of rock, concrete, etc).

The rules for whether a given point in the domain behaves in a drained or an undrained manner are summarized in the table below.

Figure 5: Material behavior as function of Drainage and Time Scope.

4.2 Choice of primary variables

Whenever pore pressures are involved, a sharp distinction must be made between total and effective stresses. Compared to the situation where no pore pressures are present or generated during the analysis, the equilibrium equations are unaltered and still involve the total stresses and the total body forces. The static boundary conditions are similarly unaffected although it should be noted that boundary tractions due to water pressures comprise part of the total tractions (see Section 1). The real difference is that the constitutive relations involve the effective rather than the total stresses. This situation leads to at least three different possibilities for implementation:

  • Equilibrium and boundary conditions cast in terms of total stresses, constitutive relations in terms of effective stresses, and the effective stress principle embedded as a set of constraints. This is the option shown in Table 1.
  • All governing equations cast in terms of effective stresses and excess pore pressures. These are shown in Table 2.
  • All governing equations cast in terms of total stresses and excess pore pressures.

OPTUM G2 makes use of the second of these options.

4.3 Undrained total stress analysis

In some particular cases, it is possible to conduct the analysis entirely in terms of total stresses. In particular, for undrained analysis using the Mohr-Coulomb or Drucker-Prager models, it is under certain circumstances possible to conduct the analysis in terms of total stresses using undrained material parameters that differ from but are related to the original, drained, parameters.
As an example, consider the linear elastic/perfectly plastic Drucker-Prager model. This model involves Hooke’s law with parameters and and the yield and plastic potential functions:

where , , and are the friction coefficient, dilation coefficient, and cohesion respectively and and are the mean and deviatoric stresses respectively (see the Materials Manual for details).

4.3.1 Elasticity

Assuming a finite fluid compressibility, the elastic stress-strain relations are given by

where it has further been assumed that the seepage pressures do not change, . Straightforward manipulations together with the assumption that then lead to the following relation between strains and total stresses:


In other words, in terms of strain versus total stress, the elastic behaviour is identical to that of a drained material with a Young’s modulus and Poisson’s ratio . Alternatively, the material may be thought of as having an unaltered shear modulus and a Poisson’s ratio of . In other words, the modeling may be carried out in terms of effective stresses and excess pore pressures using Eq. 12 or in terms of total stresses using Eq. 13. In both cases, we have

That is, the effective mean stress remains at its initial value throughout.

4.3.2 Plasticity

Due to the incompressibility of the medium, the change in volumetric strain is zero. Below yield, this means that the change in elastic volumetric strain is zero while, at yield, the sum of the elastic and plastic volumetric strains are zero:

where is the plastic volumetric strain rate:

The change in excess pore pressure is thus given by

Bearing in mind that is very large, any small plastic volumetric strain will lead the excess pressure to decrease significantly. For a fixed total mean stress, this implies a decrease in effective mean stress and thereby effectively an increase in shear strength. In typical elastoplastic calculations, this means that a limit load will never be attained. The exception is the case where the plastic strain rate, by virtue of the plastic potential, is zero. In that case, we have at yield, and the effective mean stress remains at its initial value.
For the Drucker-Prager criterion, we can therefore conclude that a dilation parameter of should be used for the results to be meaningful. Moreover, since the effective mean stress never changes from its initial value, the yield function can be written as:

where is an effective ‘undrained strength’ parameter.
In summary, undrained elastoplastic analyses can be carried out in terms of total stresses using the following quantities:

4.3.3 Mohr-Coulomb

The same analysis as above can be carried out with respect to the Mohr-Coulomb criterion, resulting in the following undrained quantities:

The yield and plastic potentials are here the Tresca function with strength parameter:

where and are the Mohr-Coulomb parameters measured under drained conditions, is the initial earth pressure coefficient, and is the initial effective vertical stress (positive in compression).
The parameter is referred to as the undrained shear strength. As it appears above, it is a function of the drained parameters and the initial stress state. However, in practice, one may use the a distribution of that is more in line with what is measured in the field. In OPTUM G2, the Tresca model may be used for this purpose.


Variational principles are at the heart of OPTUM G2 and all problems are cast – and solved – as such. The following sections detail the fundamental concepts of variational principle and the specifics of their use in OPTUM G2.

5.1 What are variational principles?

A variational principle can be thought of as an optimization problem that offers an alternative way of stating the governing equations of a physical system. As a simple example, consider Hooke’s law:

An alternative statement of this governing equation is given by

The proof of the equivalence between Eq. 16 and Eq. 17 follows readily by differentiating with respect to and setting the result equal to zero.
A more complex example is the linear elastic truss shown in Figure 6.

Figure 6: Linear elastic truss.

Using standard finite element terminology, the governing equations can be written as

where is the stiffness matrix and and are arrays of nodal displacements and nodal forces respectively. For the time being, let all displacements be unknown and the nodal force array specify only the externally applied loads. That is, the system is not yet supported. The above set of equations may be viewed as a multidimensional generalization of Hooke’s law and an equivalent optimization – or variational – formulation is given by

Boundary conditions can now be accounted for simply by imposing appropriate constraints. Consider first the two supports on the left. The displacements are here zero or, more generally prescribed to some values which in this case is zero. The above optimization problem should thus be extended to

where the matrix contains zeros and ones and the vector contains the prescribed displacements (for example zero).
This problem may be solved using the Lagrange multiplier technique. First define the Lagrangian which is the objective function (the quantity to be optimized) augmented by penalty terms that ensure satisfaction of the constraint:

where is a set of Lagrange multipliers. It may be shown that minimizing the Lagrangian solves the original optimization problem. In other words, the gradients of with respect to the variables at hand, and , are required to vanish:

These equations are known as the optimality conditions associated with the optimization problem Eq. 21: solving them solves the optimization problem. Physically, we see that the first set of equations are the original governing equations now augmented by the term where it is clear that the Lagrange multipliers are the reactions. The second set of equations are simply the kinematic boundary conditions imposed in the original optimization problem.
A more complex situation is the one suggested by the support at a distance from the lower part of the truss. This support calls for an inequality constraint of the type where is the vertical displacement of the node above the support. More generally, we may extend the optimization problem with general inequality constraints to arrive at:

Or equivalently, introducing so-called slack variables, :

Again, the Lagrange multiplier method may be used to solve the problem. The Lagrangian is given by

where is an additional set of Lagrange multipliers. The optimality conditions are given by

where are recognized as the reactions associated with the inequality constraints. In addition, and crucially, the non-negativity requirement on imposes the following complementarity conditions:

where is the length of the arrays. These conditions have the following physical significance. Suppose that . The corresponding inequality constraint is then satisfied with equality. In other words, the distance between the truss and the support is zero. A reaction, i.e. a value of is then possible. Conversely, if , there is still a gap between the truss and the support and .
Comparing the three types of optimization problems considered above – the unconstrained problem Eq. 19, the equality constrained problem Eq. 20 and the inequality constrained problem Eq. 22 – it is the latter problem that, due to the nonlinearity of the complementarity conditions, is the most difficult to solve. Indeed, while the two first can be solved directly, the solution of optimality conditions associated with Eq. 22 requires an iterative procedure. In OPTUM G2, this task is carried out using the general purpose optimization solver SONIC.  


In analogy with the developments outlined in the previous section, it is possible to obtain optimization formulations for a broad range of practically relevant problems. The first one considered is that of structures of rigid-plastic material. For the sake of notational convenience, the effects of pore pressures are initially neglected. The main results, with the effects of pore pressures included, are summarized in Section 2.

6.1 Governing equations

Rigid-plastic materials undergo no deformations below the point of yield whereas, at the point of yield, unlimited plastic deformation takes place. As such, the governing equations should be formulated in terms of displacement rates (velocities) and strain rates, rather than total displacements and strains. Furthermore, we will assume that the deformation up to the point of collapse are sufficiently small to ignore changes in geometry.
In terms of statics, the governing equations first of all comprise static equilibrium and boundary conditions:

Secondly, the yield condition must be satisfied:

From a kinematic point of view, the associated flow rule is assumed appropriate:

where are the plastic strain rates and are plastic multipliers. Using the small deformation assumption, these should be derived from a velocity field as

Or, combing the above two equations:

Finally, the complementarity condition should be satisfied:

meaning that yielding () can only take place if the yield condition is satisfied [].

6.1.1 Linearization

For subsequent developments, it is convenient to consider a linearized yield function rather than the original nonlinear one. In other words, is replaced by a set of linear constraints of the type:

or, in matrix form:

where and collects the contributions, and respectively, from each linear constraint. Alternatively, introducing slack variables, the yield constraints may be written as

It is self-evident that the original nonlinear yield function can be approximated to within an arbitrary degree of accuracy by increasing the number of linear planes (although it may not be practically feasible to do so).
The flow rule associated with the linearized constraints is given in terms of ‘Koiter’s rule’ by

where contains the plastic multipliers associated with each of the linear constraints. This rule follows as an obvious consequence of von Mises’s principle of maximum plastic dissipation and it is difficult to imagine any other way of defining the plastic strain rates for a composite yield surface.

6.2 Limit analysis

Consider now the following situation. A structure of rigid plastic material is subjected to a set of body forces , stemming for example from self weight. On its boundary, a set of tractions, , act. For such a scenario, the central question of limit analysis can be posed as: What is the maximum magnitude of the tractions that can be sustained without the structure suffering collapse. Or alternatively: what is the minimum magnitude of the tractions that will cause collapse.

Figure 7: Solid of volume with boundary subjected to tractions on and supported on .

6.2.1 Complete solution

Let us introduce a load multiplier such that the tractions acting on the structure are given by (see Figure 7). Suppose further that the structure is at collapse. The displacements are then infinite, so it is necessary to introduce a scaling of the velocities or a relevant work quantity, or similar. This is equivalent to what is done in many hand calculation upper bound methods where a finite magnitude of some characteristic displacement defining the collapse mechanism is considered. The governing equations are then given by:
Equilibrium and static boundary conditions:

Yield conditions:

Associated flow rule/strain-displacement compatibility:


Complementarity conditions:

where the scaling has been applied with respect to the rate of work done by the reference tractions .
It may be shown that the solution to the above equations, if it exists, is unique in terms of the multiplier . However, there may be more than one stress distribution or velocity field leading to the same value of the collapse multiplier. The above governing equations may be stated alternatively in terms of a number of variational principles that in some cases allow for the establishment of bounds to the exact collapse multiplier to be determined.

6.2.2 Lower bound principle

One possibility of stating the governing equations is in terms of the lower bound principle:

In other words, the solution to the above problem satisfies the governing equations Eq. 25–Eq. 29. The kinematic quantities, which are absent from the above problem, appear as Lagrange multipliers when solving the problem. This is similar to the situation with the reactions of the truss in Section 5.
The main strength of the lower bound principle is that it allows for a lower bound on the exact collapse multiplier to be computed, namely by constructing a stress field that satisfies the constraints without necessarily being optimal.

6.2.3 Upper bound principle

Secondly, the governing equations may be cast in terms of the following optimization problem:

This problem requires consideration of kinematic quantities and provides the possibility to compute an upper bound to the exact collapse multiplier, namely by postulating a compatible velocity field that satisfies the flow rule. This is done in such a manner that the rate of work done by the reference tractions is scaled to unity. The objective function, which comprises the internal rate of work minus the contribution from the constant body forces, is then the collapse multiplier sought.

6.2.4 Bounds

To verify that the lower and upper bound principles indeed do furnish bounds on the collapse multiplier, we may proceed in the following manner. Consider first a stress field, , satisfying the yield condition and the equilibrium and boundary conditions with load multiplier . The principle of virtual work gives

where is taken as the exact velocity field and .
Moreover, consider the exact stress field, , in conjunction with the exact velocity field:

where is the exact collapse multiplier.
Subtracting Eq. 32 from Eq. 33 and using Eq. 28 gives

where the last inequality, which hinges crucially on the associated flow rule, is illustrated in Figure 8.

Figure 8: Illustration of lower bound inequality Eq. 34.

Proceeding to the upper bound principle, consider a velocity field and a plastic multiplier field , not necessarily related to the velocity field. Furthermore, consider a stress field, , not necessarily in equilibrium, but satisfying the yield conditions . We may then define a collapse multiplier by

The exact collapse multiplier is defined via the exact stress field:

where . We thus have:

Invoking the associated flow rule, for the assumed displacement field, leads to:

where the inequality again follows from a simple geometric construction similar to the one shown in Figure 8.

6.2.5 Duality

The two problems Eq. 30 and Eq. 31 are in fact dual to each other in the sense that i) they are constructed on the basis of the same data (, , , and ) but involve different sets of variables, and ii) their solutions in terms of the objective function are identical. This latter point is verified as follows. Consider a stress field, , satisfying the yield condition , and the equilibrium and boundary conditions with load multiplier . Consider further a velocity field, , related to the plastic strain rate and plastic multiplier fields by , and satisfying . We then have:

where the duality gap vanishes if the complementarity conditions are satisfied.
Besides demonstrating the duality between the two principles, the above manipulations also formally prove that each of the two principles are equivalent to the full set of governing equations.

6.2.6 Role of associated flow rule

While the assumption of associated flow is essential to the framework of limit analysis, it has rather different implications for each of the bounding principles. Thus, in order to prove the lower bound inequality Eq. 34, it is necessary to assume that the exact solution satisfies the flow rule while no assumptions are made regarding the approximate solution. Indeed, this latter solution operates only with static variables whereas no assumptions on the kinematics are made. If the exact solution does not satisfy the associated flow rule, an equilibrium stress solution is still valid in the sense that it satisfies all the basic requirements of continuum mechanics. However, the corresponding load multiplier cannot be guaranteed to be below the collapse multiplier associated with the exact stress distribution.
In contrast, the upper bound inequality requires that the approximate solution satisfies the flow rule while no requirements are attached to the exact solution. This means that the inequality is valid even if the exact solution does not assume an associated flow rule. In other words, an upper bound to a problem with an associated flow rule will also be an upper bound to the equivalent problem with a nonassociated flow rule.
In section 4, the role of the flow rule, and its implication for the collapse load, are discussed in more detail.

6.3 Pore pressures

The effects of seepage and excess pore pressures are readily included into the lower and upper bound principles, Eq. 30 and Eq. 31 respectively.

6.3.1 Drained conditions

Assuming drained conditions, the lower bound principle reads:

while the upper bound principle is given by:

In the above it is assumed that only the external load, , is amplified to attain a state of collapse while the part of the tractions due to seepage pressures, , remains constant.

6.3.2 Undrained conditions

Assuming undrained conditions, the lower bound principle reads:

while the upper bound principle is given by:

Again, it is assumed that only the external load, , is amplified to attain a state of collapse while the part of the tractions due to seepage pressures, , remains constant.

6.4 Gravity multiplier limit analysis

In some cases, it is useful to attain a state of collapse by amplifying gravity while keeping all external tractions constant. The resulting collapse multiplier may in some applications, e.g. slope stability, be interpreted as the factor of safety. The above principles are easily modified to account for this scenario, with the relevant lower bound principle being:

We here note that both the body forces, , as well as the seepage pressures, , and the tractions due to seepage pressures, , are amplified proportionally consistent with the scenario that gravity is amplified to attain a state of collapse.

6.5 Mixed principles

The key feature of the upper and lower bound theorems is that they facilitate the computation of rigorous bounds. A number of other principles, which again reproduce the governing equations, are possible. Rather than facilitate the computation of rigorous bounds, the key advantage of these principles is that they suggest various ‘compromise’ solutions, i.e. solutions that cannot be rigorously bounded but which tend to be closer to the exact solution than either of the rigorous upper and lower bound solutions. These principles are generally referred to as mixed principles and typically involve both the stress and the displacements as primary variables.
A commonly used mixed principle is the following:

where both the velocities and the stresses appear as variables. This principle is particularly suitable as a starting point for mixed stress-displacement finite element formulations and has been used extensively for that purpose .
The gravity multiplier version of the principle is given by:

where both body forces and seepage pressures again are amplified.  


Variational principles for boundary value problems of linear and nonlinear elasticity have long been available and have often been used as a basis for the construction of finite element approximations. In the following, a number of such principles are presented.

7.1 Governing equations

The governing equations comprise:
Equilibrium and static boundary conditions:

Strain-displacement compatibility:

Constitutive law:

where and are the Helmholtz free and complementary energy functions respectively. These functions are assumed convex and are related by

For linear elasticity we thus have and , where and are the elastic compliance and stiffness moduli respectively. Moreover, with we have

7.2 Complementary energy principle

The governing equations Eq. 44–Eq. 46 can be shown to be reproduced by the following optimization problem:

This principle, which is known as the total complementary energy principle, can be viewed as the elastic analogue to the lower bound principle Eq. 30. Even though it operates only with stresses, the displacements are recovered in the process of solving the problem, namely as the Lagrange multipliers associated with the equilibrium constraints.

7.3 Potential energy principle

Alternatively, the governing equations may be cast in terms of the total potential energy principle:

This form bears some resemblance to the upper bound principle Eq. 31. It is the form most commonly used as a basis for the finite element method.

7.4 Hellinger-Reissner principle

Finally, the Hellinger-Reissner principle is often useful for constructing non-standard finite element formulations. This principle is given by

In other words, the governing equations are reproduced by minimizing the objective function with respect to the displacements and maximizing it with respect to the stresses.

7.5 Bounds

The principles Eq. 49 and Eq. 50 provide bounds on the energy of the system. This property may be used to either calculate upper and lower bounds on properties of interest or as a means of gauging the proximity of approximate solutions to the exact solution.
Before proceeding, recall that a convex function is one where the function lies above all of its tangents:

for all and (see Figure 9).

Figure 9: Convex function.

Considering Eq. 49, any admissible, though not necessarily optimal, stress distribution satisfies the principle of virtual work:

where is the exact displacement field and . Similarly, the exact stress distribution, , satisfies:

Subtracting Eq. 54 from Eq. 53 and using Eq. 46 and Eq. 52 gives


thus demonstrating the lower bound nature of Eq. 49. Generally speaking, we may say that since the internal energy is over-estimated, so is the external energy and for fixed loads, the displacements implied by the principle of virtual work are larger than the exact displacements.
Proceeding to Eq. 50, we consider the principle of virtual work for an admissible displacement field in conjunction with the exact stresses:

and the exact displacements in conjunction with the exact stresses:

Subtracting Eq. 56 from Eq. 55 and using Eq. 46 gives

Using Eq. 52, we have

which verifies the upper bound nature of Eq. 50. This inequality implies that approximate solutions generally speaking will be ‘too stiff’ in the sense that the displacements, for a given set of external forces, are underestimated. This is a well known property of the standard displacement finite element method based on the principle of minimum potential energy.

7.6 Duality

Finally, we show that the problems Eq. 49 and Eq. 50 are dual to each other in the sense they are constructed on the basis of the same data and have the same solution. We consider two admissible, but not necessarily related, fields and . Subtracting the objective function of Eq. 49 from that of Eq. 50 gives

where the final inequality follows from the results of the previous section. We see that the gap vanishes when both the displacements and the stresses are the exact fields.  


In the following, a general thermomechanical formulation of elastoplasticity is first considered after which a number of useful variational principles are derived and discussed. For the sake of notational convenience, the effects of pore pressures are initially neglected. The main results, with the effects of pore pressures included, are summarized in Section 3.

8.1 Thermomechanical formulation of elastoplasticity

In this section a thermomechanical formulation of elastoplasticity is briefly summarized following primarily the exposition of and .
From the first and second laws of thermodynamics and assuming isothermal conditions, the following central identity can be derived :

where is the rate of internal work, is the Helmholtz free energy function and is the dissipation function. Following we will assume that the Helmholtz free energy is a function of the elastic strains, , and a set of strain-like hardening variables, :

The time derivative of is given by

By comparison with Eq. 58 for a purely elastic process (, ), it can be verified that the stresses are given by

The set of variables conjugate to are referred to as stress-like hardening variables and are denoted by :

Next, the complementary Helmholtz free energy function, , is defined via the Legendre transformation:

from which it follows that the elastic strains and the strain-like hardening variables can be expressed as

From straightforward manipulations of the above relations the dissipation follows as

which is a well-known result [except that a change of sign of is often made ]. Furthermore, the rate of internal work may be expressed in the following form:

Following , we seek to derive constitutive models by maximizing the internal work rate subject to yield conditions:

From the expression Eq. 61 of the internal work rate, we can conclude that

The relevant maximization principle thus reduces to

The procedure is now to postulate a relevant potential after which the constitutive equations follow as the first-order optimality conditions associated with the above maximization problem. In this connection, we note that von Mises’ principle of maximum plastic dissipation appears for the particular choice of .
In order to derive more general constitutive equations, the time derivative of is first expanded as

Using this expansion, the constitutive equations associated with Eq. 63 follow as:

It is convenient to introduce the effective moduli:

so that the constitutive equations can then be expressed in the following format:

The moduli and here have the same physical significance as in an conventional elastoplastic model of the kind summarized in Section 1: is an elastic compliance modulus and is an array of hardening functions. In addition, the governing equations include a new constitutive modulus, , referred to as a coupling modulus. This additional modulus is useful in the construction of constitutive models for soils and other granular materials .

8.1.1 Finite-step formulation

A finite-step version of the incremental variational formulation discussed above is considered.
The power of deformation may be approximated in time as

where subscript refers to an initial, known, state. Using the relations of the previous section together with a fully implicit evaluation of the dissipation, the finite-step power of deformation may be expressed in terms of and as


The finite-step version of the material point principle Eq. 62 thus reads:

where is assumed known.

8.2 Variational principles

Following the treatment of rigid-plasticity, a number of similar variational principles may be derived.

8.2.1 Hellinger-Reissner principle

The most straightforward variational principle appears by extending the local material-point principle Eq. 67 to the entire domain. The resulting principle reads:

For linear elasticity, this principle reduces to the well-known Hellinger-Reissner principle.

8.2.2 Lower bound principle

An equivalent statement of the governing equations is the lower bound principle:

where we note that the rigid-plastic lower bound principle Eq. 30 appears for .
As a special case, consider a linear elastic-perfectly plastic material with linear (or linearized yield functions). Further assume that the elastic modulus is time-independent (). We then have:

where , and thereby:

where again are slack variables. The similarity of this principle to the rigid-plastic principle Eq. 30 is apparent. Indeed for (rigid elastic behaviour) the latter principle is recovered exactly.

8.2.3 Upper bound principle

Finally, a dual, upper bound version of Eq. 69, may be derived. In the general case, the result is not particularly transparent and involves quantities that do not follow immediately from the main quantities introduced so far (the potential and the yield function ). Instead, consider the special case Eq. 70. The dual to this principle is given by:

where are to be interpreted as the elastic strains, are the displacements, and are plastic multipliers. The similarity to the rigid-plastic upper bound principle Eq. 31 is again apparent. Indeed, for (rigid elastic behaviour), the minimization nature of the problem will imply .

8.3 Pore pressures

The effects of seepage and excess pore pressures are readily included into the principles derived above.

8.3.1 Drained conditions

Assuming drained conditions, the lower bound principle Eq. 69 extends to:

while the modified Hellinger-Reissner principle reads:

In the above it is assumed that only the external load, , is amplified to attain a state of collapse while the part of the tractions due to seepage pressures, , remains constant.

8.3.2 Undrained conditions

Assuming drained conditions, the lower bound principle is given by:

while the modified Hellinger-Reissner principle reads:

Again, it is assumed that only the external load, , is amplified to attain a state of collapse while the part of the tractions due to seepage pressures, , remains constant. We note that the excess pressure, , enters the lower bound principle as an additional variable while the upper bound principle includes an additional constraint ensuring incompressibility.  


The variational principles derived in the preceding sections all rely crucially on the concept of associated flow. For some materials, notably metals, this assumption is in reasonable agreement with experiments and the principles are thus directly applicable. For frictional materials, on the other hand, the assumption of associated flow is usually in contradiction to experiments. Indeed, the dilation implied by the flow rule associated with the Mohr-Coulomb, Drucker-Prager or other relevant failure criteria is usually far greater than actually observed in experiments. A nonassociated flow rule is therefore often used, usually by specifying a flow potential of the same functional form as the yield potential .

Figure 10: Soil behaviour and capabilities of associated and nonassociated Mohr-Coulomb models. Note: the nonassociated response assumes that no localization takes place.

Using a nonassociated flow rule, it is possible to model the real behaviour of soils fairly well, even with a linear elastic/perfectly plastic model. The situation is outlined schematically in Figure 10. For loose sands and normally consolidated clays, the response in typical drained triaxial or biaxial test is such that the shear stress increases monotonically up to some ultimate limit state. During this process, the soil will undergo a compaction, the rate of which will eventually tend to zero. With the Mohr-Coulomb model implemented in OPTUM G2, it is possible to account both for the initial approximately linear stress-strain response and for the ultimate strength. In the intermediate range, the response of the Mohr-Coulomb model tends to be somewhat too stiff due to the assumption of perfect, i.e. non-hardening, plasticity. These characteristics are independent of whether an associated or a nonassociated flow rule are used. However, the deformations – here shown in terms of volumetric strain versus shear strain – are highly dependent on the flow rule. Indeed, for the present example, the associated flow rule predicts a dilation that is far too large. The nonassociated flow rule, on the other hand, specifies the amount of dilation directly and it is thus possible to account for the deformations quite accurately.
For dense sands and overconsolidated clays the situation is, qualitatively, much the same. The stress-strain response is independent of the flow rule while the associated flow rule over-predicts the amount of dilation significantly. We also note the appearance of a peak followed by a softening in the stress-strain response. This coincides with a tendency for the dilation to decrease, eventually reaching a value of zero.
On the basis of these observations, it would be tempting to conclude that the yield function and the flow rule can be adjusted independently – the former to match an observed strength the latter to match an observed deformation behaviour. While there is some truth to this, it is unfortunately somewhat naive. In fact, the introduction of a nonassociated flow rule gives rise to a number of phenomena that do not reveal themselves from the kind of simple observations made above. This is discussed in the following section.

9.1 Consequences of nonassociated flow

The question of how the response of plastic or elastoplastic boundary value problems are affected by the flow rule is an old one that still in many ways has not been resolved. Since the effect of the flow rule on the deformations is quite obvious, at least qualitatively, the focus has mostly been on how the flow rule affects the limit load. On the basis of the kind of behaviour illustrated in Figure 10, it could be concluded that there is no effect at all: the ultimate limit load is independent of the flow rule. This view is still upheld in some circles. A slight moderation of it is that the flow rule only affects problems with a ‘high degree of kinematic constraint’. However, exactly how the ‘degree of kinematic constraint’ is quantified has never been elucidated other than occasional statements that it for most soil mechanics problems is ‘low’. Such statements are only rarely backed up by quantitative results, for example in the form of finite element analyses for problems with different friction and dilation angles. And if they are, the problem and/or the set of material parameters chosen are usually of little practical significance. On the other hand, where more rigorous and practically relevant numerical analyses have been undertaken, the conclusion has inevitably been that the flow rule has a rather significant effect on the ultimate limit load. For example , and report reductions in bearing capacity of up 45% for strip and circular footings on sand. Such reductions in bearing capacity find abundant theoretical support as do a number of other perhaps unexpected phenomena that are observed in typical numerical analyses. Some of these are discussed below.

9.1.1 Consequence 1: Localization and non-uniqueness of limit load

An implicit assumption usually made when solving elastoplastic boundary value problems is that the response, or at least the limit load, is unique. This assumption does indeed hold when the flow rule is associated as was shown in Sections 2 and 3 for rigid-plastic and elastoplastic materials respectively. However, when the flow rule is nonassociated, the response will generally be non-unique. In other words, there may be multiple solutions all satisfying the governing equations, but each implying a different limit load. This fact and the tools with which to analyze constitutive models with respect to uniqueness have long been established , but have found relatively little resonance in the geomechanics community.

Figure 11: Response of biaxial test prone to localization.

The type of non-uniqueness experienced for nonassociated materials is intimately tied to the appearance of shear bands, i.e. localized bands of intense deformation that tend to form rather abruptly and usually lead to an apparent softening in the load-displacement curve. The phenomenon can be illustrated as follows. Consider a biaxial test as shown in Figure 11. The test is modeled by means of a constitutive model that makes use of a nonassociated flow rule. Initially, the response will be homogeneous, i.e. the stress and strain states are constant throughout the sample. Then at some point, a bifurcation is possible. This means that the response may continue along a path where the deformation remains homogeneous or the response may change fundamentally by the sample developing a shear band. This change of deformation mode will lead to a drop in the axial force that can be sustained as the loading is continued. As such the response is non-unique: there are two solutions, a homogeneous and a localized, that both satisfy all of the governing equations.
Furthermore, suppose that the sample had remained homogeneous beyond the first bifurcation point. This does not prevent bifurcation from taking place at some other point during the loading. And the geometry of the localized solution (here given by the shear band inclination) need not be the same as for the first bifurcation point. In fact, if the sample remains homogenous beyond the first bifurcation point, there will usually be a range of possible localized solutions.
It may be shown that a suitable measure of the tendency for a constitutive model to lead to localization can be gauged by the eigenvalues of the acoustic tensor defined by

where is the elastoplastic constitutive modulus and is the usual projection matrix Eq. 4 with the normal being the normal to the shear band (see Figure 11). Thus, for a given stress state, , and an assumed shear band orientation, , localization is possible if the determinant of is less than or equal to zero. If not, the sample is stable and the response will continue along the homogenous path. This and similar results are discussed in detail by among others.

9.1.2 Consequence 2: Reduction in limit load

The second consequence of nonassociated flow is that the limit load usually observed in numerical experiments is significantly lower than the corresponding limit load obtained under the assumption of associated flow. Generally speaking and with reference to the Mohr-Coulomb model, the reduction in limit load depends on the difference between the friction and dilation angles and on the magnitude of these angles. Thus, the reduction in limit load increases as the difference between the friction and dilation angles increases. And for a fixed difference, the reduction usually increases with increasing absolute magnitude of the friction and dilation angles. In other words, the parameter set usually leads to a greater reduction in bearing capacity that the set .
The causes of this reduction are two-fold. Firstly, the tendency of nonassociated flow rules to induce localization will in far most cases mean that the ultimate limit state will be governed by highly localized states of stress and strain. Secondly, the flow rule imposes certain constraints on the kinematics of the shear bands that in the nonassociated cases leads to a reduction in bearing capacity.

Figure 12: Rigid block on a rigid frictional surface.

The second effect can be illustrated by the simple example shown in Figure 12. We here consider a rigid block on a rigid frictional surface. The block is subjected to a constant normal force and the objective is to determine the maximum tangential force, , that can be sustained. Regarding the exact geometry of the problem, we imagine an arbitrarily thin interface between the block and the surface in which the relevant yield condition is that of Coulomb:

where is the friction angle.
The flow potential, on the other hand, is given by

That is, a non-dilative behaviour is assumed.
The plastic strain rates are thus given by

Furthermore, the strain-displacement relations give

where is the velocity jump across the interface, i.e. the difference in between the upper and lower parts of the discontinuity. Similarly, is the velocity jump along the interface, i.e. the difference in between the two ends of the interface. Letting and requiring a finite magnitude of the largest strain rate, we find that . In other words, the infinitely thin interface does not stretch longitudinally.
Returning to the flow rule, this imposes the constraint that

Substituting this constraint into the yield condition gives:

from which we arrive at the solution:

rather than the standard solution of which indeed is the solution assuming associated flow, .
Introducing an ‘effective’ friction angle by , we may write the solution as

For we thus have which is similar to the ratio between the characteristic and design friction angles specified by many codes of practice, for example Eurocode 7. As such, the effect of nonassociated flow can hardly be ignored, even for the least ‘kinematically constrained’ of problems.

Figure 13: Effect of nonassociativity according to Eq. 76: in terms of friction angles (left) and in terms of friction coefficients (right).

The procedure described above can be generalized to arbitrary three-dimensional solids undergoing localized deformation . The result is that the effective strength domain is defined by a function identical to the original yield function but involves reduced material parameters, the magnitude of which depends on the friction and dilation angles. For the Mohr-Coulomb criterion we thus have an effective strength domain given by

where the effective strength parameters are given by

These relations were first derived by for the special case of and by for the general case. The influence of the degree of nonassociativity on the effective Mohr-Coulomb friction angle assuming different dependencies between the friction and dilation angles is shown in Figure 13.
Similarly, for the Drucker-Prager criterion the effective strength domain is given by :


9.1.3 Consequence 3: ‘Numerical problems’

In finite element analyses of for example a biaxial test, one will often observe that some kind of imperfection is necessary to induce localization. This may come in the form of a slight material inhomogeneity or a slight imperfection in the geometry or boundary conditions. Furthermore, it is often observed that the response tends to be rather dependent on the finite element mesh. In particular, if regular meshes are used, there is a tendency for the localization to take place along the edges of the mesh. In view of the non-uniqueness of the orientation of the localization bands – and thereby the overall force-displacement response – this is to be expected.
Secondly, and more seriously, it has frequently been reported that numerical solutions to boundary value problems involving nonassociated constitutive models are much more difficult to obtain than in the case where the flow rule is associated . These complications have a tendency to be more pronounced for high (but realistic) values of the friction angle and the degree of nonassociativity. Similarly, for fixed material parameters, one usually observes a degradation of the performance as the number of finite elements in the model is increased.

Figure 14: Load-displacement response for strip footing on a weightless soil with and [after ].

While often reported as ‘numerical problems’, these characteristics are not entirely unexpected. Since there may be a range of possible solutions, each associated with a different pattern of localization and all of which are entirely valid, it can be expected that any numerical solution will be very sensitive to both physical imperfections as well as round-off errors and the exact sequence in which the procedures defining the solution scheme are carried out. In the end, the result is a load-displacement response that tends to be rather oscillatory. An example is shown in Figure 14 which is concerned with the analysis of a strip footing on a cohesive-frictional and weightless soil (the so-called problem). Compared to the response of the corresponding associated problem we note several characteristics: i) the stiffness of the associated problem is somewhat higher, except in the very initial stages where the response is purely elastic, ii) nonassociativity leads to a significant decrease in ultimate strength (in this case by almost 40%), and iii) the load-displacement response beyond a certain level of displacement is somewhat oscillatory for the nonassociated problem. These oscillations are a consequence of the non-uniqueness of the boundary value problem and correspond physically to a switching between different modes of failure beyond the point at which the load carrying capacity of the structure first becomes exhausted. This is illustrated in Figure 15 where it is seen that the failure modes change quite significantly between different load steps beyond the point at which the limit load apparently is reached (in this case, at a displacement of ).

Figure 15: Rates of plastic shear strain at different stages of loading .

Under such conditions, standard Newton-Raphson based procedures are likely to fail. Indeed, the critical assumption behind such procedures is that a good estimate of the solution is available. Since the only real possibility for estimating the unknown solution is by the last converged solution and since the mode of stress and deformation may change very significantly and abruptly from one step to the next, convergence is not likely. This feature and the oscillations in the load-displacement curve are often grouped together and labeled ‘numerical problems’. However, whereas the former definitely is due to a shortcoming of the method of solution, the latter is a consequence of the constitutive model and is to be expected, despite perhaps being somewhat counterintuitive.

9.2 Variational formulation

Motivated by the shortcomings of conventional Newton-Raphson based solution procedures developed a new optimization based solution scheme specifically aimed at nonassociated elastoplasticity. The basic principles of this scheme are described in the following.

9.2.1 Friction and plasticity

It is well established that constitutive models of the type described in Section 1 do not permit a variational formulation unless , that is, unless the flow rule is associated. In this case, the constitutive equations can be cast in terms of a variational statement similar to von Mises’s principle of maximum plastic dissipation. Moreover, if the yield function is convex, the governing equations may be cast in terms of a convex mathematical program. Such a formulation allows for a straightforward analysis of properties related to the existence and uniqueness of solutions. In the general nonassociated case, such a formulation is not possible. Although this shortcoming does not pose any fundamental obstacles to developing conventional Newton-Raphson based solution methods analogous to those of the associated case, the desirable mathematical properties of associated plasticity are lost. Moreover, whereas associated plasticity involves a symmetric tangent modulus, a nonassociated flow rule generally gives rise of an unsymmetric set of discrete finite element equations. Finally, although problems of associated plasticity can be solved very efficiently using methods of modern mathematical programming, such formulations are not possible in the nonassociated case.
Motivated by the relative efficiency and robustness of numerical algorithms for associated plasticity, a numerical formulation that retains the desirable properties of associated computational plasticity, but which is applicable to general nonassociated models, is presented in the following. It should be noted, however, that the characteristics of nonassociated plasticity in terms of localization of deformations are retained. This includes the apparent global softening often observed in boundary value problems involving a perfectly plastic nonassociated model. Similarly, the non-uniqueness of solutions implied by most nonassociated models also persists and is manifested by a strong sensitivity to the finite element mesh, boundary conditions, and so on.

9.2.2 Micromechanics of friction

The basic idea behind the formulation derives from the structure of the internal dissipation associated with constitutive models of the type presented in Sections 1. Let us assume a yield function of the type

where and are some measure of deviatoric and mean stress respectively, is a friction coefficient and is the cohesion. Furthermore, let the plastic potential be given by

where is a dilation coefficient.
In space, the plastic strain rates are given by

where and are the volumetric and deviatoric strains conjugate to and respectively. The internal dissipation associated with the model Eq. 80 is then given by

Figure 16: Microscopic origins of friction as plastic shearing of asperities. A higher confining pressure implies a higher degree of interlocking of the asperities and thereby a higher apparent shear strength.

This expression for the internal dissipation reveals several interesting, albeit well known, features. Firstly, for an associated material (), the dissipation is proportional to the internal cohesion, . As such, no internal dissipation takes place for a purely frictional material () which is in obvious contrast to experimental evidence. Secondly, for , the dissipation is proportional to an apparent cohesion which comprises two terms: the internal cohesion and a contribution which stems from the prescribed degree of nonassociativity. The interpretation of the latter term as an apparent cohesion is consistent with the viewpoint that friction results from the mechanical interaction of microscopic asperities on the surfaces of the solids in contact . With the stresses at the scale of the asperities being much greater than the elastic limit of the material, it is primarily plastic deformations at the microscale that govern the macroscopically observed frictional resistance. This point is illustrated in Figure 16 which shows two rough surfaces under different levels of confining pressure. The plastic shearing may be assumed to be of the ductile, purely cohesive kind. For brittle materials such as sand grains, this assumption is justified by the very high stress level at the scale of the asperities which effectively renders the otherwise brittle material ductile. The apparent shear strength of each assembly thus derives exclusively from the geometric changes induced by the confining pressure and with Coulomb friction implying a linear relation between apparent shear strength and confining pressure. This interpretation motivates rewriting the yield function Eq. 78 as


is the apparent, pressure dependent, cohesion. This material parameter embodies all the complexities of the actual micromechanics of the frictional interfaces and their evolution in response to the applied loads. As such, its relative simplicity is surprising, but nevertheless found to be appropriate for a very broad range of materials although there are also a number of noteworthy exceptions as discussed for example by .

9.2.3 Time discrete formulation

Suppose now that the apparent cohesion, , is known. The associated flow rule then produces the desired result, namely the plastic strain rates Eq. 80. In the solution of boundary value problems, the apparent cohesion is of course not known a priori as it is directly proportional to the pressure that is to be determined as part of the solution. However, assuming that such problems are solved incrementally via a sequence of pseudo-time steps, some parts of the yield function may, in principle, be evaluated implicitly while other parts may be evaluated explicitly. Assume that the state at time is known. The yield condition imposed at may then be approximated as


Again, the associated flow rule produces the desired time discrete result, namely:

However, the explicit evaluation of the apparent cohesion means that the original yield function may be exceeded for the new stress state at . Similarly, the approximation may imply plastic yielding for stress states that would otherwise be deemed purely elastic (see Figure 17). However, for small enough increments, i.e. for , it would be expected that the error introduced by the explicit evaluation of the apparent cohesion would tend to vanish. This supposition is confirmed by numerical experiments .
Any problem of nonassociated plasticity can thus be transformed into an equivalent approximate problem of associated plasticity. Effectively, each time step involves the solution of an approximate associated problem. As such, methods previously developed for associated plasticity are applicable with little modification. Moreover, variational formulations that subsequently are resolved using either general or more specialized methods of mathematical programming are applicable. Such methods are adopted in OPTUM G2.

Figure 17: Explicit evaluation of apparent cohesion: original and approximate yield functions. The area indicated by (a) is non-permissible according to the original yield function but permissible according to the approximate yield function. Similarly, the area indicated by (b) is non-permissible according to the approximate yield function but is within the elastic domain according to the original yield function.

9.2.4 Mohr-Coulomb

With the above considerations, the plane strain Mohr-Coulomb yield function actually implemented takes the following form:


with subscript referring to the current, known, state. In other words, only the regular cone is modified whereas the tension and compression cut-offs remain unaltered.

9.2.5 Drucker-Prager

Similarly, the Drucker-Prager yield function is given by


with subscript referring to the current, known, state. Again, it is only the regular cone that requires modification while the tension and compression cut-offs remain unaltered.

9.3 Status of limit analysis

Classic limit analysis relies crucially on the concept of associated flow. Indeed, as discussed in Section 2, the upper and lower bound theorems are only valid under the assumption that the flow rule is associated. For nonassociated flow, the only available result is that an upper bound for a material with is also an upper bound for the corresponding material with . In other words, the limit load calculated on the basis of an associated flow rule is on the unsafe side. How much on the unsafe side depends both on the degree of nonassociativity (the difference between and ) and on the particular problem. As described in Section 4, it is possible to derive an ‘effective strength domain’, defined by the effective material parameters and (for Mohr-Coulomb) or and (for Drucker-Prager).
It has long been proposed that these effective parameters can be used in a standard associated plasticity framework to compute ‘nonassociated limit loads’ . For the simple block example of Section 4, the effective material parameters would certainly result in the correct solution. For general problems, solved for example by means of the finite element method, the use of the effective material parameters in an associated framework also provide reasonable estimates of the magnitude of the ultimate limit load resulting from a full elastoplastic calculation with the correct friction and dilation angles. However, the non-uniqueness implied by the introduction of a nonassociated flow rule should be borne in mind – as should the fact that the exact value that results from a finite element analysis will be highly dependent on the geometry, boundary conditions, imperfections, as well as the nature of the mesh in terms of its ability to capture the expected localized modes of deformation. This is illustrated in the example below.

Figure 18: Uniaxial compression test setup (left) and possible mode of localization (right).
Figure 19: Schematic load-displacement behaviour of uniaxial compression test.
Figure 20: Determinant of acoustic tensor (normalized by determinant of elastic modulus) as function of shear band inclination angle for associated and nonassociated flow rules at point A in Figure 19.
Figure 21: Limit load for nonassociated material as function of shear band inclination angle according to Eq. 84. The associated limit load is indicated by the dashed line.

We consider a prismatic sample of a cohesive-frictional material as shown in Figure 18. The material model is Mohr-Coulomb with elastic parameters and and cohesion and friction angle both greater than zero. The top platten, which is assumed smooth, is gradually displaced downwards while the normal and shear stresses on the vertical sides are maintained at zero. Regardless of the flow rule, i.e. of the value of , the load-displacement behaviour up to the yield point (A) will be as shown in Figure 19. At this point, the vertical stress is:

corresponding to a continuing state of homogeneous deformation, again independent of the flow rule (although the homogeneous plastic strain field of course will depend on the flow rule). For the associated flow rule this is also the exact (and unique) limit load.
To gain insights into possible localized modes of deformation, the acoustic tensor is examined for the stress state at point A. The result is shown in Figure 20 for an associated flow rule with and for a nonassociated flow rule with and . As seen, the associated flow rule implies the possibility of a shear band inclined at corresponding to the well-known slip line solution that results in the collapse load Eq. 83.
textThe nonassociated flow rule, on the other hand, allows for the possibility of localizing along a range of directions. It may be shown that this range is . For the angles in question, this means that any localized solution defined by a shear band inclination in the interval is possible. Following the analysis of Section 1, it may be shown that the corresponding collapse load is given by

where and are given by Eq. 76. The minimum values is attained for leading to a collapse load of

For the entire range of , the nonassociated limit load Eq. 84 is less than the associated limit load Eq. 83 while for the two coincide (see Figure 21). This is also indicated in Figure 19, where the strain magnitude in principle is zero, but in finite element calculations will depend on the thickness of the shear band which is limited from below by the resolution of the mesh.
Interestingly, Figure 21 also reveals that the limit load in the nonassociated case is below that of the associated case outside the range of possible shear band inclinations, i.e. for . However, the governing equations do not allow for a switch from a homogeneous and a localized mode of deformation in that range.
In finite element calculations, the above analytical observations are observed by firstly a reduction in limit load and secondly a high degree of sensitivity to imperfections. An example is shown in Figures 22-23.

Figure 22: Load-displacement curves for uniaxial compression test.
Figure 23: Finite element mesh (six-node elements) and different localized modes of deformation for uniaxial compression test at (see Figure 22).

Figure 22 here shows the load-displacement curves of various simulations of the uniaxial compression test with . As expected the associated calculation leads to the analytical limit load Eq. 83. Next, the dilation angle is set to and the coordinates of the finite element mesh are perturbed randomly by a magnitude of times the width of the sample. This leads to a family of load-displacement curves that, while following the expected trend of displaying a reduction in strength, are somewhat different, both in terms of their evolution after the peak as well as in terms of the final residual load. This illustrates the non-uniqueness discussed above.
Regarding the magnitude of the limit load, it can be seen that it is somewhat above the minimum magnitude predicted theoretically Eq. 85. This is due to an inability of the mesh to capture the minimum solution defined by a shear band inclination of . Indeed, the localized modes of deformation shown in Figure 22 correspond to and thus, according to Eq. 84, to a limit load of .
The trends observed in this example are quite representative of typical problems of practical interest . In particular, concerning the applicability of limit analysis, the use of the original and lead to an overestimate of the bearing capacity while the use of the effective parameters and leads to an underestimate. While the very nature of nonassociated plasticity in principle makes it meaningless to speak of the limit load, the effective parameters often provide a reasonable estimate to what is obtained in a full elastoplastic analysis using the correct , and .

9.4 Undrained conditions

While the dilation angle under drained conditions may have some influence of the bearing capacity, its effect under undrained conditions is much more dramatic. Indeed, the Mohr-Coulomb, Drucker-Prager and similar models under undrained conditions must be used with a dilation angle or coefficient equal to zero to obtain reasonable results. Indeed, any finite dilation will lead to an infinite bearing capacity. This result is sometimes presented as a numerical artefact, but is in fact in complete accordance with these models and there is nothing suspect or mysterious about it.
To see this, consider the relation between excess pore pressure and volumetric strain in incremental form:

Clearly, an increase in volumetric strain, , leads to an increase in excess pore water pressure, . For a fixed total stress, this in turn leads to a decrease in effective mean stress:

For frictional yield criteria, this effectively means that a material point undergoes an apparent ‘hardening’ once yielding commences and dilation results (see Figure 24). Consequently, an ultimate limit state is never attained. A nonassociated flow rule implying zero volumetric straining must therefore be used. Thus, with the Mohr-Coulomb model it is necessary to use a nonassociated flow rule with . Similarly, with the Drucker-Prager model, must be used.

Figure 24: Consequence of excess pore generation as a result of plastic dilation.



Many problems in geotechnics involve the flow of water through the soil. Typical problem include the flow through earth dams and in and out of excavations. Seepage problems may be divided into two categories: confined and unconfined. In former case, all boundary conditions are known a priori while in the latter case some of the boundary conditions are to be found as part of the solution. This situation is sketched in Figure 25. The pressures on the submerged boundaries are here known as is the no-flow condition at the base of the dam. But the location of the free surface is not known a priori and is to be determined as part of the solution. This complicates matters considerably and an iterative procedure must generally be used.

Figure 25: Unconfined seepage problem.

Alternatively, a more general saturated/unsaturated approach may be taken whereby the governing equations are valid throughout the domain, but involve material parameters that vary significantly with the degree of saturation. In the following, this is the approach taken.

10.1 Governing equations

Variably saturated flow through porous media can be described by the mass balance equation

supplemented with the generalized Darcy’s law


  • = Porosity
  • = Degree of saturation
  • = Fluid velocity [m/day]
  • = Saturated hydraulic conductivity modulus [m/day]
  • = Relative hydraulic conductivity which is a function of degree of saturation
  • = Vertical coordinate
  • = Unit weight of water (= 9.8 kN/m)
  • = Pressure
  • = Head (= )

The effective hydraulic conductivity is the product of the relative conductivity, , and the saturated conductivity, :

where and are the saturated hydraulic conductivities in the and directions respectively.
Combining Eq. 86 and Eq. 87 leads to what is sometimes (especially in 1D) called Richards equation:

For the time being, we are only concerned with steady state solutions so that the above governing equation simplifies to

It should be noted that both the relative hydraulic conductivity, , and the degree of saturation, , are highly nonlinear functions of the pressure . For further details and relevant relationships, please refer to the Materials Manual.

10.2 Boundary conditions

Boundary conditions come in the form of standard Dirichlet conditions:

and Neuman conditions:

where , and are prescribed pressures, heads and fluxes respectively, and is the outward normal to the boundary in question. With references to Figure 25, the Dirichlet condition is imposed on boundaries AB and FG while the Neuman condition is imposed along AG (). At the free surface BE (the exact location of which unknown), both the pressure and the normal flux are zero.
Furthermore, along the seepage face EF, the pressure is zero while, at the same time, the normal flux is required to be non-negative. Such segments are referred to as seepage faces. In OPTUM G2, all external boundaries are seepage faces by default and remain so unless another boundary condition is imposed.


The most basic consolidation problem is that of a fine grained fully saturated soil being subjected to a rapidly applied load that is then kept constant. This will lead to some immediate settlements and to the generation of excess pore pressures. With time, these pressures will dissipate and further settlements take place. The situation is sketched in Figure 3.
In simples cases, the classic theory of Terzaghi may used while, in the more general case, the theory of Biot is applicable. The latter, which is the basis of the Consolidation analysis in OPTUM G2, contains the former as a special case.

11.1 Governing equations

The governing equations of elastoplastic consolidation (assuming linear elasticity/perfect plasticity) are as follows.

Strain-displacement relations:

Hooke’s law:

Flow rule:

Yield and complementarity conditions:

Pore fluid conservation:

Darcy’s law:

where are the fluid velocities, is the permeability matrix and is the vertical coordinate. After replacing the volumetric strain rate by the finite difference approximation , these governing equations may be cast in terms of the following variational principle:

Assuming that the seepage pressures and their associated boundary fluxes are available, this simplifies to:

which generalizes the Hellinger-Reissner type principles discussed in previous sections. The Euler-Lagrange equations are:

where it has been assumed that the following boundary conditions are satisfied:

Following the procedure described in detail in previous sections, the above principle is discretized by replacing the continuous variables (, and ) by their discrete counterparts using finite element shape functions. Furthermore, it is assumed that the steady state seepage pressures are available. The following discrete principle then appears:


where a new set of variables, , has been introduced. It is noted that the above formulation reproduces the usual undrained problem for .
Regarding the exact discretization, it is prudent, for stability issues, to operate with a low-order approximation of the excess pore pressures. In OPTUM G2, for the analysis type Consolidation, a linear and continuous interpolation is always used regardless of the order of the stress and displacement interpolations. Similarly, the auxiliary variables are always assumed constant within each element. An example of such an element is shown in Figure 26.
Finally, it is noted that nonassociated flow rules, nonlinear elastic laws, hardening, etc can be used in OPTUM G2. The implementation of these features follows that described in the previous sections.

Figure 26: Mixed consolidation element (6-node Gauss).



It is commonly recognized that typical soil parameters of relevance to design display a not insignificant spatial variation. In traditional design this is accounted for by using spatially independent “cautious estimates" of the relevant parameters. For example, if a given parameter has been measured at a number of different locations, one may decide on the final value by using a weighted average skewed towards the unfavourable end of the measured parameter range.
As an alternative to traditional deterministic analysis, OPTUM G2 offers the possibility to conduct stochastic analysis. The spatial variability is here included in the analysis such that not only one value of bearing capacity, settlement, strength reduction factor is obtained, etc is obtained, but rather a probability distribution of these.
In OPTUM G2, spatial variation of the parameters can be taken into account by generating random fields for specific parameters and subsequently running the analysis via a series of Monte-Carlo simulations. The end result is probability distributions of the key results of the analysis, for example bearing capacity (Limit Analysis), settlement (Elastoplastic analysis) or strength reduction factor (Strength Reduction analysis).

12.1 OPTUM G2 input

To generate a random field of a given material parameter, four input parameters are required:

  • The mean value of the parameter (30 kPa in Figure 27). Gradient, Profile and Map can be used to specify the spatial distribution of the mean value of the parameter.
  • The coefficient of variation of the parameter (35% in Figure 27).
  • The horizontal correlation length of the parameter (50 m in Figure 27).
  • The vertical correlation length (1 m in Figure 27).
Figure 27: Specification of parameters for stochastic analysis. The mean value (or spatial distribution of mean value) is specified under the Space tab (middle) and the remaining parameters are specified under the Random tab (right).

The parameters are specified with respect to an underlying probability density function (Normal or Lognormal).

Figure 28: Examples of random fields of undrained shear strength. In all cases, the mean value and coefficient of variation are kPa and respectively.

The spatial correlation length describes the distance over which the spatially random values will tend to be significantly correlated. A large value will thus imply a smoothly varying field while a smaller value will imply a more ragged field. Some examples are shown in Figure 28. It is usually assumed that the spatial correlation structure of all engineering properties of rocks and soils are identical. That is, the correlation structure of, say, the friction angle is assumed to be identical to that of the Young’s modulus. It is commonly recognized that the horizontal correlation length is rather larger than the vertical correlation length. conducted a literature review and concluded that the horizontal correlation length is of order 40-60 m while the vertical correlation length is in the range of 2-6 m.
Other than the physical parameters, further parameters related to stochastic analysis are available in Project under Stochastic Parameters. These include:

  • The number of Monte Carlo runs (default = 1,000).
  • The seed used to generate the random field. The seed increases by 1 for each new Monte Carlo run. Thus, in order to recover the problem of the ’th run in a Monte Carlo simulation, the seed can be set to and a single Monte Carlo run performed.

12.2 OPTUM G2 output

The output from analyses with random parameters are the probability density and cumulative density functions for various quantities of interest depending on the analysis type. For Limit Analysis the collapse multiplier is recorded, for Strength Reduction analysis the strength reduction factor is recorded while for Elastoplastic analysis key data such as displacements at Result Points are recorded. Additionally, as described above, to recover the full solution of the ’th Monte Carlo run, the seed (Project/Stochastic Parameters) should be set equal to and a single Monte Carlo run performed.
If more information than what is available by default is required, the command line version of the program, OPTUM G2Cmd, is useful. A desired number of input files can first be generated, each with a different seed and set to perform a single Monte Carlo run. These can then be run using OPTUM G2Cmd and the results further processed via an external script.

12.3 Random field theory

The following section gives a brief overview of the concept of random fields and the generation of these fields.
Let be a random field, where defines the physical space and defines the probability space. The correlation structure of a random field is modeled by the covariance function, denoted by , where , are bounded, symmetric and positively defined.
The first issue in the application of random field theory is the selection of a probability distribution for the field. Although many probability distributions can be used, the lognormal distribution has some distinct advantages and is commonly used to model the engineering properties of soils and rocks. The lognormal distribution offers the advantage of simplicity in that it is arrived at by a simple nonlinear transformation of the classic normal Gaussian distribution. Furthermore, the lognormal distribution guarantees that the random variable is always positive.
Secondly, regarding the covariance function, there are several possibilities with the exponential covariance function often being the function of choice. Besides the mean and standard deviation, the covariance function also involves the spatial correlation length.

12.3.1 Karhunen-Loeve expansion

Several random field generation methods are available . In OPTUM G2, the Karhunen-Loeve expansion method is used. This method is convenient as it provides analytical solutions for the exponential covariance function. Using Mercer’s theorem, the covariance function can be decomposed according to

where and are, respectively, the eigenvalues and eigenfunction of the .
Since the above sum has to be truncated to a finite number of terms, a significant concern is that the simulated variance will be reduced. In order to control this reduction, the eigenvalues are sorted in descending order and the number of terms, , is decided on by the eigenvalues having decayed sufficiently to satisfy the condition

where typically is set to .

12.4 Risk

Most stochastic analyses of the kind that are possible to carry out in OPTUM G2 are aimed at determining the probability of failure, or of the settlement exceeding a certain magnitude, or similar. Moreover, the full probability density function is also determined, thus providing valuable information on the sensitivity of the problem.
The concept of risk takes this type of analysis one step further by not only considering the probability of a certain event (e.g. failure) but also quantifying the consequences of such events. While the consequence of a certain event obviously is highly problem dependent, a general measure that often would be of interest in connection with failure is the volume of the soil mobilized. This is particularly relevant for slope stability where a deep seated failure usually would imply a higher consequence than a more shallow one. Consequently, in OPTUM G2, the volume of soil mobilized in each Monte Carlo run is recorded. Following , the risk, , may then be defined as

is the volume of mobilized soil for the ’th – assuming that this run implies failure, otherwise .
In OPTUM G2, slope stability analysis would usually be carried out using either Strength Reduction analysis or Limit Analysis with Multiplier = Gravity. In both cases, strength reduction factor/collapse multiplier less than 1 implies failure while a value above 1 implies that the system is stable. Since the volume of mobilized soil is recorded for all runs, a manual filtering is necessary for the risk to be calculated according to the above.


In this section the solid elements available in OPTUM G2 are detailed. These include elements that lead to rigorous upper and lower bounds on the exact solution as well as ‘mixed’ elements that often are more accurate but which do not lead to rigorously bounded solutions.

13.1 Lower bound element

Consider the elastoplastic lower bound principle:

As shown in Section 3, any stress field that satisfies the constraints of this problem leads to a lower bound estimate of the objective function – which in the case of rigid plasticity () is the collapse multiplier. With the equilibrium equations involving first derivatives and considering that the yield function must be satisfied everywhere, the obvious candidate for a lower bound element is a triangle with a linear variation of the stresses between the corner nodes. This element, shown in Figure 30, is capable of satisfying all the above constraints provided that unit weight, , is constant over the element, that the seepage pressures, , vary linearly within each element, and that the yield conditions are enforced at the three corner nodes.

Figure 30: Lower bound element.

We introduce the following finite element approximations for the stresses and the excess pore pressures:

where , , , and are the nodal effective stresses, hardening variables, steady state pore pressures and excess pore pressures respectively and , , and contain the respective linear shape functions. Typically, the hardening variables and excess pore pressures are interpolated in the same fashion as the stresses, i.e. via linear shape functions.
The discrete equilibrium equations may then be written as




Figure 31: Lower bound elements joined by two zero-thickness elements to produce a statically admissible stress discontinuity.

The above equilibrium equation can here be multiplied by the element area, , to obtain a set of equations that are valid even for a zero-thickness element. This property can be utilized to construct statically admissible stress discontinuities between elements. The typical situation is shown in Figure 31 where two ‘continuum elements’ are connected by a patch of zero-thickness elements of exactly the same type. The four equilibrium equations established in the zero-thickness patch effectively impose continuity in the normal and shear stresses at each end of the interface and thus along the entire interface.
Evaluating the quadratic term in the objective function of Eq. 94 by standard Gauss quadrature such that and , the final discrete optimization problem to be solved reads:

where the yield conditions are enforced at each of the stress points.
The above formulation is a generalization of the standard lower bound element originally conceived for limit analysis applications .

13.2 Mixed elements

Consider the elastoplastic Hellinger-Reissner principle:

The discretization of this principle proceeds by specifying approximations for the variables involved, namely the effective stresses, pore pressures and displacements. Using standard finite element terminology, these are given by

Although there in principle is some freedom in choosing the various shape functions, we will make the following assumptions: i) the displacements are continuous between elements, ii) the stresses, hardening variables, and excess pore pressures are approximated in the same manner and are discontinuous between elements, and iii) the polynomial degree of the former is one order higher than that of the latter. With these assumptions, standard finite elements originating from direct application of the principle of virtual work and involving only the displacements as unknowns can be reproduced as a special case.
Inserting the above approximations into Eq. 96 leads to the following discrete problem:


The problem Eq. 97 may be solved with respect to the incremental displacements to yield the final problem:

The above formulation is similar to the one proposed by and is a generalization of earlier formulations by .

13.2.1 Gaussian and Lagrangian families

As mentioned above, the stresses are interpolated using shape functions that are one polynomial order lower than those used for the displacements. Moreover, the displacement shape functions are always Lagrange polynomials with the nodes located at the so-called Lagrange points. The location of the stress points, on the other hand, is less obvious. One possibility is the Gauss points while another is the Lagrange points. In OPTUM G2, each of these possibilities are considered and the corresponding elements are labeled accordingly (Gauss or Lagrange).
An example is shown in Figure 32. The displacements here vary quadratically while the stresses vary linearly. For the Gauss elements, the stress points are the Gauss points with Barycentric coordinates while the Lagrange element uses the corner nodes as interpolation points.

Figure 32: Quadratic displacement/linear stress triangles of the Gauss and Lagrange types.

13.2.2 Upper bound elements

Upper bound elements may be constructed directly on the basis of upper bound principles such as Eq. 31. Alternatively, they may be be derived as special cases of mixed elements. This latter approach is in many ways preferable as it avoids any explicit use of the dissipation function (which may not be straightforward to derive although it in principle can be done once the yield function is available).
Two upper bound elements are possible: a linear displacement/constant stress triangle and a quadratic displacement/linear stress triangle first proposed by . The former appears as the only mixed element possible in this case. The latter appears as a special case of a mixed element where the stresses are interpolated on the basis of the Lagrange points and the Lagrange points are used as integration points, with equal weights, for the matrices and in Eq. 98. The formulation of upper bound elements as special cases of mixed elements is discussed by .

Figure 33: Selection of elements available in OPTUM G2.

13.2.3 Overview of elements

An overview of some of the elements available in OPTUM G2 is given in Figure 33.

13.2.4 Discontinuities

Discontinuities may be introduced by collapsing patches of regular continuum elements to zero thickness . Such discontinuities are included by default between all solid domains.  


The association between optimization and computational limit analysis is long established. Indeed, all computational limit analysis procedures lead directly to optimization problems. In contrast, the application of optimization methods as a means of solving elastoplastic problems is more unusual. However, the idea of applying optimization algorithms to such problems is an old one going back at least to the work of Giulio Maier and his coworkers in the 1960s .
The application of the finite element method to nonlinear constitutive models involving nonlinear elasticity, plasticity, hardening, etc generates systems of nonlinear equations. Typically, these must be solved in the course of a number of load steps (as in static elastoplastic analysis) or time steps (as in consolidation analysis). The most common approach to the solution of nonlinear finite element equations is the classic Newton-Raphson method. The nonlinear equations are here linearized and a problem that resembles one of linear elasticity is solved in a sequence of iterations until convergence, i.e. until the original nonlinear equations are satisfied to within some specified tolerance.
In OPTUM G2 a different approach is used. Instead of solving the nonlinear equations directly, an equivalent optimization problem is solved instead, the solution to which satisfies the original equations. Consider as an example, the nonlinear equation:

Application of the Newton-Raphson method to this equation is straightforward, eventually leading to the solution .
The Optum approach would instead consider the optimization problem:

To solve this problem, we differentiate the objective function (the function being minimized) and set the result equal to zero:

which recovers the original nonlinear equation Eq. 100. We can then proceed with any method of choice, including Newton-Raphson, to solve the resulting equation – in optimization terminology, the optimality condition.
For this particular problem, the two approaches are completely identical. For more complex problems, however, the optimization approach is usually technically superior to the traditional one in terms of efficiency and robustness. These include both problems of elastoplasticity and of limit analysis where the optimization approach is the only feasible one.

14.1 Algorithms

The earliest algorithms applied to limit analysis were designed for linear programming (LP). That is, algorithms capable of solving linear programming problems of the type:

Both upper and lower bound limit analysis can be cast in this form. In lower bound limit analysis the objective function would comprise the load multiplier, the equality constraints would represent the equilibrium equations and the inequality constraints would account for the yield condition. Since most practical yield conditions are nonlinear, a linearization has to be performed before the problem can be set up as a linear program. While such linearizations (internal or external to the original yield surface as desired) are straightforward to formulate, early LP algorithms suffered from the trait that their running time were exponentially dependent on the number of constraints. It was therefore hailed as a major breakthrough when announced the discovery of a polynomial time algorithm in 1984 (not just in academic circles but also in the mainstream media, see Figure 34). While there was significant excitement, even hype, around the new algorithm and its potential for commercial exploitation, it was quickly discovered that it in fact was closely related to the so-called Sequential Unconstrain Minimization Techniques proposed earlier by . A notable feature of these techniques is that they are equally applicable to linear and nonlinear programming whereas earlier methods distinguished sharply between the two. These insights came to define much of the research in the field in the decades that followed. Labeled interior-point methods, the class of algorithms of which Karmarkar’s is one is commonly viewed as having revolutionized the field of numerical optimization .
Interior-point methods are applicable to problems that can be written as

where and usually are assumed convex.
In a further development, realizing that much could be gained in terms of devising efficient and robust algorithms for more restricted classes of problems than the most general convex one, algorithms for so-called conic programming were proposed . Conic programs are often written in the following standard form:

where is a cone. Examples of cones include the quadratic cone:

and the positive semidefinite cone in which a square matrix of variables is required to be positive semidefinite

Both these cones are of interest in plasticity – the first for the modeling of quadratic yield criteria such as Drucker-Prager while second can be used to model the Mohr-Coulomb criterion.

Figure 34: News clippings related to Karmarkar’s algorithm (collected by RJ Vanderbei,

A further specialization of general conic programming is second-order cone programming (SOCP) which considers only quadratic constraints, i.e. the quadratic cone Eq. 104 and the rotated quadratic cone:

The algorithms SONIC/Serial and SONIC/Parallel included in OPTUM G2 are SOCP algorithms. A partial list of papers dealing with the application of SOCP (as well as general conic programming) to plasticity problems is as follows:

  •  : plane strain Mohr-Coulomb limit analysis using SOCP.
  •  : limit analysis and elastoplasticity for plane strain Mohr-Coulomb problems, limit analysis and optimal design of reinforced concrete plates using SOCP.
  •  : 3D Mohr-Coulomb limit analysis using semidefinite programming.
  •  : drained and undrained elastoplastic analysis of plane strain Modified Cam Clay problems using SOCP.
  •  : nonassociated elastoplastic analysis using SOCP.
  •  : dynamic large deformation analysis and granular flows governed by Mohr-Coulomb using SOCP.

14.2 Approximation of non-quadratic problems

OPTUM G2 makes use of SOCP only. This poses somewhat of a challenge as many of the problems dealt with are not quadratic in nature. The Hoek-Brown criterion, for example, is not quadratic. Similarly, both the elastic law and the hardening law of Modified Cam Clay are not quadratic. However, it is relatively straightforward to approximate these problems in terms of quadratic problems and then solve the resulting approximate problems in a sequence of fixed-point iterations.
To illustrate the typical procedure, consider again the problem:

The second term is here non-quadratic and cannot be dealt with directly using SOCP. Instead, consider the approximation given by the second-order Taylor expansion about a point :

We then have

which can be solved using SOCP. Starting with , we obtain the solution

Next, setting and solving again gives the solution

Using this solution for gives

and so on until convergence to the correct solution of .


A central question faced by user of finite element software is to what extent the results can be taken at face value – to what extent do they represent reality, to what extent are well-known analytical solutions reproduced, etc. The basic idea behind the concepts of Verification and Validation (V&V) is to aid in answering these questions in a systematic manner. V&V has been used extensively in product design and software engineering and have in recent years also found its way into the mechanics of materials, computational mechanics, and related fields.

Figure 35: From question to decision using numerical analysis.

The use of numerical analysis begins by posing a relevant question. For example: what is the bearing capacity of a given foundation, what should the dimensions of a sheet pile wall supporting a given excavation be, etc. The process of a arriving at an answer to such questions, and thereby make the relevant decisions, is usually approached via the sequence of steps sketched in Figure 35.
The first task consists of establishing a conceptual model, i.e. an abstraction and eventually simplification of physical reality. In the conceptual model, only the subset of physical phenomena that is thought to have an influence on the problem that we eventually want to describe is included. In addition, the establishment of the conceptual model involves collecting and organizing relevant observations, for example on the response of the material in laboratory tests, site geometry, ground water conditions, etc.
Secondly, on the basis of observations and possibly more ‘fundamental’ theories, a mathematical model is formulated. By a mathematical model is meant a set of equations that, as a minimum, captures the available observations and that further is expected to make useful predictions outside the region of experience. For solid mechanics problems, these would typically include equilibrium equations, strain-displacement relations and stress-strain relations.
The third task is to solve the governing equations. This may be done by approximating the continuous variables (displacements, stresses, etc) of the original governing equations in terms of discrete counterparts. The finite element method, where the continuous variables are approximated on the basis of element shape functions and the discrete nodal variables, offers one way of discretizing the governing equations. The resulting set of discrete equations are referred to as the computational model and must in practice be solved numerically.
Finally, by solving the discrete governing equations, we arrive at a solution, i.e. a set of stresses, displacements, etc, on the basis of which further decisions can be made.
The idea behind V&V is to gauge the usefulness of this solution in a systematic manner. Though often confused, verification and validation are two independent procedures. In the present context, we may say that the two procedures involve the following questions:

  • Validation: is the mathematical model suitable in terms of capturing the physical reality?
  • Verification: are the equations that comprise the mathematical model solved correctly?


  • Validation: are the correct equations solved?
  • Verification: are the equations solved correctly?

It is important to deal with these two questions separately and, more generally, to recognize that the final solution is influenced by errors at several different stages in the process of arriving at the numerical solution.

15.1 Verification

The process of verification and validation begins from the end, i.e. with the computed solution. This must first and foremost be verified as actually satisfying the discrete governing equations. Although this task is relatively straightforward (the solution is simply inserted into the discrete equations), the task of ensuring that the correct solution is computed is not. Indeed, the discrete equations associated with typical geomechanics problems tend to be highly nonlinear, will usually involve certain types of singularities, and are often ill-posed. These difficulties are in addition to standard numerical problems such as truncation and round-off errors. It is therefore not uncommon that a significant error results from the process of solving the discrete governing equations alone. And although most finite element programs – including OPTUM G2 – provide information about the final solution accuracy and have in-built procedures for detecting errors of an unacceptable magnitude, it is ultimately the responsibility of the user to ensure that the computed solution satisfies the discrete governing equations to within an acceptable tolerance.
A second, and usually more common and serious error is that which results from the discretization of the governing equations. Using standard finite element methods, the only real information about the inherent error is that the computed solution is on the unsafe side in the sense that the bearing capacity is overestimated and the settlements are underestimated. No robust and general means of estimating the magnitude of the error is available. Typically, one therefore carries out a convergence analysis where the mesh is gradually refined until it is observed that some key quantities (displacements, reactions, etc) attain steady values. Although this procedure in many cases is quite appropriate, it has certain limitations, the most important of which is that the number of finite elements required to attain convergence may be prohibitively large. Similarly, care must be taken that convergence really is attained, i.e. that the rate of change from one mesh to the next really is sufficiently small for the solution to be deemed converged.
As an alternative to the this approach, it is in some cases possible to construct finite element approximations where it can be guaranteed that certain quantities are either underestimated or overestimated. Such elements are implemented in OPTUM G2 and may be used to compute upper and lower bounds on the bearing capacity or factor of safety of geostructures. Similarly, bounds on the elastic energy may be computed such that solutions that are either ‘too stiff’ or ‘too flexible’ can be obtained. As such, using OPTUM G2, the verification of the computational model against the mathematical models consists of refining the mesh until the gap between the upper and lower bound solutions is deemed sufficiently small. Denoting the lower bound solution by and the upper bound solution by , the exact solution, , is bounded by

Alternatively, the exact solution can be approximated by the mean value between the upper and lower bounds:

We then have

where and are the absolute and relative errors, respectively, which are bounded by


As an example, consider a strip footing of width subjected to a central vertical load and resting on a purely frictional material with unit weight and friction angle . For this problem, the exact bearing capacity is . Suppose that a lower bound of and an upper bound of have been computed. The mean value between the upper and lower bounds is then

giving a relative error of

In other words, the exact solution is given by

As indicated by this example, the actual error in the mean value between the upper and lower bounds is much smaller than suggested by the worst-case scenario. Indeed, with a mean value of and an exact solution of , the actual error is only about 1.7%). In practice, the upper and lower bounds tend to converge at the same rate to the exact solution, i.e. , such that often is a very good estimate of the exact solution, even for very coarse meshes.

15.2 Validation

The second part of the V&V procedure aims to first of all validate the mathematical model against the conceptual model. For this purpose, one usually has to rely on incomplete observations such as the response of the representative soil samples in laboratory experiments and possibly measurements of deformations, pore pressure changes, etc during the construction phase. In addition, certain types of expected response should be validated. For example, the surface settlements at some distance from a vertically loaded foundation should tend to zero. Some of the more common model errors made in translating the conceptual model into a mathematical model include:

  • Shortcomings in the material model
  • Inappropriate initial stress distributions
  • Inappropriate initial pore pressure distributions

Secondly, the conceptual model itself must be scrutinized and held up against physical reality. Relevant questions include:

  • Have certain observations been misinterpreted?
  • Is the conceptual model too simple?
  • Have all relevant phenomena been considered?

These and related questions are often not easily answered unequivocally and will inevitably involve a substantial amount of judgment that is not easily quantified.

15.3 Proper use of VV

The key feature of the Verification and Validation framework described above is that it facilitates a rigorous analysis of the possible errors generated at each step on the way from physical reality to numerical solution. In particular, it stresses that it is not the final numerical solution that should be compared directly with physical reality. Indeed, it may well be that errors made at the different stages counteract each other. For example, using a substandard finite element discretization in conjunction with an equally inappropriate constitutive model might well lead to a ‘better’ result than if a higher quality finite element discretization was employed with the same constitutive model. Proper use of V&V would identify this scenario as being fundamentally unsatisfactory as a given solution can only be deemed acceptable if the errors involved in each of the stages towards obtaining it are acceptable.  

Back to OPTUM G2 main page