## Abstract

A new three-dimensional limit analysis formulation that uses the recently developed discontinuity layout optimization (DLO) procedure is described. With DLO, limit analysis problems are formulated purely in terms of discontinuities, which take the form of polygons when three-dimensional problems are involved. Efficient second-order cone programming techniques can be used to obtain solutions for problems involving Tresca and Mohr–Coulomb yield criteria. This allows traditional ‘upper bound’ translational collapse mechanisms to be identified automatically. A number of simple benchmark problems are considered, demonstrating that good results can be obtained even when coarse numerical discretizations are employed.

## 1. Introduction

The formal theorems of plastic limit analysis provide the theoretical framework necessary to allow direct evaluation of the load required to cause collapse of a body or structure, without the need for intermediate calculation steps. While early work in the field focused on the development of hand-type limit analysis methods, which are still widely used in engineering practice, more recent research has focused on the development of computational methods, such as finite-element limit analysis [1–3]. Computational limit analysis methods have successfully been used to obtain bounds on the collapse load for plane strain problems, providing a rapid and robust means of evaluating safety. However, application of limit analysis techniques to three-dimensional problems has been more limited. This is unfortunate since engineers are increasingly seeking to model the real three-dimensional geometry of a given problem, rather than a two-dimensional idealization.

Research on three-dimensional problems has focused primarily on exploiting the potential of finite-element limit analysis formulations, making use of optimization techniques to obtain upper or lower bound solutions. The effectiveness of a given formulation is normally closely linked to that of the optimization technique used to obtain solutions, and the yield criterion is key in determining which optimization technique is appropriate for a particular formulation.

Park & Kobayashi [4] were among the first to develop a limit analysis formulation for three-dimensional problems, referred to as the ‘rigid-plastic finite-element method’. Their approach involved relaxing the flow rule to obtain a direct estimation of the collapse load. Rigid-plastic finite elements have successfully been used to investigate a wide range of metal-forming applications. However, the method appears only to have been applied to von Mises or Drucker–Prager materials, and the solutions obtained do not have clear upper or lower bound status within the context of the fundamental theorems of plastic limit analysis.

Lyamin & Sloan [5,6] proposed a three-dimensional finite-element limit analysis formulation that used nonlinear programming to obtain solutions. Their formulation could be applied to problems involving any yield function that is everywhere differentiable. However, this precludes direct use of the Mohr–Coulomb yield criterion owing to the singularity that exists at the apex of the Mohr–Coulomb cone. Consequently, Lyamin & Sloan [5] proposed the use of a hyperbolic smoothing technique to create an approximated yield surface, differentiable everywhere. More recently Krabbenhøft *et al.* [7] and Martin & Makrodimopoulos [8] have used semidefinite programming (SDP) to create formulations capable of handling the Mohr–Coulomb yield function directly. However, at present SDP algorithms are relatively immature, limiting the size of problems that can be solved.

Chen *et al.* [9] used fully rigid finite elements to obtain upper bound solutions for three-dimensional slope stability problems. A distinguishing feature of their formulation is that, since the elements themselves are not free to deform, deformations can only take place along discontinuities lying at predefined element boundaries. An advantage of this approach is that the form of the critical failure mechanism is clear, aiding interpretation by engineers. However, a major disadvantage is that the topology of the finite-element mesh used will often have a major influence on the form of the mechanism identified, which may be very different in form to the true critical mechanism. A principal aim of the research described in the present paper was to develop a method that overcomes this obvious drawback, allowing a much wider range of failure mechanisms to be modelled, and without having to turn to SDP algorithms to obtain a solution. The development of such a method is described, and its efficacy demonstrated through application to various benchmark problems considered in the literature.

## 2. Solution of three-dimensional rigid finite-element problems

The three-dimensional rigid finite-element (or ‘rigid block’) method applied to Mohr–Coulomb materials by Chen *et al.* [9] used quadratic programming to obtain solutions. However, such a problem can alternatively be recast as a second-order cone programming (SOCP) problem. Optimization problems taking this form can be efficiently be handled by modern interior-point algorithms, similar to those used to solve linear programming problems [10].

