Application of discontinuity layout optimization to plane plasticity problems

Colin Smith, Matthew Gilbert


A new and potentially widely applicable numerical analysis procedure for continuum mechanics problems is described. The procedure is used here to determine the critical layout of discontinuities and associated upper-bound limit load for plane plasticity problems. Potential discontinuities, which interlink nodes laid out over the body under consideration, are permitted to crossover one another giving a much wider search space than when such discontinuities are located only at the edges of finite elements of fixed topology. Highly efficient linear programming solvers can be employed when certain popular failure criteria are specified (e.g. Tresca or Mohr–Coulomb in plane strain). Stress/velocity singularities are automatically identified and visual interpretation of the output is straightforward. The procedure, coined ‘discontinuity layout optimization’ (DLO), is related to that used to identify the optimum layout of bars in trusses, with discontinuities (e.g. slip-lines) in a translational failure mechanism corresponding to bars in an optimum truss. Hence, a recently developed adaptive nodal connection strategy developed for truss layout optimization problems can advantageously be applied here. The procedure is used to identify critical translational failure mechanisms for selected metal forming and soil mechanics problems. Close agreement with the exact analytical solutions is obtained.


1. Introduction

Pioneering theoretical developments in the field of plasticity in the middle of the last century led to the development of simple yet practical and powerful analytical tools which could be used to rapidly estimate the limit loads of simple bodies and structures. Chen (1975), for example, describes a wide range of methods of identifying upper- and lower-bound solutions for common geotechnical engineering problems, many of which can be partially automated. Of these, the method of characteristics (Sokolovskii 1965) has been successfully used to generate highly accurate solutions for several problem types. However, while this method is a powerful tool, it suffers from several limitations. For example, it provides incomplete lower-bound solutions and considerable insight from the operator is needed to identify the global form of a solution. Solutions thus tend to be tied to particular problem types. Engineers, however, increasingly demand more generally applicable methods, which normally require the use of generic numerical rather than analytical methods. Unfortunately, and despite their potential usefulness in practice, relatively few generally applicable numerical methods have to date been explored for this application. Of the methods which have been explored for plane plasticity problems, finite element limit analysis has proved popular with researchers over the past few decades (e.g. Lysmer 1970; Sloan 1988; Kobayashi 2005; Makrodimopoulos & Martin 2006). However, at the time of writing finite element limit analysis, which analyses the collapse state directly, has still to find widespread usage in engineering practice, with engineers instead being forced to rely either on analytical methods of limited applicability or on iterative elastic–plastic analysis methods, for example nonlinear finite element analysis. While the latter is certainly a powerful tool, it can suffer from numerical stability problems and significant operator expertise and time is required to prepare the requisite input data and to validate models.

Modern formulations of finite element limit analysis typically involve discretization of a body using both solid and interface elements, the latter placed between solid elements to permit jumps in the stress or strain rate fields (such formulations may therefore be considered as hybrid continuous–discontinuous analysis methods). For the solid elements, suitable finite element shape functions are used together with the desired yield criteria to ensure that the internal stresses are everywhere statically admissible (equilibrium formulation), or that the flow rule is everywhere satisfied (kinematic formulation). The failure surfaces used in limit analysis are generally nonlinear, although these may be linearized to permit the problem to be solved using linear programming (LP). Alternatively, problems involving certain nonlinear yield surfaces can be treated using efficient convex programming techniques (e.g. second order cone programming; Makrodimopoulos & Martin 2006). Unfortunately, the solutions obtained using finite element limit analysis are often highly sensitive to the geometry of the original finite element mesh, particularly in the region of stress or velocity singularities. Although meshes may be tailored to suit the problem in hand this is clearly unsatisfactory since advance knowledge of the mode of response is then required. Adaptive mesh refinement schemes can potentially overcome this problem (e.g. Lyamin et al. 2005), although it might be argued that the resulting analysis procedure is overly complex considering the simple rigid-plastic material idealization involved.

While finite element analysis is traditionally concerned with formulation and solution of a continuum mechanics problem, it is alternatively possible to consider a potentially simpler discontinuous problem. This involves directly identifying the discontinuities which form at failure (e.g. slip-lines, which transform a planar continuum into discontinua). However, previous attempts to develop discontinuous limit analysis formulations have been largely unsuccessful, principally because previous formulations have permitted only a severely limited range of possible failure mechanisms to be considered. For example, Alwis (2000), following on from a study by van Rij & Hodge (1978), prescribed that discontinuities could only coincide with the boundaries of rigid elements. This means, for example, that fan zones in failure mechanisms will not be identified unless the mesh takes that form at the outset (modern hybrid continuous–discontinuous finite element limit analysis formulations also prescribe that discontinuities can only coincide with element boundaries, though such formulations partly compensate for this by allowing deformations of the elements themselves).

In fact, a successful discontinuous limit analysis procedure must be able to identify the critical arrangement of discontinuities from a wide, preferably near-infinite, number of possibilities. With this in mind, an alternative approximation procedure to the traditional finite element method might involve discretization of a given planar body using a suitably large number of nodes laid out on a grid, with the failure mechanism comprising the most critical subset of potential discontinuities interconnecting these nodes. Posed in this form, the problem is now strikingly similar to the problem of identifying the optimum layout of discrete bars in trusses (e.g. ‘Michell’ trusses), a problem for which a well-established numerical solution procedure already exists (Dorn et al. 1964). This is perhaps not entirely unexpected as the analogy between certain plane plasticity and truss optimization problems was formally identified in the middle of the last century by workers such as Hemp (1958) and Prager (1959). Then, theoretical developments from the field of limit analysis were transferred to the field of optimal layout design of trusses. However, the computational layout optimization tools subsequently developed in the field of optimal layout design appear not to have been transferred back to the field of limit analysis. One objective of the present study was to seek to rectify this situation.

