Application of discontinuity layout optimization to three-dimensional plasticity problems

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.


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][2][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 threedimensional problems has been more limited.  where 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) p tan φ − n = 0 (2.3) and p ≥ s 2 + t 2 , (2.4) where 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 DLObased formulation that can be applied to three-dimensional problems involving Tresca and Mohr-Coulomb materials.

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], subject to and 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 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.  . 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).

(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 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 3b, the following must also be true along edge OO : 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 on July 20, 2018 http://rspa.royalsocietypublishing.org/ Downloaded from  where 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 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: where 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: and where 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: where 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 diag(tan φ 1 , tan φ 2 , . . . , tan φ m ). 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 i c da, 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.13) subject to and where f D and f L are vectors containing, respectively, the specified dead and live loads. Here d contains displacement jumps across the discontinuities, 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 3l × 3m 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 Li 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: where 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; -set-up problem, using equations (3.13)-(3.17); and -solve the resulting SOCP problem.

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.      [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. 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 Failure 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 where q is the average pressure on the anchor. Merifield et al. [23] used a lower bound finiteelement 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,

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.

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 v 1 and v 2 , respectively. The displacement jump v ABC across discontinuity ABC

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 (B 1) 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 13a, two overlapping discontinuities are shown. These are divided in figure 13b so that the original overlap is eliminated, but two overlapping discontinuities − −−− → I 1 I 3 I 4 I 2 ABC and − −−− → I 1 I 3 I 4 I 2 DEF 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