To explain how this is achieved, it is convenient to use an orthogonal coordinate system local to each discontinuity, comprising axes **n**, **s** and **t**, where **n** is a unit vector normal to the discontinuity and **s** and **t** are unit vectors in the plane of the discontinuity. Considering translational mechanisms, the Mohr–Coulomb yield criterion can now be enforced for stress resultants acting on the plane of the discontinuity by equations (2.1) and (2.2)
2.1and
2.2where *c* and *ϕ* are the material cohesion and angle of shearing resistance and *a* is the face area of the discontinuity. Here *N*, *S* and *T* denote the normal force and shear traction components along the **n**, **s** and **t** axes, respectively, and *P* is the maximum shear traction on the discontinuity, as indicated in figure 1. Similarly, the associative flow rule for a Mohr–Coulomb material can be expressed by equations (2.3) and (2.4)
2.3and
2.4where *n*, *s* and *t* denote the component of the relative jump in displacement rate across the discontinuity in the **n**, **s** and **t** directions, respectively, and *p* is a plastic multiplier, as indicated in figure 1. (Note that henceforth, ‘energy dissipation’ and ‘displacement’ will be used as shorthand for ‘rate of energy dissipation’ and ‘rate of displacement’, respectively.) Using equation (2.4), it is now possible to formulate the problem as an SOCP problem.

In the formulation outlined earlier, individual three-dimensional elements of material are rigid, and not free to deform. It therefore seems natural to use a formulation where the problem is posed entirely in terms of discontinuities. The discontinuity layout optimization (DLO) procedure, outlined for two-dimensional plane plasticity problems by Smith & Gilbert [11], is posed entirely in terms of discontinuities, and its potential application to three-dimensional problems is therefore considered here. An advantage of DLO is that it can inherently model singularities, such as those that exist at the edge or corner of an indentor or foundation footing. Conversely, finite elements, both rigid and deformable, generally require refinement of the mesh in the regions of singularities [12], or alternatively the use of higher order shape functions [13] in order to obtain reasonable results. The aim of this paper is to develop a DLO-based formulation that can be applied to three-dimensional problems involving Tresca and Mohr–Coulomb materials.

## 3. Discontinuity layout optimization

### (a) Background

DLO is a direct limit analysis method, which means that it can directly provide an estimate of the maximum load sustainable by a solid body, without intermediate calculation steps. DLO can be used to generate upper bound plasticity solutions and associated collapse mechanisms.

The DLO procedure for plane strain problems is outlined in figure 2. Firstly, the initial problem is discretized using nodes distributed across the body under consideration. Potential discontinuity lines (e.g. ‘slip-lines’), along which jumps in displacements **d** can occur, are set up by linking each node to every other node, and optimization is used to identify the subset of discontinuities active in the critical failure mechanism. Provided a sufficiently large number of nodes are employed, this allows a very wide range of potential mechanisms to be considered.

For plane strain problems, the corresponding optimization problem is given in equations (3.1)–(3.5), posed entirely in terms of potential discontinuities [11], 3.1subject to 3.2 3.3 3.4 and 3.5

In equation (3.2), compatibility is enforced explicitly at each node via a suitable compatibility matrix **B** containing direction cosines. Equation (3.3) enforces the flow rule on each discontinuity by introducing a vector of plastic multipliers **p** and a suitable flow matrix **N**, while equation (3.4) ensures the work done by the live loads (**f**_{L}) equals unity (where here the term ‘live loads’ is taken to encompass all loads that are to be varied as part of the optimization process). Equation (3.1) is the objective function, where the goal is to obtain the minimum load factor *λ* by minimizing the work done by dead loads (**f**_{D}) and internal plastic energy dissipation **g**^{T}**p** (where **g** contains the corresponding dissipation coefficients).

Linear programming is then used to obtain the minimum value of *λ*, seeking the optimal values of the problem variables in **d** and **p**. The critical subset of potential discontinuities, defining the critical mechanism, is also identified, as indicated in figure 2*d*.

In this formulation, various intersections between discontinuities (or ‘crossover points’) occur throughout the problem domain. However, the displacement along each discontinuity remains unchanged either side of a ‘crossover point’, and hence compatibility is inherently maintained.

### (b) Three-dimensional formulation