The intervening decades have seen considerable improvements in both available computing power and the efficiency of the mathematical programming solvers required to solve such problems. However, a further stimulus for the present work has been the development of an efficient adaptive nodal connection scheme, which has led to a significant increase in the size of tractable truss layout optimization problems (Gilbert & Tyas 2003). In the second part of this paper, the scheme is applied to plane limit analysis problems.

The paper begins by exploring the nature of the analogy between optimum trusses and optimum layouts of discontinuities, focusing particularly on approximate-discretized formulations which can be solved using LP.

2. Analogy between layout optimization of truss bars and discontinuities in Tresca (cohesive) material

(a) Background

The arrangement of slip-line discontinuities in plane strain metal plasticity problems is known to mirror the arrangement of bars in optimal Michell trusses; the geometry of both turn out to be Hencky–Prandtl nets, which are orthogonal curvilinear coordinate systems (Strang & Kohn 1983). Since a Michell truss contains an infinite number of infinitesimal truss bars they are sometimes referred to as ‘structural continua’, and hence may be treated using variational methods. However, Michell trusses are more commonly and conveniently treated as a series of discrete bars. In contrast, except when using simple hand analysis methods, plane plasticity problems are seldom treated as a series of discrete slip-lines (that, in turn, define a compatible rigid block mechanism). However, since the analogy holds at a fundamental level, such treatment must be possible, i.e. an analogy must also exist between approximate-discretized problem formulations for each problem type. It transpires that by posing both ‘equilibrium’ truss layout and ‘kinematic’ discontinuity layout LP formulations, the analogy between these becomes evident from inspection.

(b) Layout optimization of trusses: LP formulation

Consider a planar design domain which is discretized using n nodes and m potential connections (truss bars). The classical equilibrium plastic truss layout optimization formulation for a single load case is defined in equation (2.1) as follows:Embedded Imagesubject toEmbedded ImageEmbedded Image(2.1)where V is the total volume of the structure; Embedded Image; and Embedded Image are the tensile and compressive internal forces in bar Embedded Image; Embedded Image, where li and σi are, respectively, the length and the yield stress of bar i. B is a suitable 2n×2m equilibrium matrix and Embedded Image, where Embedded Image and Embedded Image are the x and y components of the external load applied to node Embedded Image. The presence of supports at nodes can be accounted for by omitting the relevant terms from f, together with the corresponding rows from B.

This problem is in a form which can be solved using LP, with the bar forces in q being the LP variables.

(c) Discontinuity layout optimization: LP formulation

The kinematic slip-line discontinuity layout optimization (DLO) formulation for the plane strain analysis of a quasi-statically loaded, perfectly plastic cohesive body discretized using m nodal connections (slip-line discontinuities), n nodes and a single load case is defined in equation (2.2) as follows:Embedded Imagesubject toEmbedded ImageEmbedded Image(2.2)where E is the total internal energy dissipated due to shearing along the discontinuities, where Embedded Image and Embedded Image are the relative shear displacement jumps between blocks along discontinuity Embedded Image; Embedded Image, where li and ci are, respectively, the length and the associated cohesive shear strength of discontinuity i. B is a suitable Embedded Image compatibility matrix and Embedded Image where Embedded Image and Embedded Image are the x and y components of the (virtual) displacement imposed at node Embedded Image. Rather than representing absolute Cartesian movements, these displacements correspond in a work sense to equivalent internal nodal forces, which can be altered to signify the presence of an external live load (refer to §4). In the absence of dead loads, energy balance requires that the internal energy E will also equal the work done by the external live load(s), thereby enabling the failure load (or failure load factor λ) to be determined. The discontinuity displacements in d are the LP variables. Note that for convenience the terms ‘energy dissipation’ and ‘displacement’ are used here as shorthand for ‘rate of energy dissipation’ and ‘displacement rate’, respectively.

By comparing equation (2.1) and (2.2) it is evident that the two problems are closely related and key features of the analogy are summarized in table 1. Figure 1a,b clarifies how equilibrium and compatibility conditions, respectively, are enforced at a specific node, also highlighting the relationship between the two problems. Furthermore, figure 1c provides an indication that there is scope to prescribe alternative compatibility criteria, for example if a dilatant material is present in the limit analysis problem.

View this table:
Table 1

Features of analogy between primal truss and slip-line discontinuity problem.

Figure 1

Equilibrium and compatibility conditions at a node (where Embedded Image and Embedded Image are, respectively, x-axis and y-axis direction cosines). (a) Equilibrium at (unloaded) node in a truss, (b) compatibility at node in slip mechanism (translation only) and (c) compatibility at node in slip mechanism (translation and dilation).

