Royal Society Publishing

Application of discontinuity layout optimization to three-dimensional plasticity problems

Samuel Hawksbee , Colin Smith , Matthew Gilbert


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 [13]. 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) Embedded Image2.1and Embedded Image2.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) Embedded Image2.3and Embedded Image2.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.

Figure 1.

Mohr–Coulomb yield criteria: (a) yielding discontinuity between two rigid blocks showing normal force N, resultant shear force P and its components S and T; (b) conic yield surface where Embedded Image is the normal axis (tension positive) and Embedded Image and Embedded Image are orthogonal shear axes. The displacement jump orthogonal to the yield surface is shown, where n and p are the normal displacement and shear plastic multiplier, respectively.

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.

Figure 2.

Stages in DLO procedure, after [14]: (a) starting problem (surcharge applied to block of material close to a vertical edge); (b) discretization of domain using nodes; (c) interconnection of nodes with potential discontinuities; (d) identification of critical subset of potential discontinuities using optimization (giving the layout of slip-lines in the critical failure mechanism).

For plane strain problems, the corresponding optimization problem is given in equations (3.1)–(3.5), posed entirely in terms of potential discontinuities [11], Embedded Image3.1subject to Embedded Image3.2 Embedded Image3.3 Embedded Image3.4 and Embedded Image3.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 (fL) 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 (fD) and internal plastic energy dissipation gTp (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 2d.

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, vAB to vEA. The prisms are separated by rectangular discontinuities OOXX′, where X=A, BE. Each discontinuity OOXX′ has a normal nX and a relative displacement jump ΔvX 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 Embedded Image3.6

Figure 3.

Solid bodies meeting at a common edge, (a) prior to, and (b) after movement, with prisms separated by discontinuities Embedded Image, where X=A, BE. For clarity, O′ has not been shown, but lies in the plane ABCDE′.

A sign convention for determining the directions of nA to nE and ΔvA to ΔvE is presented in appendix A. Using this sign convention and the vertex ordering in figure 3b, the following must also be true along edge OO′: Embedded Image3.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 Embedded Image3.8where Ti is a 3×3 transformation matrix converting local to global displacement jumps. Ti is chosen such that one axis aligns with ni, the unit column vector in the normal direction using the ‘right-hand screw rule’; the two remaining orthogonal axes, si and ti, are in the plane of discontinuity i such that Ti={nisiti}. Also di={nisiti}T, where ni, si and ti are the local displacement jumps across discontinuity i in the ni, si and ti directions, respectively.

Compatibility can now be enforced for each edge j (j=1,2,…,l) as follows: Embedded Image3.9where Sj is a subset of the total number of discontinuities m that contains the discontinuities meeting at edge j, and Bij is a local compatibility matrix equal to kijTi, where kij 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: Embedded Image3.10and Embedded Image3.11where pi 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: Embedded Image3.12where dn 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 Embedded Image. 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.

View this table:
Table 1.

Discontinuity flow rule conditions: applicable constraints and variables.

(iii) Dissipation function

The energy dissipated on a given discontinuity i is given by gipi, where gi is a dissipation coefficient equal to Embedded Image, 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 gi=aici, where ai and ci 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: Embedded Image3.13subject to Embedded Image3.14 Embedded Image3.15 Embedded Image3.16 and Embedded Image3.17where fD and fL are vectors containing, respectively, the specified dead and live loads. Here d contains displacement jumps across the discontinuities, dT={n1,s1,t1,n2,s2,t2,…,tm}, where ni is the displacement jump normal to discontinuity i,si and ti are the displacement jumps within the plane of discontinuity i and g is a vector of dissipation coefficients. Here B is a suitable 3l×3m compatibility matrix, N is a suitable m×m flow matrix and dn (a subset of d) is a vector containing the normal displacement jumps, dTn={n1,n2,…,nm}. Also p is a vector of plastic multipliers, pT={p1,p2,…,pm}, where pi is the plastic multiplier for discontinuity i given by equation (3.17). The optimization variables are the displacement jumps in d (and dn) and the plastic multipliers in p. The objective function and the first three constraints are linear. The final constraints on the plastic multipliers pi 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 fTD and fTL in equations (3.13) and (3.16) are now defined such that Embedded Image and Embedded Image, where Embedded Image, Embedded Image, Embedded Image and Embedded Image, Embedded Image, Embedded Image are, respectively, the dead and live loads acting in the ni, si and ti 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 fDi and fLi 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 fDi and fLi 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: Embedded Image3.18where Wi is a 1×3 row vector containing the components in the ni, si and ti 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;

  • — set-up problem, using equations (3.13)–(3.17); and

  • — 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.

Figure 4.

Compression of a block: problem geometry, nodal spacings (Δx, Δy and Δz) and the 1/16 of the volume modelled, owing to symmetry.

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

View this table:
Table 2.

Compression of a block: comparison with benchmark solutions. Upper bound (UB) solutions from Martin & Makrodimopoulos [8] have been used to benchmark the present solutions, and their lower bound (LB) solutions are also listed for information; Δ=Δxyz.

Figure 5.

Compression of a block: typical failure mechanism for ϕ=0 case (Δx, Δy, Embedded Image).

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.

View this table:
Table 3.

Compression of a block: use of different nodal grid spacings. Upper bound (UB) solutions from Martin & Makrodimopoulos [8] have been used to benchmark the present solutions, and their lower bound (LB) solutions are also listed for information; Δx=1.

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

View this table:
Table 4.

Punch indentation: comparison with benchmark solutions. Upper bound (UB) solutions from Vicenta da Silva & Antão [21] have been used to benchmark the present solutions, and lower bound (LB) solutions from Salgado et al. [17] are also listed for information.

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 Nc=6.102 (difference 0.84%).

Figure 6.

Punch indentation: problem geometry, nodal spacing Δ, and the 1/8 of the volume modelled, owing to symmetry.

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.

Figure 7.

Punch indentation: representative failure mechanism (Δ=1/2, ϕ=0) (dashed lines indicate extent of domain modelled).

(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, Δxy=1, Δz=H/4 and varying H.

Figure 8.

Anchor: problem geometry, nodal spacings (Δx, Δy and Δz) and the 1/8 of the volume modelled, owing to symmetry.

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 Nc0 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 as Embedded Image4.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 Nc* is independent of the embedment depth H. For a given ratio γ/c, the deep mechanism becomes critical at depths greater than or equal to Hcr. If both deep and shallow mechanisms are considered, N must be less than or equal to Nc*. Merifield et al. [22] established a lower bound on Nc*≈11.9. For γ/c=0, Merifield et al. [22] found Hcr≈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<Hcr. 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 Nc and N were not influenced by the extent of the domain modelled.

Figure 9.

Anchor: cohesive soil break-out factors.

(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 Embedded Image4.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, Embedded Image4.3Using nodal grids of the form shown in figure 8 and fixing B=2, W=10, Δxy=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, Δxy=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).

Figure 10.

Anchor: frictional soil break-out factors using nodal grids GRID 1 and GRID 2.

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 11a shows two triangular prisms separated by a discontinuity ABC moving with absolute displacements v1 and v2, respectively. The displacement jump ΔvABC across discontinuity ABC must equal either v1v2 or v2v1. 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 11b,c for discontinuity ABC. ΔvABC 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.

Figure 11.

Sign convention: (a) two triangular prisms separated by a discontinuity ABC, the dashed outline shows the original positions of the prisms; (b) discontinuity vertex ordering Embedded Image; (c) discontinuity vertex ordering Embedded Image.

Appendix B. Convention for selecting the vertex ordering parameter kij

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: Embedded ImageB1and Embedded ImageB2where Sj is the subset of all discontinuities meeting at edge j. kij=±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 Di of the total nodes n, where m is the total number of discontinuities and Di contains the vertices of discontinuity i. Discontinuity i has a boundary formed by a subset Ki of the total number of edges l. Each edge j (j=1,2,…,l) is defined by subset Ej of the total nodes n, where Ej contains the two vertices of edge j. The following is a suitable convention for selecting kij 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 Di. The two vertices of each edge j (jKi) 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 ni (see appendix A).

  • (ii) For each edge j (j=1,2,…,l), define a particular ordering of the vertices Ej as positive.

  • (iii) For each edge j (j=1,2,…,l) and discontinuity i (iSj) with common vertices {X,Y }=DiEj, compare the orderings of X and Y in the positive orderings of Ej and Di defined in (i) and (ii), respectively, noting that the ordering of Di in (i) is cyclical. If these orderings coincide, kij=1; otherwise, kij=−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 12a. In figure 12b, 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.

Figure 12.

Crossovers: (a) discontinuity Embedded Image, displacement jump ΔvABCD, intersecting discontinuity Embedded Image, displacement jump ΔvEFGH; (b) Embedded Image split into Embedded Image and Embedded Image, both with a displacement jump ΔvABCD; Embedded Image split into Embedded Image, Embedded Image, Embedded Image and Embedded Image, all with a displacement jump ΔvEFGH.

(b) Overlaps

Intersections also occur between discontinuities in the same plane, resulting in ‘overlaps’. As with crossovers, compatibility is implicitly ensured. In figure 13a, two overlapping discontinuities are shown. These are divided in figure 13b so that the original overlap is eliminated, but two overlapping discontinuities Embedded Image and Embedded Image 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 ΔvI on the overlapping region I1I3I4I2 equals ΔvABC+ΔvDEF.

Figure 13.

Overlaps: (a) a discontinuity Embedded Image, displacement jump ΔvABC, overlapping a discontinuity Embedded Image, displacement jump ΔvDEF; (b) Embedded Image split into Embedded Image, Embedded Image and Embedded Image, each with a displacement jump ΔvABC; Embedded Image split into Embedded Image, Embedded Image and Embedded Image, each with a displacement jump ΔvDEF.

Calculating pABC, pDEF and pI from equation (3.11), it is clear that pABC+pDEFpI 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 ΔvABC, ΔvDEF and ΔvI. Note that the load on I1I3I4I2, LI, cannot be the sum of the loads on Embedded Image and Embedded Image, Embedded Image and Embedded Image, respectively (i.e. Embedded Image). Load distributions resulting in Embedded Image will always be consistent with the work done.

  • Received January 4, 2013.
  • Accepted March 28, 2013.
Creative Commons logo

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


View Abstract