In the following sections a three-dimensional formulation is developed, using a three-dimensional grid of nodes and polygonal discontinuities. While any simple polygonal shape, or combination of simple polygonal shapes, may be used for the discontinuities, triangular discontinuities provide the most flexibility and will be used in the numerical examples. However, for the sake of clarity when describing the method, discontinuities of both rectangular and triangular shape will be used. In general, when using triangular discontinuities there will be a total of *n*(*n*−1)(*n*−2)/6 potential discontinuities, where *n* is the total number of nodes used to discretize the problem (cf. a total of *n*(*n*−1)/2 potential line discontinuities in a plane strain problem). As with plane strain problems, this results in a rich field of potential translational mechanisms, and reduces the need for refinement in the region of singularities.

#### (i) Compatibility

In plane strain DLO, compatibility is enforced at nodes. In three-dimensional DLO, compatibility is conveniently enforced along shared edges. Figure 3 shows a set of triangular prisms sharing a common edge *OO*′ and with absolute displacements, **v**_{AB} to **v**_{EA}. The prisms are separated by rectangular discontinuities *O*′*OXX*′, where *X*=*A*, *B*…*E*. Each discontinuity *O*′*OXX*′ has a normal **n**_{X} and a relative displacement jump **Δ****v**_{X} that denotes the difference in absolute displacement of the prisms meeting at this discontinuity. It follows from this definition that equation (3.6), which involves summing all relative displacement jumps around edge *OO*′, must hold
3.6

A sign convention for determining the directions of **n**_{A} to **n**_{E} and **Δ****v**_{A} to **Δ****v**_{E} is presented in appendix A. Using this sign convention and the vertex ordering in figure 3*b*, the following must also be true along edge *OO*′:
3.7

Compatibility can be similarly enforced along the remaining edges, using the sign convention given in appendix B. Moreover, equation (3.7) can be reformulated for each edge. Before proceeding, it is first necessary to use coordinate systems local to each discontinuity *i*
3.8where **T**_{i} is a 3×3 transformation matrix converting local to global displacement jumps. **T**_{i} is chosen such that one axis aligns with **n**_{i}, the unit column vector in the normal direction using the ‘right-hand screw rule’; the two remaining orthogonal axes, **s**_{i} and **t**_{i}, are in the plane of discontinuity *i* such that **T**_{i}={**n**_{i} **s**_{i} **t**_{i}}. Also **d**_{i}={*n*_{i} *s*_{i} *t*_{i}}^{T}, where *n*_{i}, *s*_{i} and *t*_{i} are the local displacement jumps across discontinuity *i* in the **n**_{i}, **s**_{i} and **t**_{i} directions, respectively.

Compatibility can now be enforced for each edge *j* (*j*=1,2,…,*l*) as follows:
3.9where *S*_{j} is a subset of the total number of discontinuities *m* that contains the discontinuities meeting at edge *j*, and **B**_{ij} is a local compatibility matrix equal to *k*_{ij}**T**_{i}, where *k*_{ij} is defined in appendix B.

The DLO procedure results in intersections and overlaps between discontinuities that do not coincide with nodal connections; however, compatibility at these is enforced implicitly, as demonstrated in appendix C.

#### (ii) Flow rule

The local coordinate system described allows the associative flow rule for a Mohr–Coulomb material to be enforced using equations (2.3) and (2.4). Taking advantage of SOCP, the flow rule on a discontinuity *i* can be restated as follows:
3.10and
3.11where *p*_{i} is a plastic multiplier and where *ϕ*_{i} is the angle of friction on discontinuity *i*. Equation (3.10) is a linear constraint and equation (3.11) is a second-order cone. In matrix form, equation (3.10) can be stated as follows:
3.12where **d**_{n} is a subset of **d** containing only the displacement jumps normal to discontinuities *i* (*i*=1,…,*m*), **p** is a global vector containing plastic multipliers and *m* is now the total number of discontinuities in the problem. **N** is a global *m*×*m* matrix enforcing the flow rule and equal to . However, depending on the material properties, it may not be necessary to apply constraint equations (3.10) and/or (3.11), as indicated in the first half of table 1.

#### (iii) Dissipation function