At this point it should be noted that a feature of formulations given in equations (2.1) and (2.2) is that truss bars and discontinuities, respectively, are free to crossover one another at locations other than nodes. While equilibrium and compatibility are not explicitly checked at crossover points, these are implicitly satisfied. This is because a bar (or discontinuity) which crosses over others at some point X may instead be considered as two separate bars (or discontinuities) which meet at X. Since these two bars (or discontinuities) share the same orientation and force (or slip displacement) value, they simply cancel each other out in the equilibrium (or compatibility) constraint equations. Thus, solutions containing crossover are perfectly valid (however, this is not to say that explicit addition of a node at crossover points would not lead to the generation of a more optimal solution). In the case of limit analysis problems, crossover can be considered to be particularly important as this effectively allows the domain to be divided into an extremely large number of potential sliding blocks, considerably more than when a conventional finite element mesh is used.

(d) Example problems

To further illustrate the analogy a number of closely related example problems will now be considered. Consider first a truss layout optimization problem. Figure 2a shows a 9×4 unit design domain and also prescribed support and load conditions. Nodes are positioned at unit intervals in the x and y directions and each node is connected to every other node, leading to what is often referred to as ‘a fully connected ground structure’ (though for the sake of clarity only discontinuities up to 3.5 units in length are shown in the figure). Pre-existing truss bars of infinite capacity are assumed to be present along the top edge of the domain, at y=4. Since the external load is vertically aligned, horizontal restraints have been omitted for the sake of simplicity. LP can now be used to identify the minimum volume truss layout(s). Figure 2b shows one of the several layouts which share the same optimal volume (V=64). The internal forces in the bars present in the optimal layout are also shown.

Figure 2

Analogy between truss and discontinuity layout optimization. (a) Truss layout problem: loading and supports, ground structure and locations of pre-existing bars (at y=4), (b) optimal layout of truss bars for problem (a), (c) discontinuity layout problem: imposed displacement conditions and optimal layout of slip-lines, i.e. as (b), (d) truss problem fully analogous to (c).

Consider next a planar limit analysis problem. Figure 2c shows the problem specification for a variant of the well-known Prandtl punch problem, considered by Hill (Chakrabarty 2006). This leads to the same critical layout of slip-line discontinuities as truss bars in figure 2b. It is assumed that there is zero cohesion along the boundary at y=4, coinciding with the locations of the pre-existing bars in the truss problem. The total internal energy dissipated (E=64) equals the energy dissipated by the external applied load, i.e. E=qBd, where q is the unknown limiting bearing pressure, B is the breadth of the footing and d is the imposed displacement. Thus, the limiting stress at failure q=E/Bd=64/(4×3)=5.33. This overestimates the exact solution of 2+π by less than 4% (with an increased nodal density this overestimate has been reduced to less than 0.04% by Gilbert et al. (2003), who considered the analogous truss problem). The slip displacements shown in the figure are equal in magnitude to the internal bar forces in the optimum truss solution shown in figure 2b. However, though the resulting layouts and bar forces/slip displacements are identical, the problems defined in figure 2a,c are not strictly speaking equivalent since the corresponding LP problem matrices are not identical. In fact the truss layout design problem which corresponds exactly to the limit analysis problem shown in figure 2c has rather more complex and unusual loading and boundary conditions, as shown in figure 2d.

Finally, figure 3 shows the hodograph (velocity diagram) for the slip-line discontinuity problem, which corresponds to the Maxwell force diagram for the truss.

Figure 3

Hodograph and Maxwell force diagram for the problem shown in figure 2.

(e) Commentary

The analogy between truss layout optimization and plane strain limit analysis problems has now been demonstrated; the DLO procedure generates a kinematically admissible discontinuous velocity field (or compatible rigid block mechanism) which represents a valid upper-bound solution according to the theorems of plasticity. However, it should be borne in mind that the direct analogy with truss optimization applies only when translational failure mechanisms and a weightless cohesive (Tresca) material are involved. While attention will be limited to translational failure mechanisms in this paper, there is a clear need to treat other common constitutive models (e.g. the Mohr–Coulomb failure criteria) and also to be able to include self weight.

Furthermore, it is not especially convenient to specify external loads and other boundary conditions in terms of nodal displacements. A more convenient and potentially powerful revised limit analysis formulation will therefore be outlined in §3.

3. Discontinuity layout optimization: revised kinematic formulation

(a) Extension to cohesive-frictional materials

The contribution of each potential discontinuity i to the global compatibility constraint equation in equation (2.2) can be written as follows:Embedded Image(3.1)where Bi, di and ui are, respectively, the local compatibility matrix, discontinuity and nodal displacement vectors. Alternatively, this may be expressed in expanded form for a cohesive (Tresca) material asEmbedded Image(3.2)where αi and βi are, respectively, the x-axis and y-axis direction cosines for discontinuity i, connecting nodes A and B; and the actual shear displacement at the discontinuity Embedded Image.

Now consider a cohesive-frictional material governed by the Mohr–Coulomb failure criteria. In this case, the associated flow rule can be written for discontinuity i asEmbedded Image(3.3)where ϕi is the angle of friction (dilation) and ni is the normal displacement accompanying the sliding, i.e. as shown in figure 1c. (It should be noted that while for real materials (e.g. sand) the frictional resistance cannot be entirely attributed to dilatancy, the so-called associative friction model adopted here leads to a good estimate of the limit load for problems which have few kinematic constraints.)

Since equation (3.3) stipulates that the normal displacement ni is simply a linear function of the shear displacement, this means that a combined compatibility and flow matrix can be derived and used when a cohesive-frictional soil is involved. Thus, equation (3.1) can be replaced with a new combined compatibility and flow constraint which, when expressed in the same form as equation (3.2), can be written asEmbedded Image(3.4)

However, the linear dependence of the normal displacement on the shear displacement makes this a rather special case. If more complex yield criteria are involved (e.g. if a no-tension condition were to be added) then it becomes necessary to adopt an alternative approach. With this in mind, a more generally applicable limit analysis formulation will now be presented, with new LP variables and with compatibility and flow constraints decoupled.

Thus, now using LP variables si and ni in di, to represent, respectively, the shear and normal displacement at discontinuity i, the local compatibility constraint becomesEmbedded Image(3.5)

The flow rule is then enforced by introducing the following constraint:Embedded Image(3.6)where Ni is a local plastic flow matrix; pi is a vector containing plastic multipliers Embedded Image, Embedded Image, where Embedded Image, Embedded Image. Each plastic multiplier corresponds to a linear yield surface constraint in the dual equilibrium problem formulation. Since a given discontinuity connecting two nodes obviously has known, fixed orientation, only two plastic multiplier variables are required to fully describe the Mohr–Coulomb cone.

(b) Specification of dead loads

In the LP problem formulation presented in equation (2.2) the objective was to minimize the internal energy E. This implicitly also leads to the minimization of the work done by the external live loads. However, in many problems dead loads (including body forces) are also present, requiring that additional terms are added to the work balance equation. Thus, an expanded work balance equation may be written as follows:Embedded Image(3.7)where Embedded Image and Embedded Image; and Embedded Image and Embedded Image, represent, respectively, the shear and normal dead and live loads applied locally at discontinuity Embedded Image.

It is convenient to distinguish between loads applied at external boundaries and those applied within a body itself. This is because whereas the displacements in d are absolute at external boundaries, the displacements are relative elsewhere. Thus, for a loaded boundary discontinuity i the local dead load vector fDi comprises simply the resolved shear and normal components of the applied load. For dead loads applied within a body, the contents of fDi can be obtained by summing up the total overlying dead load, excluding boundary dead loads. Consider for example a problem where the dead load results solely from self-weight effects. Here the contribution made by discontinuity i to the Embedded Image term in equation (3.7) can be written as follows:Embedded Image(3.8)where Wi is the total weight of the strip of material lying vertically above discontinuity i.

(c) Specification of live loads

In the LP problem formulation presented in equation (2.2) live loads had to be specified by directly modifying the nodal displacement values in u. These nodal displacements correspond in a work sense to equivalent internal nodal forces which need to be altered to account for the presence of the prescribed live load. This indirect approach is not especially convenient. However, live loads can instead be applied directly to individual discontinuities. Hence, now taking u=0, all live loads can be specified by imposing the following unit displacement constraint:Embedded Image(3.9)The contents of fLi, the local live load vector for discontinuity i, can be derived for boundary and within-body loads in the same way as described in §3b for fDi. If a rigid loading plate is present then additional constraints can be included to equalize the imposed displacements (e.g. if the applied load is applied normal to discontinuities lying on a straight line joining nodes A, B, C and D then the following constraints can be added: Embedded Image).

(d) Boundary conditions

It is convenient to define boundaries by simply adjusting the properties of boundary discontinuities. For the sake of simplicity it is assumed here that potential discontinuities spanning multiple nodes on boundaries are not present. Some common scenarios are now considered.

(i) Free boundary (including boundary subject to an applied load)

This may be achieved by setting each potential discontinuity i lying on a free boundary to have zero cohesion ci and also not applying flow rule constraint equation (3.6), thereby leaving ni and si unconstrained.

(ii) Fixed boundary

This is the default boundary condition, which requires no special treatment (i.e. the associated cohesion ci and angle of friction ϕi values for potential discontinuity i lying on such a boundary should be the same as those in non-boundary discontinuities, and the flow rule constraint equation (3.6) should be imposed as usual).

(iii) Line of symmetry

This may be achieved by setting each potential discontinuity i lying on the line of symmetry to have zero cohesion ci and angle of friction ϕi. The latter condition effectively constrains the normal displacement ni to be zero according to equation (3.6).

(e) Revised LP formulation

A revised primal kinematic problem formulation for the plane strain analysis of a quasi-statically loaded, perfectly plastic cohesive-frictional body discretized using m nodal connections (slip-line discontinuities), n nodes and a single load case can therefore be stated as follows:Embedded Imagesubject toEmbedded Image(3.10)where Embedded Image and Embedded Image are vectors containing, respectively, specified dead and live loads; d contains displacements along the discontinuities, Embedded Image, where Embedded Image and Embedded Image are the relative shear and normal displacements between blocks at discontinuity i; Embedded Image, where Embedded Image and Embedded Image are, respectively, the length and the cohesive shear strength of discontinuity i. Embedded Image is a suitable Embedded Image compatibility matrix; Embedded Image is a suitable Embedded Image flow matrix; and p is a 2m vector of plastic multipliers. The discontinuity displacements in Embedded Image and the plastic multipliers in p are the LP variables.

(f) On the efficiency of the basic method