The energy dissipated on a given discontinuity *i* is given by *g*_{i}*p*_{i}, where *g*_{i} is a dissipation coefficient equal to , the integral of the cohesive strength over the area *a* of discontinuity *i*. In the case of uniform cohesive strength across discontinuity *i*, the dissipation coefficient *g*_{i}=*a*_{i}*c*_{i}, where *a*_{i} and *c*_{i} are, respectively, the area and cohesion of the discontinuity.

On overlapping regions, the upper bound status of the solution is maintained, as demonstrated in appendix C.

### (c) Mathematical formulation

A three-dimensional kinematic formulation for a cohesive-frictional body discretized using *m* polygonal discontinuities and *l* edges can be summarized as follows:
3.13subject to
3.14
3.15
3.16
and
3.17where **f**_{D} and **f**_{L} are vectors containing, respectively, the specified dead and live loads. Here **d** contains displacement jumps across the discontinuities, **d**^{T}={*n*_{1},*s*_{1},*t*_{1},*n*_{2},*s*_{2},*t*_{2},…,*t*_{m}}, where *n*_{i} is the displacement jump normal to discontinuity *i*,*s*_{i} and *t*_{i} are the displacement jumps within the plane of discontinuity *i* and **g** is a vector of dissipation coefficients. Here **B** is a suitable 3*l*×3*m* compatibility matrix, **N** is a suitable *m*×*m* flow matrix and **d**_{n} (a subset of **d**) is a vector containing the normal displacement jumps, **d**^{T}_{n}={*n*_{1},*n*_{2},…,*n*_{m}}. Also **p** is a vector of plastic multipliers, **p**^{T}={*p*_{1},*p*_{2},…,*p*_{m}}, where *p*_{i} is the plastic multiplier for discontinuity *i* given by equation (3.17). The optimization variables are the displacement jumps in **d** (and **d**_{n}) and the plastic multipliers in **p**. The objective function and the first three constraints are linear. The final constraints on the plastic multipliers *p*_{i} are second-order cones, so that the formulation is amenable to solution using SOCP. The problem can also be posed in an equilibrium form, established using duality principles [15].

### (d) Boundary conditions and loads

Many common boundary conditions can readily be modelled simply by using a reduced number of variables and constraints, as indicated in the second half of table 1. This for example shows that a discontinuity on a rigid boundary is dealt with in exactly the same manner as a discontinuity in the interior of a domain.

Dead and live loads **f**^{T}_{D} and **f**^{T}_{L} in equations (3.13) and (3.16) are now defined such that and , where , , and , , are, respectively, the dead and live loads acting in the **n**_{i}, **s**_{i} and **t**_{i} directions on discontinuity *i*. Areas of flexible loading can be applied directly to the discontinuities with no special treatment. Rigid loads can be applied using a discontinuity covering the whole loaded area and reducing the degrees of freedom of underlying discontinuities appropriately.

At an external boundary, after taking account of any overlapping regions, displacement jumps must equal the absolute displacement of that boundary. Hence **f**_{Di} and **f**_{Li} are simply the local dead and live loads, respectively, on discontinuity *i*, resolved in the direction defined by the local coordinate system when applied at boundary discontinuity *i*. For discontinuities within a body, the contents of **f**_{Di} and **f**_{Li} can be obtained by summing up the total overlying dead or live loads, excluding boundary loads. For example, for the case where dead loads are due to self weight only, and assuming this is applied in the negative *z* direction, the contribution to the summation made by discontinuity *i* is as follows:
3.18where **W**_{i} is a 1×3 row vector containing the components in the **n**_{i}, **s**_{i} and **t**_{i} directions of the total weight of the column lying vertically above discontinuity *i*.

### (e) Summary of procedure

Steps in the DLO procedure for three-dimensional problems can be summarized as follows:

— discretize the problem using nodes;

— connect nodes to create edges;

— join edges to create polygonal discontinuities;

— solve the resulting SOCP problem.

## 4. Numerical examples

To evaluate the potential of the method, various three-dimensional examples are now considered. Unless stated otherwise, all computations were performed using a workstation equipped with 2.6 GHz AMD Opteron 6140 processors, 8 GB RAM and running 64 bit Scientific Linux. The Mosek interior-point solver with SOCP capability was used [16]. The default settings of the optimizer were used, including the pre-solve feature. The CPU times reported are for the optimizer only, including the time taken for the pre-solve routine to execute, but excluding the time taken to read in and set up a given problem.