To obtain accurate solutions each node in the problem should preferably be connected to all other nodes by a discontinuity (though overlapping discontinuities can be omitted, discontinuities which crossover one another should be included). However, when full connectivity is specified the resulting LP problem quickly becomes very large. Thus, the original and revised LP formulations described previously can be used on current generation personal computers to rapidly solve problems containing perhaps several hundred nodes. The accuracy so obtained will often be sufficient for routine engineering purposes. However, more accurate solutions can be obtained using an increased number of nodes. In fact problems with the order of several tens of thousand nodes can be treated using an adaptive connection adding scheme similar to that already used in truss layout optimization (Gilbert & Tyas 2003). Central to the scheme is consideration of the dual, equilibrium problem formulation. Since the equilibrium formulation also provides valuable additional insights into the DLO procedure, this is described in §4.

4. Discontinuity layout optimization: equilibrium formulation

(a) LP formulation

Duality principles can be used to derive the dual of equation (3.10), which is an equilibrium formulation. Thus, for a planar body discretized using m nodal connections (slip-line discontinuities) and n nodes this may be stated as follows:Embedded Imagesubject toEmbedded ImageEmbedded Image(4.1)where Embedded Image and Embedded Image and Embedded Image can be interpreted as x and y direction equivalent nodal forces acting at node Embedded Image, corresponding in a work sense to Embedded Image and Embedded Image, respectively; and q is a vector of shear and normal forces acting on discontinuities, i.e. Embedded Image, where Embedded Image and Embedded Image represent, respectively, the shear and the normal force acting on discontinuity Embedded Image. Therefore, the LP variables are Embedded Image, Embedded Image, Embedded Image and Embedded Image and the live load factor Embedded Image. The objective is thus to maximize Embedded Image while ensuring that the yield condition is not violated along any potential discontinuity.

The required equilibrium constraint can alternatively be written for a potential discontinuity i as follows:Embedded Image(4.2)or in expanded form asEmbedded Image(4.3)

Equation (4.3) can alternatively be derived from equilibrium considerations (e.g. figure 4 shows the forces involved, for clarity excluding external dead and live loads). Figure 5 illustrates the relationship between the discontinuity forces and the equivalent nodal forces for a simple problem.

Figure 4

Slip-line discontinuity forces in equilibrium with equivalent nodal forces (no external load case).

Figure 5

Relationship between discontinuity forces and equivalent nodal forces for a simple problem (2×1 unit block composed of weightless material possessing unit cohesive strength). (a) Problem definition: loaded block fixed to a base, (b) discretization using 6 nodes and 11 discontinuities, (c) computed equivalent nodal forces t at failure, (d) discontinuity forces (S, N) which can be derived from (c) (NAB uncorrected for λ); yielding lines highlighted.

The required yield constraint can also be written for a potential discontinuity i as follows:Embedded Image(4.4)or in expanded form for the Mohr–Coulomb yield condition asEmbedded Image(4.5)noting that here tensile forces are taken as positive.

Duality principles mean that values for Embedded Image, Embedded Image, Embedded Image and Embedded Image are available even if the primal problem is actually formulated and solved. However, it should be noted that Embedded Image and Embedded Image are free to take on arbitrary (though constrained) values except in yielding regions.

(b) Adaptive nodal connection procedure

To obtain more accurate results and to enable problems containing perhaps tens of thousands of nodes to be solved, a modified procedure employing adaptive addition of nodal connections, in this case discontinuities, is required. This is because the number of potential connections, and hence the number of LP variables, rises dramatically with the number of nodes, n. In fact the total number of possible connections, including overlapping connections, can be shown to be Embedded Image, and hence for a problem involving 10 000 nodes, there will be approximately 50 million potential connections (discontinuities), and a correspondingly large number of LP variables. Although overlapping connections can be omitted, the problem size will typically remain of the same order (the precise number of connections will then also depend on the relative spatial locations of nodes). Such a large problem cannot be solved directly using current generation personal computers and LP solvers.

Hence, it is not generally feasible to make all Embedded Image connections at the outset; instead it is more efficient to begin with minimal initial connectivity and to then add further connections as required as part of an iterative scheme. Thus, alternatively suppose that initially nodes are only partially connected (e.g. nodes are only connected to adjacent nodes so that m connections are represented in the LP problem, where Embedded Image), and a solution is then obtained. As the range of failure mechanism geometries that can be identified will be severely limited, clearly this solution is likely to significantly overestimate the true collapse load factor. What is required at this stage is a means of identifying which of the Embedded Image potential nodal connections can be added with a view to improving the solution.

Now, it transpires that the required yield constraint equation (4.5) can also be easily checked for a potential discontinuity i which is not presently represented in the current LP problem. This is achieved by firstly rearranging equation (4.3) so that the shear and normal force acting on the potential discontinuity between nodes A and B can be obtained from the solution of the last LP problem, using values of the internal nodal forces Embedded Image, Embedded Image, Embedded Image, Embedded Image and load factor Embedded Image.Embedded Image(4.6)where Embedded Image and Embedded Image are identical to Embedded Image and Embedded Image, respectively, except that they are not LP variables. Thus, while values for Embedded Image and Embedded Image are directly available for the m discontinuities currently represented in the LP problem, values for Embedded Image and Embedded Image can be computed for all Embedded Image other potential discontinuities. Furthermore, using the newly computed values of Embedded Image and Embedded Image, the yield constraint equation (4.5) can be checked for violation for the potential discontinuity. This checking process can be repeated at each iteration for all potential discontinuities, with any violating discontinuities becoming candidates for admission to an expanded LP problem at the next iteration. Using mathematical programming terminology the process of expanding the problem matrix is referred to as column generation when variables are added to the primal kinematic problem, or cut generation when constraints are added to the dual equilibrium problem.