Prior to solving, various measures were taken to condition and/or reduce the size of a given problem. Firstly, the coefficients in the objective function (3.13) and unit displacement constraint (3.16) were scaled as necessary to ensure the problem was well posed. Secondly, as overlapping edges do not provide any extra degrees of freedom, these were removed. Thirdly, discontinuities covering areas that could be reconstructed by combining several discontinuities covering smaller areas were removed. Fourthly, noting that the formulation naturally results in 3(*n*−1) linear dependencies in the constraint matrix, where *n* is the total number of nodes, such linear dependencies were identified and removed prior to passing the problem to the optimizer. Lastly, while the basic DLO procedure involves positioning nodes on a Cartesian grid, the use of other nodal arrangements is possible, and, where clearly indicated, regular grids with differing *x*, *y* and *z* spacings are used in this paper. However, it should be noted that an adaptive solution procedure capable of dramatically reducing problem size (e.g. of the sort described in [11] for plane strain problems) was not used in the numerical studies described herein. Thus the size of a problem increases rapidly as nodal resolution is increased, and consequently only relatively small problems are considered here.

### (a) Compression of a block

The unconfined compression of a square block with shear strength parameters *c*, *ϕ* between two perfectly rough rigid platens has previously been considered by Martin & Makrodimopoulos [8], who have obtained upper and lower bound solutions for this problem. This problem will therefore be used as a benchmark for the proposed three-dimensional DLO procedure. For the geometry shown in figure 4, the objective is to find the ratio of the average bearing pressure *q* to cohesive shear strength *c*.

Symmetry means that only 1/16 of the block needs to be modelled, as indicated in figure 4. Nodes were initially positioned on a Cartesian grid (i.e. equal nodal spacings in the *x*, *y* and *z* directions) and the solutions for various nodal spacings sought; these are presented in table 2. For *ϕ*=0^{°}, the best new solution presented is close (within 0.39%) to the benchmark upper bound solutions of Martin & Makrodimopoulos [8]. A representative collapse mechanism is shown in figure 5. In the case of *ϕ*=30^{°}, the best new solution compares less favourably with the benchmark (difference 10%).

Noting that the active discontinuities in figure 5 radiate from the centre, it is of interest to investigate the effect of using different nodal spacings in the *x* direction compared with those in the *y* and *z* directions. Results for the *ϕ*=0^{°} case for various nodal spacings are presented in table 3.

Firstly, it is evident that a solution matching the best reported value of 2.314 presented previously could be obtained in only 0.22 s, compared with 9700 s previously. This is because only a small subset of the discontinuities present previously are now included. Secondly, it is evident that it has been possible to improve on the benchmark upper bound solution (the best DLO solution of 2.304 is slightly less than the value of 2.305 obtained by Martin & Makrodimopoulos [8]).

### (b) Punch indentation

The bearing capacity of a perfectly rough square indenter resting on the surface of a purely cohesive Tresca material is now considered. The value of interest is once again the ratio *q*/*c*, otherwise known as the bearing capacity factor *N*_{c}. Salgado *et al.* [17] have established upper and lower bounds for a variety of indenter embedment depths and geometries using finite-element limit analysis. Gourvenec *et al.* [18] have modified the upper bound mechanism proposed by Shield & Drucker [19] for a smooth indenter to model the effect of a rough indenter. Michalowski [20] and Vicente da Silva & Antão [21] have also established upper bounds for a number of indenter geometries bearing onto the material surface. The best reported upper and lower bounds for a square indenter are included in table 4.

The problem geometry used is shown in figure 6. While both Michalowski [20] and Gourvenec *et al.* [18] have obtained mechanisms with only two planes of symmetry, the best upper and lower bound solutions, owing to Vicente da Silva & Antão [21] and Salgado *et al.* [17], respectively, have taken advantage of the four planes of symmetry inherent in the problem geometry. Therefore, only 1/8 of the problem has been modelled, as indicated in figure 6. Table 4 presents new solutions for three nodal spacings, with all nodes positioned on a Cartesian grid (i.e. equal nodal spacings in the *x*, *y* and *z* directions). These solutions compare well with the best reported upper bound, especially considering the comparatively low nodal resolutions employed. Also, following the initial study, access to a workstation with 16 GB RAM allowed consideration of a finer numerical discretization (with Δ=1/8), resulting in an improved *N*_{c}=6.102 (difference 0.84%).

A representative failure mechanism is shown in figure 7. It should be noted that the critical mechanisms for all three nodal grids extend up to the fixed boundaries. However, extending the problem domain relative to the foundation quickly leads to impractically large problem sizes, so this issue was not investigated further.

### (c) Anchor in a purely cohesive soil (*ϕ*=0)

Consider a perfectly rough anchor of width *B* embedded at a depth *H* in a purely cohesive Tresca soil, as shown in figure 8. Immediate breakaway (i.e. no suction or transmission of tensile stresses) is assumed between the anchor base and the soil. Taking advantage of symmetry, only 1/8 of the problem needs to be modelled, as shown in figure 8. Various *H*/*B* ratios have been considered by fixing *B*=2, *W*=10, Δ*x*=Δ*y*=1, Δ*z*=*H*/4 and varying *H*.

Merifield *et al.* [22] used a lower bound finite-element limit analysis procedure to establish bounds on the break-out factor at various embedment depths. The break-out factor *N*_{c0} for an anchor in a weightless cohesive soil is defined as average bearing pressure *q* divided by the cohesive strength *c*. For a soil with unit weight *γ*, Merifield *et al.* [22] defined the break-out factor *N*_{cγ} as
4.1Failure mechanisms can be classified as ‘shallow’, where the failure mechanism extends up from the anchor to the soil surface, and ‘deep’, where the mechanism involves only localized deformations around the anchor, and where the break-out factor *N*_{c*} is independent of the embedment depth *H*. For a given ratio *γ*/*c*, the deep mechanism becomes critical at depths greater than or equal to *H*_{cr}. If both deep and shallow mechanisms are considered, *N*_{cγ} must be less than or equal to *N*_{c*}. Merifield *et al.* [22] established a lower bound on *N*_{c*}≈11.9. For *γ*/*c*=0, Merifield *et al.* [22] found *H*_{cr}≈7.

Figure 9 shows break-out factors for shallow failures at various embedment depths. These show good agreement with those reported by Merifield *et al.* [22] when *γ*/*c*=0 and *H*<*H*_{cr}. Furthermore, the results for *γ*/*c*=1 and *γ*/*c*=2 are consistent with the relationship described in equation (4.1). The mechanisms developed did not extend to the lateral fixed boundary, suggesting that *N*_{c} and *N*_{cγ} were not influenced by the extent of the domain modelled.

### (d) Anchor in a purely frictional soil

Now consider a perfectly rough anchor of width *B* embedded at depth *H* in a purely frictional soil with an angle of shearing resistance *ϕ*=30^{°} and a unit self weight *γ*. Assuming an associative flow rule, only mechanisms that extend to the surface are possible. The break-out factor *N*_{γ} in this case can be expressed as
4.2where *q* is the average pressure on the anchor. Merifield *et al.* [23] used a lower bound finite-element analysis procedure to establish bounds on the break-out factor for various angles of friction and embedment depths. Also, Murray & Geddes [24] proposed an upper bound solution given in equation (4.3), obtained by assuming a simple rigid block mechanism,
4.3Using nodal grids of the form shown in figure 8 and fixing *B*=2, *W*=10, Δ*x*=Δ*y*=1, Δ*z*=*H*/4 (denoted GRID 1) and varying *H*, a very close fit to equation (4.3) has been found for 0.5≤*H*/*B*≤2.5 (figure 10). Above this range, the mechanism reaches the lateral fixed boundary and no solution was obtained. The critical mechanism is relatively simple, and very similar results can be obtained using an alternative grid, where *B*=2, *W*=16, Δ*x*=Δ*y*=1, Δ*z*=*H* (denoted GRID 2). GRID 2 also has been used to determine *N*_{γ} for 2.5<*H*/*B*≤5, which also compare well with equation (4.3).

## 5. Discussion