However, as described in Gilbert & Tyas (2003) in the context of truss layout optimization, for maximum efficiency the number of discontinuities added at a given iteration should be limited to prevent the problem size from increasing too rapidly. This is achieved by identifying and then adding only the potential discontinuities where the dual inequality constraint is most violated. The iterative adaptive nodal connection procedure proceeds until no potential discontinuity violates equation (4.5). Once this condition is met the last solution obtained must also represent a globally optimal solution for the fully connected problem.

The adaptive nodal connection procedure therefore comprises the following steps: (i) discretize the problem using n nodes; (ii) minimally connect these nodes using a relatively small number of discontinuities; (iii) set up the corresponding LP problem; (iv) solve the LP problem; (v) for all discontinuities not represented in the current LP problem, check for violation of the yield condition using equation (4.6) and then equation (4.5); (vi) if violation is detected, add representation of the discontinuities where violation is greatest to the LP problem and repeat from step (iv); (vii) end.

5. Case study problems

To verify the accuracy of the solutions obtainable using the procedure, it will now be applied to a range of plane strain metal forming and geotechnical engineering problems, for some of which known analytical solutions already exist. All the problems chosen have pronounced singularities which can create difficulties when using other numerical solution methods.

To obtain the solutions a 2.4 GHz AMD Opteron PC equipped with 4 GB of memory and running 64 bit Scientific Linux was employed. Mosek v. 4.0, the commercially available interior point LP optimizer which uses the homogeneous and self-dual algorithm was used (Mosek 2006). The problem was initially passed to the optimizer in memory, using the supplied subroutine library; subsequently, it was only necessary to pass changes to the current problem to the optimizer rather than the entire revised problem. The presolve feature of the optimizer was enabled and default tolerances were used. Although the adaptive nodal connection procedure is amenable to parallelization, and a parallel version of the Mosek optimizer is available, a single processor was used for all computations.

Unless stated otherwise, the adaptive nodal connection scheme was used to solve all problems. For problems involving purely cohesive media (i.e. using the Tresca failure criteria), pairs of nodes not more than a distance of Embedded Image apart were initially connected (assuming unit spacing of nodes in the x and y directions). However, for problems involving frictional media (i.e. using the Mohr–Coulomb failure criteria) it was found that solutions could be obtained more quickly when the initial maximum connection distance was increased to Embedded Image. In fact when the angle of friction reached or exceeded 45° an initial solution could not otherwise be obtained for the bearing capacity problem considered, due to volumetric locking. It was also specified that not more than 5% of the number of connections present in the initial, minimal connectivity, problem could be added at each iteration. Even though changes to the LP problem at each iteration might be relatively modest, with the interior point optimizer used it was not possible to use the solution from a previous iteration as a starting point for the next optimization (i.e. it was not possible to perform a ‘warm start’).

Quoted CPU times include the time required to make all connections between the nodes and to set up and solve the LP problem(s). When the adaptive nodal connection scheme was used, the CPU times quoted also include the total time required to identify candidate connections for admission at the next iteration.

Finally, it should be noted that uniformly spaced nodes have been used for all problems studied. This meant that overlapping connections could easily be filtered out by checking that the greatest common divisor of Embedded Image and Embedded Image did not exceed unity (where Embedded Image, Embedded Image, and where (Embedded Image, Embedded Image) and (Embedded Image, Embedded Image) are the coordinates of nodes A and B which are potentially to be connected, also assuming unit spacing of nodes in the x and y directions).

(a) Metal block compressed between perfectly rough rigid platens

Consider a rectangular metal block compressed vertically between two parallel platens. Suppose that the platens are sufficiently rough to ensure that any slip at the platen–block interface occurs only when the shear stress reaches the prescribed cohesion. Chakrabarty (1991, 2006) has previously derived analytical solutions for this problem. The problem is of particular interest in the context of the present paper since the boundary conditions are such that the optimal layout of discontinuities obtained for half of the width of the block is identical to the optimal layout of truss bars in a well-known cantilever truss problem (e.g. refer to Lewinski et al. (1994)). The analogy was observed for this type of problem by Johnson (1961). While exact analytical solutions have been derived independently in each field, it is now evident that the mean pressure at failure q which can be applied to the compressed block can be related to the non-dimensional weight Embedded Image of the equivalent cantilever truss as follows:Embedded Image(5.1)where Embedded Image is the non-dimensional span of the cantilever (equal to half the width : height aspect ratio of the equivalent compressed block) and assuming unit limiting bar stress and block cohesion.

Taking advantage of symmetry, only a quarter of the block needs to be modelled. Solutions for block width : height ratios of 2.0, 3.64 and 6.72 are given in table 2. The locations of nodes and optimal layout of discontinuities for the 6.72 aspect ratio case are also given in figure 6. In the figure, the optimal layout of discontinuities has also been reflected about the horizontal line of symmetry to give a visual representation of the half block solution. Boundary conditions are also indicated on the reflected part of the figure. Note that for the sake of clarity only discontinuities with slip displacements Embedded Image have been plotted. Also note that the nodes were extended one division beyond the edge of the platen to model the influence of the remaining rigid block material.

View this table:
Table 2

Results for compressed metal block problem.

Figure 6

Layout of discontinuities at failure: compressed metal block between rough platens (169×25 nodal divisions; boundaries: F, Free; S, Symmetry plane; L, free with Live load).