As can be seen in §4, the three-dimensional DLO procedure is capable of obtaining results comparable with those found in the literature. Also in one case an improved upper bound solution was found. A key benefit of the procedure is that it can automatically identify traditional ‘upper bound’ mechanisms, obviating the need for these to be found manually, as has hitherto been necessary. Such mechanisms are generally easy to visualize and can be checked by hand if required. Thus, a potentially useful application of the method is in helping designers identify critical mechanisms for subsequent use in hand calculations. While only translational mechanisms are currently identified, experience gained applying plane strain DLO to practical problems has indicated that solutions of acceptable accuracy for engineering purposes can often be obtained, even for problems where rotational modes are critical. Additionally, mechanisms that include continuous deformations can be captured by using high concentrations of active discontinuities.

However, refining the nodal grid even slightly leads to a large increase in the number of potential discontinuities, which means that problems can quickly become impractically large using currently available computational power. Thus enhancements to improve performance are required. One potential approach is to implement an adaptive solution procedure, similar to that proposed for plane strain problems by Smith & Gilbert [11]. With the latter procedure a subset of the total number of potential discontinuities is used in the initial optimization. The solution is then assessed, and discontinuities deemed most likely to improve it are added. A further optimization is then performed and the process repeated until a rigorous termination criterion is met, ensuring that the numerical solution is the same as would have been obtained had all potential discontinuities been included in the optimization at the outset. For all problems considered in this paper, the percentage of discontinuities with non-zero displacements was small, suggesting that a similar adaptive solution procedure might be effective, ensuring that the number of discontinuities modelled in the optimization problem is kept to a minimum. Developing such a procedure is therefore a topic for future research.

## 6. Conclusions

— A new three-dimensional plasticity formulation has been described. The formulation uses the DLO procedure to directly identify critical translational collapse mechanisms. Planar triangular discontinuities, which interconnect nodes distributed across the domain, are employed and SOCP is used to identify the optimal subset of discontinuities, which define the geometry of the collapse mechanism.

— A key benefit of the DLO procedure is that it can automatically identify traditional ‘upper bound’ collapse mechanisms, obviating the need for these to be found manually, as has hitherto been necessary.

— Despite the use of comparatively coarse nodal discretizations, reasonably good correlation with literature benchmark solutions has been observed. Furthermore, the best upper bound solution for the compression of a purely cohesive block between two perfectly rough platens problem reported in the literature has been improved upon.

— The use of an adaptive solution procedure, similar to that successfully used in plane strain DLO, is likely to allow finer nodal discretizations to be treated. This has therefore been identified as a topic for future research.

## Appendix A. Sign convention

In DLO, the optimization problem is formulated purely in terms of displacement jumps across a set of potential discontinuities. This requires a consistent sign convention to be adopted. Figure 11*a* shows two triangular prisms separated by a discontinuity *ABC* moving with absolute displacements **v**_{1} and **v**_{2}, respectively. The displacement jump **Δ****v**_{ABC} across discontinuity *ABC* must equal either **v**_{1}−**v**_{2} or **v**_{2}−**v**_{1}. Using the ‘right-hand screw rule’, the normal **n** to the discontinuity can be determined for a particular ordering of the vertices, as demonstrated in figure 11*b*,*c* for discontinuity *ABC*. **Δ****v**_{ABC} is simply defined as **v**_{+}−**v**_{−}, where **v**_{+} and **v**_{−} are the velocities of the blocks on the **n** positive and **n** negative sides of the discontinuity, respectively, for the chosen ordering of the vertices.

## Appendix B. Convention for selecting the vertex ordering parameter *k*_{ij}

A sign convention is necessary for equations (3.7) and (3.9) to be valid along edge *j* (where *j*=1,2,…,*l* and where *l* is the total number of edges used to discretize the problem). Taking account of the sign convention, these equations can be rewritten as follows:
B1and
B2where *S*_{j} is the subset of all discontinuities meeting at edge *j*. *k*_{ij}=±1 depends on the relative vertex orderings of discontinuity *i* and edge *j*.