It is clear from the results for this problem that close approximations of exact analytical solutions (less than 0.15% error) can be obtained reasonably rapidly (less than 700 s) using the DLO procedure.

(b) Embedded anchor in clay soil problem

Consider the problem of a horizontal anchor embedded in a homogeneous weightless clay. The goal is typically to determine the maximum pullout force for various embedment ratios (anchor plate depth : breadth). There appear currently to be no exact analytical solutions available for this problem, although numerical upper- and lower-bound solutions are available (Merifield et al. 2001). It is convenient to express results in terms of a dimensionless breakout factor Embedded Image where the ultimate pullout force Embedded Image in a weightless cohesive (i.e. Tresca) soil is given by Embedded Image, where B is the plate width and Embedded Image the undrained shear strength of the soil.

For each analysis, a line of symmetry was defined along one vertical edge, and a rectangular void of height one nodal division was positioned directly beneath the location of the loaded rigid anchor plate. Free boundaries were prescribed along the base and one edge of the void. It was also prescribed that discontinuities could not cross the void. The breakout factor was computed for an anchor embedment ratio of 8 for a range of nodal densities. The results are presented in table 3. Figure 7 shows the upper-bound mechanism for the 128×288 nodal division case (129×289=37 281 nodes, giving 415 807 859 potential discontinuities after removing overlaps). The optimal layout of discontinuities has been reflected about the vertical line of symmetry to give a visual representation of the full solution. Boundary conditions are also indicated on the reflected part of the figure. For the sake of clarity only discontinuities with slip displacements Embedded Image have been plotted in the figure.

View this table:
Table 3

Results for embedded anchor problem.

Figure 7

Layout of discontinuities at failure: uplift of embedded anchor in clay problem (128×288 nodal divisions; boundaries: F, Free; R, Rigid; S, Symmetry plane; L, free with Live load).

The results for all nodal densities appear to demonstrate an improvement on the results of Merifield et al. (2001). While Merifield et al. (2001) do not provide exact values, an examination of their approximate curve-fitted equations and graphs would indicate that their upper-bound solution for Embedded Image was above 7.6 and lower bound was below 7.1 for an embedment ratio of 8. It is also clear from the visual output (figure 7) that the new procedure is capable of providing new insights into the mode of response.

(c) Soil bearing capacity problems

(i) Frictional weightless soil with surface surcharge

Now consider the bearing capacity of a strip footing acting on a frictional weightless soil subject to a uniform surcharge pressure q. The bearing pressure at failure may be expressed as Embedded Image, where Embedded Image is referred to as a bearing capacity factor. Exact analytical solutions for Embedded Image were established by Reissner (1924). Here the new procedure was used to determine Embedded Image using various nodal densities. In each case, the nodes were laid out on a rectangular domain of aspect ratio 3 (length:height). A symmetry plane was defined along one vertical edge and a rigid footing placed adjacent to this along the top surface of the soil. The footing width was reduced for increasing angles of friction to account for the increasing extent of the failure mechanism. Problems with low, medium and high nodal densities were set up which comprised, respectively, 48×16, 96×32 and 192×64 divisions between nodes (i.e. 49×17, 97×33 and 193×65 nodes). After filtering out overlapping connections between nodes a total of 210 768, 3 116 344 and 47 839 456 potential discontinuities were present in these problems. The adaptive nodal connection procedure was used for all runs, although for the 48×16 nodal divisions problem the available memory was sufficient to permit the full initial connectivity problem also to be solved. Results are presented in table 4.

View this table:
Table 4

Results for bearing capacity on a surcharged frictional weightless soil problem.

From table 4 it is evident that reasonable approximations of the bearing pressure were obtained even when a low nodal density was used, and that closer estimates could be obtained when the nodal density was increased (albeit at considerable extra computational cost). It is also clear that, as well as permitting larger problems to be solved in the first place (due to considerably reduced memory requirements), the adaptive nodal connection procedure also permits solutions to small problems to be obtained more quickly. Finally, as expected it is found that the accuracy of the solution degrades with increasing angle of friction. This is at least partly a consequence of the larger extent of the failure mechanism in a soil with a high angle of friction, which means for example that the important zone below the footing becomes increasingly sparsely populated with nodes, when the nodal density is uniform throughout the problem domain.

Figure 8 shows a typical solution to a problem when 192×64 divisions between nodes were specified, in this case when the angle of friction was 25°. For the sake of clarity only discontinuities with slip displacements Embedded Image have been plotted in the figure. As expected, discontinuities within the body of different sign meet at an angle of approximately 90±25°, rather than approximately orthogonally, as would be the case when a cohesive (Tresca) material is involved.

Figure 8

Layout of discontinuities at failure: Embedded Image footing problem (ϕ=25°; 192×64 nodal divisions).

(ii) Cohesive-frictional soil with self weight

Martin (2003) has recently made available a software program (ABC) which uses the method of stress characteristics to obtain highly accurate (partial) lower-bound solutions specifically for bearing capacity problems. This means that it is now possible to obtain accurate benchmark solutions for practically relevant problems which involve cohesive-frictional soil possessing self weight. For example, consider a 1 m wide rigid strip footing on a soil with cohesion c=5 kPa, angle of friction ϕ=30° and soil unit weight γ=20 kN m−3. Assuming a frictionless interface between the footing and soil, the computed ABC limiting footing pressure for this problem is 268.312 kPa with zero surcharge.

A DLO solution obtained using 288×96 divisions between nodes is given in figure 9. Note that for the sake of clarity only discontinuities with slip displacements Embedded Image have been plotted in the figure. The computed limiting footing pressure for this example was 269.186 kPa, obtained in 25 012 s. While this is 0.3% greater than the ABC solution, which was obtained almost instantaneously on a desktop PC, the DLO procedure can of course be applied to a very much wider range of problems. Furthermore, on the assumption that the ABC solution is a true lower bound, this indicates that DLO is capable of obtaining answers very close to the true solution.

Figure 9

Layout of discontinuities at failure: footing in cohesive-frictional soil with self weight (288×96 nodal divisions).

6. Discussion

The real novelty of the DLO procedure lies in the expression of the limit analysis problem formulation entirely in terms of lines of discontinuity, rather than in terms of elements (as is usually the case). Using DLO a large number of potential discontinuities are set up and these take up many different orientations and are free to cross each other arbitrarily. In contrast, with element-based formulations discontinuities are typically restricted to lie only at the edges of elements positioned in a fixed mesh, and, since overlapping elements are not normally allowable, discontinuities are not permitted to cross. This means that when a critical mechanism is being sought the search space is very much more restricted.

Additionally, since discrete discontinuities are explicitly identified by the DLO procedure, the output is very different from that obtainable using continuum finite elements. In fact, one of the most attractive features of the new procedure is the way in which the form of hitherto poorly understood failure mechanisms can be readily visualized. For problems of potentially broad interest (e.g. the embedded horizontal anchor in clay problem), the newly identified layouts can potentially be used to help others to derive new analytical solutions. Alternatively, for more specific practical problems, the very clear visual output for example helps engineers to easily distinguish between rigid and yielding zones. Furthermore, as is evident in figure 8 for example, the procedure can automatically identify singularities in the displacement (velocity) field, something which in contrast can be difficult when using conventional finite elements.

The proposed procedure has been found to be very robust in use and can be guaranteed to determine the most critical layout of discontinuities for a specified arrangement of nodes. While the proximity of the numerical solution to the true solution in general remains unknown, performing a nodal density study provides a good indication of the improvement in accuracy which can be expected as the nodal spacing tends towards zero. Particularly when cohesive (Tresca) materials are involved the accuracy of the solutions will typically be found to be extremely good, even when relatively low nodal densities are used.

The use of uniformly spaced nodes in this paper contrasts with usual practice followed by most researchers involved in the development of finite element limit analysis methods, who have tended to use highly non-uniform meshes, often painstakingly tailored to suit the particular problem in hand (thereby potentially giving a false impression of the probable accuracy of the method when used in routine practice). However, it should be pointed out that it is not necessary to arrange nodes on uniformly spaced grids; the density of nodes can be increased in zones of interest, either at the outset by the user or, perhaps in the future, automatically during the course of the solution procedure as part of an adaptive nodal refinement scheme. Additional boundary nodes can also be added to more accurately model the shape of a complex body.

Furthermore, problems containing distinct zones which might have differing values of unit weight, cohesion and angle of friction can be treated. In the case of frictional media, all potential inter-zone discontinuities should use the highest angle of friction encountered to maintain the upper-bound status of the solution. Arbitrarily varying cohesion and self weight can be treated by, respectively, computing the total cohesive resistance and total weight of soil overlying a given potential discontinuity.

Finally, while indicative CPU times have been given, since the technology is relatively undeveloped (e.g. cf. the finite element method), it is probable that order of magnitude level reductions in CPU time can be achieved by diligent refinement of the formulation and adaptive nodal connection scheme heuristics.

7. Conclusions

  1. The analogy between approximate-discretized formulations for truss layout optimization and slip-line discontinuity layout optimization has been established.

  2. This has led to the development of a potentially important new analysis procedure, discontinuity layout optimization, which has been described here in the context of limit analysis of plane plasticity problems. The procedure involves identification of the subset of inter-node discontinuities which are present in the critical failure mechanism. The accuracy of the solution can be improved by increasing the density of nodes covering the body under consideration, which in turn increases the number of discontinuities available for possible inclusion in the critical mechanism.

  3. In common with conventional finite element limit analysis, the DLO procedure is flexible enough to handle the kind of generic limit analysis problems which would challenge, for example, the method of characteristics. However, unlike finite element limit analysis but in common with the method of characteristics, the procedure can handle stress/velocity singularities with ease and can also provide output that is amenable to immediate visual interpretation by the user.

  4. Highly accurate upper-bound translational mechanism solutions to problems involving cohesive (Tresca) material have been obtained using the procedure. Problems involving frictional media (and the Mohr–Coulomb failure criteria) were found to be more challenging, although promising results for problems with and without self weight were nevertheless obtained.

  5. The DLO procedure seems to retain much of the inherent simplicity of the traditional ‘upper-bound’ analytical solution method used in classical plasticity. This contrasts with conventional finite element limit analysis, which is necessarily more complex.


The financial support of EPSRC is gratefully acknowledged by the second author (Advanced Research Fellowship grant GR/S53329/01). The authors would also like to thank Dr Andy Tyas for general encouragement and entering into valuable discussions on this topic; Mr Wael Darwich for co-developing the cross-platform software framework employed when solving the case study problems.


    • Received October 10, 2006.
    • Accepted May 30, 2007.


View Abstract