Each discontinuity *i* (*i*=1,2,…,*m*) is defined by a subset *D*_{i} of the total nodes *n*, where *m* is the total number of discontinuities and *D*_{i} contains the vertices of discontinuity *i*. Discontinuity *i* has a boundary formed by a subset *K*_{i} of the total number of edges *l*. Each edge *j* (*j*=1,2,…,*l*) is defined by subset *E*_{j} of the total nodes *n*, where *E*_{j} contains the two vertices of edge *j*. The following is a suitable convention for selecting *k*_{ij} for all edges *j* (*j*=1,2,…,*m*) and all discontinuities *i* (*i*=1,2,…,*m*).

(i) For each discontinuity

*i*(*i*=1,2,…,*m*), define a positive ordering of the vertices in*D*_{i}. The two vertices of each edge*j*(*j*∈*K*_{i}) must be adjacent in this ordering. The ordering defined is cyclical (i.e. the first vertex is adjacent to the last) and is also used to define**n**_{i}(see appendix A).(ii) For each edge

*j*(*j*=1,2,…,*l*), define a particular ordering of the vertices*E*_{j}as positive.(iii) For each edge

*j*(*j*=1,2,…,*l*) and discontinuity*i*(*i*∈*S*_{j}) with common vertices {*X*,*Y*}=*D*_{i}∩*E*_{j}, compare the orderings of*X*and*Y*in the positive orderings of*E*_{j}and*D*_{i}defined in (i) and (ii), respectively, noting that the ordering of*D*_{i}in (i) is cyclical. If these orderings coincide,*k*_{ij}=1; otherwise,*k*_{ij}=−1.

## Appendix C. Crossovers and overlaps

**(a) Crossovers**

Intersections between discontinuities can occur, much as in plane strain DLO problems. The intersections between discontinuities on different planes are line segments or ‘crossovers’. Compatibility is implicitly ensured at these crossovers rather than being explicitly enforced. This can be demonstrated using the two intersecting rectangular discontinuities, as shown in figure 12*a*. In figure 12*b*, the discontinuities are split into smaller discontinuities so that the crossover is eliminated. In DLO, it is assumed that the displacement jump across each of these smaller discontinuities is equal to that of its parent, taking account of the conventions in appendices A and B. It is clear that equation (B1) is satisfied for the new edges and the original crossover itself must also be compatible.

**(b) Overlaps**

Intersections also occur between discontinuities in the same plane, resulting in ‘overlaps’. As with crossovers, compatibility is implicitly ensured. In figure 13*a*, two overlapping discontinuities are shown. These are divided in figure 13*b* so that the original overlap is eliminated, but two overlapping discontinuities and remain. Discontinuities located on different planes and sharing edges *AB*, *BC*, *EF* and *FD* are similarly divided (not shown in figure 13). The displacement jump for each new discontinuity equals that of its parent discontinuity, taking account of the conventions in appendices A and B. Using equation (B.1), it is simple to prove that compatibility is maintained and that the displacement jump **Δ****v**_{I} on the overlapping region *I*_{1}*I*_{3}*I*_{4}*I*_{2} equals **Δ****v**_{ABC}+**Δ****v**_{DEF}.

Calculating *p*_{ABC}, *p*_{DEF} and *p*_{I} from equation (3.11), it is clear that *p*_{ABC}+*p*_{DEF}≥*p*_{I} since these terms are not added vectorially in the SOCP solver. It can therefore be concluded that overlaps will overestimate both the dilation and energy dissipated. This is equivalent to the overlapping region having a cohesive strength *c* and angle of friction *ϕ* greater than those of the original discontinuities, thus maintaining the upper bound character of the solution. The solver will tend to avoid such overlaps where these result in more energy being dissipated than necessary. Therefore, the significance of these overlaps is likely to reduce with increasing nodal resolution.

The load distribution across discontinuities *ABC* and *DEF* is unknown and can be interpreted in any manner consistent with the work done by **Δ****v**_{ABC}, **Δ****v**_{DEF} and **Δ****v**_{I}. Note that the load on *I*_{1}*I*_{3}*I*_{4}*I*_{2}, *L*_{I}, cannot be the sum of the loads on and , and , respectively (i.e. ). Load distributions resulting in will always be consistent with the work done.

- Received January 4, 2013.
- Accepted March 28, 2013.

© 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.