Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 2, pp. 417-433, Warsaw 2018 DOI: 10.15632/jtam-pl.56.2.417 A NUMERICAL UPPER BOUND FORMULATION WITH SENSIBLY-ARRANGED VELOCITY DISCONTINUITIES AND ORTHOTROPIC MATERIAL STRENGTH BEHAVIOUR Mingjing Li, Josef Füssl, Markus Lukacevic, JOSEF EBERHARDSTEINER Vienna University of Technology, Vienna, Austria e-mail: mingjing.li@tuwien.ac.at Christopher M. Martin University of Oxford, Oxford, UK Numerical limit analysis allows for fast estimates of the collapse loadof structures exhibiting ideal plastic material behaviour. In numerical upper bound formulations, the description of the unknown velocity field can be extended by introducing velocity discontinuities between finite elements. Through these additional degrees of freedom, localised failure modes may be approximatedmore accurately and better upper bounds can be obtained. In the existing formulations, such discontinuities are typically introduced between all elements and the de- scription is restricted to isotropic failure behaviour. In this work, a general 3D upper bound formulation is briefly proposed, allowing the consideration of both isotropic and orthotropic yield functions within finite elements as well as at velocity discontinuities. The concept of “projecting” a stress-based orthotropic yield function onto a certain discontinuity is presen- ted, giving a traction-based yield function which allows for a consistent description of the material strength behaviour across the interface. The formulation is verified by means of two classical examples, the rigid strip footing and the block with asymmetric holes. Fur- thermore, based on the computation of potential orientations of plastic flow localisation, a simple concept for a sensible arrangement of velocity discontinuities is proposed. It is shown that this concept performs very well for isotropic as well as anisotropic material strength behaviour. A feature of the present work is that, velocity jumps are allowed only across the prescribed finite element interfaces determined from the sensible discontinuity arrangement. Goodupper bounds similar to those in the existingworksare obtainedwith far fewer degrees of freedom. Keywords:numerical upper bound formulations, localised failuremodes, traction-basedyield function, sensible arrangement of velocity discontinuities, orthotropic material strength be- haviour 1. Introduction 1.1. Numerical limit analysis Themain objective of limit analysis is the determination of load bearing capacities of struc- tures exhibiting an elastoplastic material response. To achieve this, limit analysis concentrates on the critical energy dissipation rate at the time instant of structural failure, and the basic task can be expressed as follows: Find the kinematically admissible velocity field which minimi- ses energy dissipation over the set of all statically admissible stress fields which maximise the dissipated energy (Ciria et al., 2008). Statically admissible stress fields are required to be in equilibrium, fulfil the static boundary conditions and obey a plastic yield criterion at each point of the body. Kinematically admissible velocity fields are subject to compatibility, the kinematic 418 M. Li et al. boundary conditions, and fulfil an associated plastic flow rule at each point of the body.Unfortu- nately, the so-defined saddle-point problem can be solved only for simple geometric and loading situations as well as for simplematerial behaviours. Inmore complex situations, the plastic flow compatibility in the so-called static principle or the static equilibrium in the so-called kinematic principlemay be relaxed, providing lower and upper bounds on the exact load bearing capacity of a structure according to the bounding theorems by Drucker et al. (1951, 1952). However, for complex problems, the application of these bounding theorems (in the context of limit analysis) in an analytical way is very limited and often not possible. Thus, finite- -element-based formulations were first introduced in the 1970s (Belytschko and Hodge, 1970; Lysmer, 1970; Anderheggen andKnöpfel, 1972; Maier et al., 1972), and gained popularity from then on. The computational efficiency and accuracy of such numerical formulations strongly depend on the mathematical programming method used to solve the underlying optimisation problems. At the early stage, the limit analysis theorems were formulated as linear optimisation problems, by linearising the applied plastic yield functions. At the turn of the millennium, LyaminandSloan (2002a,b) proposedmoregeneral lower andupperbound formulations allowing for nonlinearyield functions,whichwere solvedusingnonlinearprogrammingconcepts.However, local smoothingwas required for yield functionswith singularities, e.g. theMohr-Coulomb yield function. Subsequently, during the past two decades, second-order cone programming (SOCP) has been proven to be an excellent alternative method byMakrodimopoulos andMartin (2006, 2007) and Ciria and Peraire (2004) for cohesive-frictional materials and Füssl et al. (2008) for composite materials, with sufficient robustness and efficiency to solve large-scale nonlinear optimisation problems of numerical limit analysis. Such implementations allow the applications of many different yield functions in their native form, since most of the commonly-used yield functions can be formulated as second-order cones. In this work, SOCP is employed to solve the nonlinear optimisation problems arising from the presented limit analysis formulations. The efficiency and accuracy of such formulations is also strongly influenced by the chosen finite elements and related shape functions. Inorder toobtain rigorousupperboundsolutions, for example, the associated plastic flow rulemust be satisfied throughout thewhole body.Basically, this can be achieved byusing constant strain triangular elements, which are often combinedwith velocity discontinuities between element boundaries (Bottero et al., 1980; Sloan and Kleeman, 1995; Lyamin and Sloan, 2002b). To improve the quality of upper bound solutions, the use of higher order interpolation functions is desired. Makrodimopoulos and Martin (2007) showed that the associated plastic flow rule could also be enforced throughout the whole body by using linear strain triangular elements, leading to a better performance than constant strain elements even without discontinuities. As a further development, the meshless method was implemented for numerical upper bound approaches in (Le et al., 2010; Liu and Zhao, 2013; Yu et al., 2016). However, in such implementations with high order shape functions, it is difficult to guarantee both compatibility and satisfaction of the associated plastic flow rule throughout each element. Additionally or as an alternative to the use of high order elements, velocity discontinuities can be implemented in upper bound formulations to increase their effectiveness. In (Chen et al., 2003; Milani and Lourenço, 2009), for example, rigid elements were used and plastic dissipation was only allowed between finite elements. Such approaches are highly dependent on the mesh and even adaptivemesh refinement cannot fully compensate for this issue. An approachwithout using classical finite elements is the so-called discontinuity layout optimisation (DLO), where velocity discontinuities are determined by using a truss layout optimisation algorithm based on a prescribed grid (Smith andGilbert, 2007; Hawksbee et al., 2013). This approach performswell for 2D problems, but the determination of complex failure mechanisms in 3D bodies requires a fine grid and large computational effort. For this reason, in theauthors’ opinion, themostpromisingapproach so far toobtain rigorous upper bound solutions still seems to be the use of solid finite elements with or without velocity A numerical upper bound formulation with sensibly-arranged velocity... 419 discontinuities. In Krabbenhøft et al. (2005) zero-thickness interface elements between constant strain elements are introduced, which performwell for a large number of applications. Another development canbe found inMakrodimopulousandMartin (2008),wherevelocity discontinuities are implemented between linear strain elements. In order to increase the efficiency of the upper bound formulations, adaptivemesh refinementwas introducedbyCiria andPeraire (2004), Ciria et al. (2008) and Martin (2011). However, a targeted arrangement of discontinuities, as will be proposed in this work, has not been introduced until now. 1.2. Objective of the paper In several previousworks, e.g. in (Füssl et al., 2017;Li et al., 2018), anisotropic yield functions have been implemented in numerical upper bound formulations. To the authors’ knowledge, the combination of anisotropic yield functions and velocity discontinuities has not previously been presented, although it could significantly improve the capability of upper bound approaches in handling localised plastic failure for anisotropic materials, like wood or fibre reinforced compo- sites. In particular, it is beneficial if the alignment of discontinuities is tuned to the direction of localised plastic failure. Thus, the main objectives of this work can be introduced as follows: 1. The formulation of 3Dnumerical upperboundapproacheswith anisotropic yield functions, quadratic shape functions for the velocity fields, and velocity discontinuities. 2. To allow for a consistent description of plastic failure also across velocity discontinuities, the derivation of a traction-based yield function which is in accordance with the stress- -based yield function assigned to the solid finite elements/bulk material. 3. Implementationof an initial concept for a sensible introductionandarrangementof velocity discontinuities only in failure regions. According to these objectives, thepaper is structuredas follows.Aquite general numerical upper bound approach is briefly proposed in Section 2, able to consider plastic energy dissipation in both finite elements and discontinuities obeying an anisotropic failure criterion. Furthermore, the process for obtaining the required traction-based yield functions for the discontinuities is described. A verification of the implemented upper bound formulations bymeans of well-known examples can be found in Section 3, as well as a discussion about the performance of velocity discontinuities. Finally, a brief summary and concluding comments are given in Section 4. 2. Upper bound approaches The upper bound theorem focuses exclusively on the kinematically admissible velocity fields u̇ = (u̇x, u̇y, u̇z) T ∈ R3, and by minimising the internal plastic energy dissipation rate Wint, which has to be equal to the work rate of the external loads Wext, the resulting failure sta- te provides an upper bound for the exact collapse load. A kinematically admissible velocity field u̇ has to satisfy the compatibility, the associated plastic flow rule, and the kinematic bo- undary conditions at each point of the considered body. Additionally, a velocity-jump field ∆u̇=(∆u̇x,∆u̇y,∆u̇z) T ∈ R3 is introduced, describing localised interface plastic failure across an prescribed interior surface. The internal energy dissipation rateWint is composed of a part referring to material failure in the continuum bodyΩ and a part related to the energy dissipation at interior surfaces Γdis, and reads Wint = ∫ Ω dmatp (ε̇) dV + ∫ Γdis ddisp (∆u̇) dA (2.1) 420 M. Li et al. with the plastic dissipation functions dmatp = sup σ∈F σ T ε̇ F = {σ| f(σ)¬ 0} in Ω ddisp =sup t∈D tT∆u̇ D= {t| f(t)¬ 0} on Γdis (2.2) where ε̇ = (ε̇xx, ε̇yy, ε̇zz, ε̇xy, ε̇yz, ε̇xz) T ∈ R6 represents the plastic strain-rate field, σ = (σxx,σyy,σzz,τxy,τyz,τxz) T ∈ R6 the stress field and t = (tx, ty, tz)T ∈ R3 the surface traction field. f(σ) ¬ 0 and f(t) ¬ 0 denote the stress-based yield function for Ω and the traction-based yield function for Γdis, respectively. The upper bound theorem can then be formulated as a nonlinear optimisation problem, reading min Wint s.t. ε̇= divu̇ in Ω u̇= u̇b on Γ ε̇= λ̇σ∂f(σ)/∂σ in Ω ∆u̇= λ̇t∂f(t)/∂t on Γ dis (2.3) inwhich the constraints enforce compatibility between the velocities and the plastic strain-rates, the kinematic boundary conditions, and the associated plastic flow rule both in the continuumΩ and at the interior surfaces Γdis. In the second constraint, u̇b refers to the prescribed velocity boundary conditions defined over the whole surface Γ = ∂Ω of the continuum body. In the last two constraints, λ̇σ and λ̇t are plastic multipliers, determining the magnitude of plastic flow within the continuum and at the discontinuities, respectively. Note that these two associated plastic flow constraints in Eq. (2.3) are valid only when the yield function is differentiable everywhere. If singular apexpoints exist, additional technology is required, and the use of SOCP in this work ensures that such points are handled naturally. As shown in Makrodimopoulos and Martin (2007), using the duality of nonlinear program- ming, amathematically equivalent optimisation problem toEq. (2.3) can be formulated, reading max Wext s.t. ∫ Ω (divu̇)Tσ dV + ∫ Γdis ∆u̇Tt dA= ∫ Ω u̇Tβg dV + ∫ Γ u̇Tβt dA in Ω f(σ)¬ 0 in Ω f(t)¬ 0 on Γdis (2.4) in which the first constraint represents weak equilibrium of the dissipated energy, and the ob- jective function is related to the external work rate given as Wext = ∫ Ω u̇ Tβg dV + ∫ Γ u̇ Tβt dA (2.5) whereβ denotes a loadmultiplier applied to the surface traction field t and the prescribed body force field g∈ R3. For the discretisation of the upper bound optimisation problem, tetrahedral linear-strain simplex elements are used, as introduced for the 2D upper bound problem under plane strain conditions inMakrodimopoulos andMartin (2007, 2008) and for the 3Dupperboundproblem in Martin and Makrodimopoulos (2008). Thus, the velocity field is approximated using quadratic interpolation functions and the plastic strain-rate field is described by linear shape functions. Worthmentioning is that each finite element has its own strain-rate evaluation nodes,whichme- ans that adjacentnodes fromdifferent elements share the samecoordinates but canhavedifferent A numerical upper bound formulation with sensibly-arranged velocity... 421 strain-rate states. The exact representation of this approximation is given inMakrodimopoulos and Martin (2008). In this work, to assess the capability of the discontinuity arrangement, ve- locity jumps are allowed only across particular prescribed finite element interfaces determined by the arrangement, on which adjacent elements have their own velocity evaluation nodes. Finally, as introduced in detail by the authors in Li et al. (2018), the discretised formulation of the dual upper bound optimisation problem, Eq. (2.4), can be written as max Wext s.t. AmatUB T q̂σ+A dis UBq̂ dis t =βA bc UBq̂ bc t ŝmat,iσ = â mat,i σ +B mat σ R mat ε̇ q̂ mat,i σ ŝmat,iσ ∈ C ŝ dis,j t = â dis,j t +B dis t R dis t q̂ dis,j t ŝ dis,j t ∈ C (2.6) with the matrices AUB obtained by applying the linear compatibility operator to the related shape functions of the velocity (or the velocity jump), within the finite elements (mat), at the discontinuities between elements (dis), and at the boundary (bc). The vectors q̂σ, q̂ dis t , and q̂ bc t collect all nodal degrees of freedom related to stress-like quantities and surface tractions. Note that, for an arbitrary vector x, the symbol x̂= ∫ Ωe x dV refers to the volume-integrated quanti- ty over each element. The remaining constraints represent a second-order cone formulation of a general quadratic yield function for both the solid material (mat) and the discontinuities (dis), with i and j ranging from 1 to the number of stress and traction evaluation nodes, respective- ly. Thereby, the matrices â and B contain strength parameters and the matrices R represent transformation operators, rotating the stress tensors into the principal material direction and the surface traction vector into the direction of the corresponding discontinuity. The external work rate in discretised form can be written as Wext = UBC∑ bc=1 6∑ i=1 βq bc,i u̇,loc T q̂ bc,i t,loc (2.7) where UBC denotes the number of 6-noded boundary surface triangular elements with a pre- scribed local traction field q̂bct,loc and q bc u̇,loc represents the related velocity degrees of freedom. In previous upper bound formulations (Sloan andKleeman, 1995; Krabbenhoft et al., 2005; Makrodimopoulos and Martin, 2008) only isotropic yield functions based on shear failure me- chanisms, e.g. the vonMises or theMohr-Coulomb yield function, were considered, leading to a straightforward definition of failure at the discontinuities between elements. On the contrary, in the upper bound formulation Eq. (2.6), the quite general orthotropic yield function according to Tsai-Wu can be implemented, reading q iT Pq i+ (1 2 F +T q i )2 − ( 1− 1 2 F −T q i )2 ¬ 0 (2.8) with i as the evaluation point of q either for the stress field (with subscript σ) in the element or for the traction field (with subscript t) at a discontinuity. The vectors F+,F− andmatrixP are related to the terms in Eq. (2.6) as follows a= [ 1 0 ] B=    − 1 2 F−T D 1 2 F+T    (2.9) 422 M. Li et al. whereD is thedecomposedproductofP=DTD.Note that thematrixdimensionsareF+σ ∈ R6, F−σ ∈ R6, Pσ ∈ R6×6, Dσ ∈ R6×6 for the stress-based yield function and F+t ∈ R3, F−t ∈ R3, Pt ∈ R3×3,Dt ∈ R3×3 for the traction-based yield function. In theabove, it is assumed that there exists a traction-based yield function for thediscontinu- ities which can also be formulated as a second-order cone. Additionally, it needs to be consistent with the stress-based Tsai-Wu criterion to allow for the description of a homogeneous strength distribution within a body. Since a surface traction state within a discontinuity cannot directly be related to a unique 3D stress state at a material point, the derivation of such a traction- based yield function is not straightforward. However, according to Wu and Cervera (2014), we can assume a 3D plastic strain-rate state to be localised with respect to a certain discontinuity if the following constraints are satisfied, for orthotropic yield functions which are differentiable everywhere (e.g. Tsai-Wu yield function): Λdismm(σ dis loc)= ε̇mm λ̇σ = ∂f(σdisloc,strength par.) ∂σdismm =0 Λdispp (σ dis loc)= ε̇pp λ̇σ = ∂f(σdisloc,strength par.) ∂σdispp =0 Λdismp(σ dis loc)= ε̇mp λ̇σ = ∂f(σdisloc,strength par.) ∂τdismp =0 (2.10) where σdisloc denotes a 3D stress state at a discontinuity with the local coordinate basis (n-m-p) with the normal vector of the discontinuity pointing in the n-direction. Note that each of the constraints in Eq. (2.10) can be formulated as a function of the local stress field σdisloc. By reformulating Eq. (2.10) using the definition of the Cauchy stress tensor, giving tdisn = σ dis nn , tdism = τ dis nm, t dis p = τ dis np , the three remaining stress tensor components σ dis mm, σ dis pp , and τ dis mp can be expressed as functions, hereafter referred to as Ldis, of tdisn , t dis m , t dis p , and certain strength parameters, reading σdismm =L dis mm(t dis loc,strength par.) σ dis pp =L dis pp (t dis loc,strength par.) σdismp =L dis mp(t dis loc,strength par.) (2.11) Therefrom, a relationship between the local stress field and the local traction field on Γdis, σ dis loc = L dis t t dis loc, under the condition of plastic strain localisation, can be derived. Finally, by making use of this relationship, it is possible to “project” the stress-based formulation of the Tsai-Wu yield function f(σdisloc)¬ 0 onto a discontinuity, delivering a consistent traction-based yield function f(tdisloc)¬ 0. Themain focus of this work is the performance assessment of this approach and to point out how such an approach could be utilised for future concepts of numerical limit analysis. For this reason, in the next Section, several numerical examples are presented and discussed in detail. 3. Numerical results In this Section, numerical results obtained using the proposed upper bound formulationwith se- lectively activated velocity discontinuities are discussed.Twobenchmarkproblemswith isotropic yield functions are used for basic verification of the presented approaches. By means of further examples, it is demonstrated that orthotropic plastic failure can also be handled appropriate- ly. Moreover, it will be shown that through the introduction of velocity discontinuities across properly-arranged prescribed interfaces, high-quality upper bound results can be obtained with relatively coarsemeshes and, thus, very efficiently. Note that, for convenience, the upper bound A numerical upper bound formulation with sensibly-arranged velocity... 423 results obtained using formulations, with andwithout velocity discontinuities, are referred to as continuous and discontinuous upper bound results, respectively. All computations presented in the following have been performed on a Linux desktopmachi- ne with an AMD Phenom(tm) II X6 1090T CPU (6 cores) and 8GB of RAM. The commercial software package Abaqus was used for mesh generation, but all other pre- and post-processing tasks as well as the assembly of the SOCP optimisation problems were carried out by self- -written codes in Fortran. The SOCP optimisation problems themselves were solved using the commercial software MOSEK (2014), which is based on the conic interior-point algorithm de- scribed in Andersen et al. (2003). 3.1. Rigid strip footing The rigid strip footing problem, as illustrated in Fig. 1a, with a weightless purely cohesive material is a commonbenchmark for limit analysis approaches.Upon the assumption ofmaterial failure according to Tresca, τ ¬ c, the ultimate load can be obtained by the classical Prandtl solutionNrefc =Plim/c=2+π (Prandtl, 1920), whereNc is the bearing capacity factor,Plim the collapse load limit, τ the principal shear stress, and c the coefficient of cohesion. Under plane strain conditions, the Tresca yield function is identical to the vonMises yield function √ J2 ¬ c, with J2 as the second deviatoric stress invariant. Fig. 1. Rigid strip footing benchmark example: (a) geometry and boundary conditions; (b) example 3D model and illustrative discretisation with 419 elements; (c) the prescribed interfaces for velocity jumps according to Prandtl’s failure mechanism (445 elements) Using the yield function formulation according toEq. (2.8), the vonMises criterion is defined through P mat σ = 1 3c2    1 −0.5 −0.5 0 0 0 −0.5 1 −0.5 0 0 0 −0.5 −0.5 1 0 0 0 0 0 0 3 0 0 0 0 0 0 3 0 0 0 0 0 0 3    F mat,+ σ =F mat,− σ =0 (3.1) A consistent traction-based yield function is easily obtained (Makrodimopoulos and Martin, 2008) by applying the shear strength as the tangential strength at discontinuities, giving P dis t = 1 c2    0 0 0 0 1 0 0 0 1    Fdis,+t =F dis,− t =0 (3.2) The geometric boundary conditions and loading are given in Fig. 1a, and an example 3D representation of the model is plotted in Fig. 1b with respect to the global coordinate basis (xyz) and an illustrative discretisation using 419 tetrahedron elements. By applying symmetric boundary conditions at the z− and z+ boundary surfaces, plane strain conditions are enforced. 424 M. Li et al. Thus, in the following, all resultswill beplotted in thexy-plane only.The rough footing interface condition is applied by setting the velocities in thex-direction to zero for all nodes in the footing region. The obtained numerical upper bound results for different fineness of discretisation or degrees of freedom (DOF) are plotted in Fig. 2a. The black curve represents the results for continuous velocity fields, and shows clear convergence behaviour. Measuring the difference between the upper bound results and the analytical referenceNrefc , as diff (%)= (N ub c −Nrefc ) ·100/(Nubc + Nrefc ), a diff of 12.92% (2790 DOF) for point a and 0.95% for point c with 188685 DOF are obtained. The corresponding CPU times are 2.09 min versus 41.85 min. Fig. 2. Numerical upper bound results for the rigid strip footing problem: (a) load limit factorN c obtained using continuous and discontinuous velocity fields as a function of DOF; (b) upper bound failure mode using discontinuities and 445 elements (range of plotted yield function values (robjv)[−3 ·10−8 : 0]); (c) upper bound failure mode using continuous velocity field with 31081 elements (robjv [−5 ·10−9 : 0]) In the next step, partitions are introduced into the model according to Prandtl’s failure mechanism (Prandtl, 1920), see Fig. 1c, to allow velocity jumpsacross prescribeddiscontinuities. Since velocity jumps are only allowed through these prescribed interfaces, the number ofDOF is not increased significantly. Theblue curve in Fig. 2a shows the related results, again for different numbers of DOF. A diff of 1.60% is obtained for point b (3147 DOF) and 0.13% for point d (196077 DOF), with corresponding CPU times of 0.56min and 41.62min. With a minimum difference of below 1% for the obtained best upper bounds compared to the analytical solution, the proposed formulations withstand this basic verification and can be assessed as performing well. Especially when velocity discontinuities are introduced, very good upper bounds can be obtained even with coarse meshes (see point b and the related failure mechanism in Fig. 2b). Of course, this is only possible if the failure mechanism is known in advance and, thus, does not yet represent an added value for general calculations. However, the potential of velocity discontinuities to capture very localised failure is evidently huge and, sensibly used, can greatly increase computational efficiency. 3.2. Block with asymmetric holes The block with asymmetric holes under tensile loading is a commonly-used benchmark for so-called direct methods (e.g. limit analysis and shakedown analysis) firstly studied by Zouain et al. (2002), as illustrated in Fig. 3a. Later on, this problemwas studied byMakrodimopoulos andMartin (2007) using an upper bound formulationwith a continuous quadratic velocity field, and the plane strain Mohr-Coulomb yield function was applied to the material. Their results are used for verification and as the reference solution in the following. Moreover, the benefit A numerical upper bound formulation with sensibly-arranged velocity... 425 of using sensibly-arranged velocity discontinuities is further discussed, and a simple strategy to find such arrangements based on preliminary upper bound results is proposed. Geometry, material properties, and boundary conditions (see Fig. 3a) are assigned as in the reference (Makrodimopoulos andMartin, 2007). Themodel is built by 3D finite elements with the global coordinate basis (xyz), similarly to Fig. 1b, and again symmetric boundary conditions at the z− and z+boundary surfaces are applied. Fig. 3. Block with asymmetric holes using the Drucker-Prager failure criterion: (a) geometry and boundary conditions; (b) upper bound failure mode (519 elements, robjv [−5 ·10−6 : 0]); (c) and (d) potential orientations of plastic strain localisation; (e) prescribed discontinuities based on (c) and (d) 3.2.1. Mohr-Coulomb failure criterion In Makrodimopoulos and Martin (2007), the Mohr-Coulomb yield function was used with a friction angle φ = 30◦. In our example, the equivalent Drucker-Prager yield function,√ J2 ¬A−BI1, is assigned to the material, defined by inserting P mat σ = 1 3A2    1 −0.5 −0.5 0 0 0 −0.5 1 −0.5 0 0 0 −0.5 −0.5 1 0 0 0 0 0 0 3 0 0 0 0 0 0 3 0 0 0 0 0 0 3    F mat,− σ = 1 A2    2AB 2AB 2AB 0 0 0    F mat,+ σ =0 (3.3) into Eq. (2.8), with the strength parameters A = 0.8321c and B = 0.1601. For the consistent traction-based yield function these terms read P dis t = 1 c2    −tan2φ 0 0 0 1 0 0 0 1    Fdis,+t =0 F dis,− t = 1 c2    2ctanφ 0 0    (3.4) In the previous example, it has been shown that velocity discontinuities can improve the upper bound significantly if they are appropriately arranged with respect to potential failure surfaces. For this reason, subsequently, the following strategy is pursued. Based on an upper bound calculation with a continuous velocity field and a relatively coarse mesh, as shown in 426 M. Li et al. Fig. 3b, potential discontinuity directions, where plastic failure could localise, are computed at each integration point where plastic flow takes place according to Eq. (2.10). The resulting orientations are plotted in Fig. 3c and Fig. 3d. With respect to the x axis, the mean orienta- tions obtained by taking the volume average over all orientations are approximately θ̄1 =−60◦ and θ̄2 = 60 ◦, as would be expected from the underlying failure criterion. Next, according to these average directions of possible localised plastic failure, and by referring to the points with maximum plastic strain-rates, the partitions (blue lines) shown in Fig. 3e are implemented into themodel. Finally, themodel is re-meshed and velocity discontinuities are introduced along the partitions. This procedure has been applied to several models with different level of mesh refinement (DOF) and compared to calculations performed without such discontinuities. The obtained upper bounds for the limit load factor Nc are plotted in Fig. 4a. Based on the upper bound result with only 3327 DOF (point a), velocity discontinuities were introduced into the model, leading to a strong improvement of the upper bound (point b) by adding only 585 DOF. The failure mechanism associated with point b is shown in Fig. 4b. To achieve a similarly good upper bound result and the related localised failure mechanism without the introduction of discontinuities, themesh needed to be refined significantly, as shown in Fig. 4c as the associated failure mechanism to point c. The CPU time required to obtain point c and the mechanism in Fig. 4c was 37.87min, whereas point bwith themechanism shown in Fig. 4b only took 0.38min. Although an adaptivemesh refinementwould probably bemore efficient than the uniformmesh refinement performed, the performance of sensibly-arranged discontinuities, even in a coarsely discretisedmodel, is excellent. Comparing the upper bound result of point b (Nc =1.062) to the best result in the reference (Makrodimopoulos and Martin, 2007) (Nc =1.063) there is almost no difference. However, the upper bound indicated by point b was obtained with 3912 DOF compared to 79955 DOF in the reference. Fig. 4. Numerical upper bound results for the block with asymmetric holes example using the Drucker-Prager yield function: (a) load limit factorN c obtained using continuous and discontinuous velocity fields as a function of DOF; upper bound results (robjv [−1 ·10−4 : 0]) (b) with discontinuities using 504 elements and (c) with continuous velocity field using 34875 elements In the following, this approach is extended to orthotropic material behaviour, using the Tsai-Wu criterion as the indicator for plastic failure. 3.2.2. Tsai-Wu failure criterion The strength parameters used for the example correspond to the typical orthotropicmaterial spruce wood and are taken from (Dorn, 2012), leading to A numerical upper bound formulation with sensibly-arranged velocity... 427 P mat σ =    2.434E-4 0 0 0 0 0 0 6.588E-2 0 0 0 0 0 0 6.588E-2 0 0 0 0 0 0 1.181E-2 0 0 0 0 0 0 2.973E-2 0 0 0 0 0 0 1.181E-2    MPa−2 F mat,+ σ =F mat,− σ = [ −6.573E-3 8.564E-2 8.564E-2 0 0 0 ]T MPa−1 (3.5) inserted into the general yield function in Eq. (2.8), definedwith respect to the local coordinate basis (LTR). Due to the orthotropic characteristics, the consistent traction-based yield func- tion is highly dependent on the orientation of the discontinuity at which it describes plastic failure. Thus, for each introduced discontinuity, a different set of strength parameters had to be computed, determined by the “projection” procedure introduced at the end of Section 2. Then, the equivalence between the material strength within elements and the strength behaviour at discontinuities is guaranteed. In the following, this orthotropic strengthbehaviour is assigned to theblockwith asymmetric holes, for two different material orientations, to assess the capability of the presented approach in handling anisotropic strength behaviour. Material orientation 1 In the first case, the longitudinal orientationL of thematerial is identical to the y-direction, as indicated in Fig. 5a. Moreover, the local material coordinate basis (LTR) is defined with L and T as the in-plane axes andR as the out-of-plane axis. Fig. 5. Block with asymmetric holes using the Tsai-Wu orthotropic failure criterion: (a) geometry, boundary conditions, and the principal material orientation indicated by the blue arrow; (b) upper bound result using 519 elements (robjv [−1.5 ·10−6 : 0]); (c) and (d) potential orientations of plastic strain localisation; (e) prescribed discontinuities based on (c) and (d) Thepreliminary upper bound calculation, based onwhich the arrangement of discontinuities will be defined, is carried out using a coarse mesh with 519 elements. The corresponding failure mechanism is displayed in Fig. 5b, and the computed orientations of potential discontinuities are shown in Figs. 5c and 5d.Averaging over these orientations results in twomean angle values of θ̄1 = −84◦ and θ̄2 = 84◦ with respect to the x-axis. According to these mean angles and starting from the points of highest plastic strain-rates (present at the boundaries of the holes), the velocity discontinuities shown in Fig. 5e are implemented into themodel. Subsequently, the model is re-meshed and a further upper bound calculation is performed. 428 M. Li et al. Fig. 6. Numerical upper bound results for the block with asymmetric holes with the principal material direction as defined in Fig. 5a: (a) collapse load using continuous and discontinuous velocity fields as a function of DOF; upper bound failure mode (robjv [−9 ·10−6 : 0]) using (b) discontinuities within a coarsemesh (486 elements) and (c) a continuous velocity field with a fine mesh (38537 elements) The numerical upper bounds on the collapse load for different levels of mesh refinement (DOF)with andwithout discontinuities are plotted inFig. 6.Again, the introduction of velocity discontinuities improves the upper bound significantly, while hardly increasing the DOF. In contrast to the isotropic casebefore, plastic failure also occurs in the solidfinite elements between discontinuities (see Fig. 6b), indicating that the orientation or arrangement of discontinuities could be improved. This can also be seen by looking at the plastic regions in Fig. 6c, which do not exactly match the definition of discontinuities above. It seems that the arrangement of discontinuities doesnotnecessarily have tobe ideal in order to achieve considerable improvement in the numerical upper bound results. Material orientation 2 In the second case, the local material orientation basis (LTR) is rotated by 30◦ in the xy- plane, as indicated in Fig. 7a. Again, based on an efficient preliminary upper bound calculation (see Fig. 7b), possible orientations of discontinuities are computed, resulting in mean angles of θ̄1 =±80◦ and θ̄2 =−58◦. The introduced velocity discontinuities are shown in Fig. 7e. Fig. 7. Block with asymmetric holes using the Tsai-Wu orthotropic failure criterion: (a) geometry, boundary conditions, and the principal material orientation indicated by the blue arrow; (b) upper bound result using 519 elements (robjv [−3 ·10−6 : 0]); (c) and (d) potential orientations of plastic strain localisation; (e) prescribed discontinuities based on (c) and (d) A numerical upper bound formulation with sensibly-arranged velocity... 429 Fig. 8. Numerical upper bound results for the block with asymmetric holes with the principal material direction as defined in Fig. 7a: (a) collapse load using continuous and discontinuous velocity fields as a function of DOF; upper bound failure mode (robjv [−2 ·10−5 : 0]) using (b) discontinuities within a coarsemesh (524 elements) and (c) a continuous velocity field with a fine mesh (34875 elements) As before, all numerical upper bounds are plotted in Fig. 8, showing the strong performance increase achieved by the selectively introduced discontinuities. Interestingly, the intensity of localisation of plastic failure is slightly different comparing the approaches with (Fig. 8b) and without (Fig. 8c) discontinuities. It seems that with discontinuities the real failure mechanism canbebetter represented, since themainplastic failuredirection iswell alignedwith theprincipal material direction (wood fibre direction) as expected. 3.3. Shear test on block The last example is used to assess the capability of this upper bound formulation regarding localised shear failure in orthotropic materials. The setup of the model is shown in Fig. 9a and is designed to represent a characteristic as well as often critical loading states in wood- -based products, like glued-laminated timber and cross-laminated timber. As in the previous examples, the boundary conditions are chosen so as to represent plane strain conditions, and the Tsai-Wu failure criterion with strength parameters representing spruce wood is assigned to the material. Again, based on a preliminary upper bound calculation using a very coarse mesh (Fig. 9c1), the orientations of possible discontinuities are determined (Fig. 9c2) and, based on that informa- tion, velocity discontinuities are implemented into themodel (Fig. 9c3). The result obtainedwith this discontinuity arrangement is shown in Fig. 9d1, where the dominant plastic failure appears within the solid elements between the introduced discontinuities but not at the discontinuities themselves, and, thus, the potential of the velocity discontinuities has not been activated suf- ficiently. For this reason, a second iteration was carried out, again computing the orientations of possible discontinuities in all plastic regions (Fig. 9d2). Based on that information, a revised discontinuity pattern was implemented as shown in Fig. 9d3. The corresponding failure mecha- nism is displayed in Fig. 10b and it can be seen that, now, localised failure occurs exclusively at the last-introduced discontinuity. The related upper bound (point b in Fig. 10a) is very good in comparison to the preliminarymodelwithout discontinuities (pointa) but uses a similar number of DOF.Moreover, the failuremechanism agrees well with that obtained using a very finemesh and no discontinuities (Fig. 10c). 430 M. Li et al. Fig. 9. Block under shear loading: (a) geometry, boundary conditions and the principal material orientation indicated by the blue arrow; (b) illustrative 3Dmodel with 222 finite elements; (c) and (d) two iteration steps for the definition of the discontinuities, with (c1) and (d1) as the upper bound result (robjv [−2 ·10−9 : 0]), (c2) and (d2) as the potential orientations of plastic strain localisation, and (c3) and (d3) as the prescribed discontinuities based on (c2) and (d2) Fig. 10. Numerical upper bound results for the block under shear loading: (a) collapse load using continuous and discontinuous velocity fields as a function of DOF; the upper bound failure mode (robjv [−4 ·10−9 : 0]) using (b) discontinuities within a coarsemesh (242 elements) and (c) a continuous velocity field with a fine mesh (32062 elements) 4. Summary and conclusions In this work, a 3D numerical upper bound formulation using a quadratic approximation of the velocity field and allowing for the implementation of orthotropic yield functions has been briefly proposed. Furthermore, the implementation of velocity discontinuities into this formulation has been presented, along with the concept of how to derive the necessary traction-based yield function which guarantees a consistent description of the material strength behaviour across discontinuities. Based on that formulation, compatibility as well as the associated plastic flow rule are fulfilled throughout the whole body and, thus, rigorous upper bounds are obtained. A numerical upper bound formulation with sensibly-arranged velocity... 431 The formulation has been verified bymeans of two classical examples, the rigid strip footing and the block with asymmetric holes. Subsequently, based on preliminary upper bound calcula- tions with very coarse meshes, the orientation of potential slip lines where very localised plastic flow may occur could be determined. Based on that information, sensibly-arranged velocity di- scontinuities were incorporated into the finite element models. In that way, upper bounds could be improved significantly with essentially no increase of degrees of freedom. It has also been shown that this concept works very well when assigning orthotropic failure behaviour to thematerial. If the first introduction of discontinuities does not improve the upper bound significantly, which means that most of the plastic dissipation still takes place in solid elements, a second iteration step can improve the situation, as shown by means of the shear test on a block. This may represent an important finding for future developments, which could lead to a general algorithm for an adaptive introduction and re-arrangement of velocity discon- tinuities, as an efficient alternative to existing adaptive mesh refinement strategies. Especially for laminated structures and orthotropic materials, where plastic failure often occurs in a very localised mechanism, as shown for wood at the microscale in (Lukacevic et al., 2014; Lukacevic and Füssl, 2016) and at the product scale in (Hochreiner et al., 2013, 2014), such an approach could have great value. Moreover, the efficiency of numerical limit analysis in combination with the accuracy of the extended finite element formulations presented in (Lukacevic et al., 2014; Lukacevic and Füssl, 2016) could lead to more flexible engineering design tools, in which the focus can be switched between accuracy and efficiency as needed. Acknowledgments Financial support for thiswork in the frameworkof thePhDSchoolDokIn’Holz fundedbytheAustrian Federal Ministry of Science, Research and Economy and the Austrian Association of Wood Industries is gratefully acknowledged. We also gratefully acknowledge the financial support of this work by the Austrian Science Fund (FWF) through the Erwin Schrödinger Fellowship J3748-N30. References 1. Anderheggen E., Knöpfel H., 1972, Finite element limit analysis using linear programming, International Journal of Solids and Structures, 8, 12, 1413-1431 2. Andersen E.D., Roos C., Terlaky T., 2003, On implementing a primal-dual interior-point method for conic quadratic optimization,Mathematical Programming, 95, 2, 249-277 3. Belytschko T., HodgeP.G., 1970, Plane stress limit analysis by finite elements, Journal of the Engineering Mechanics Division, 96, 6, 931-944 4. BotteroA.,NegreR., Pastor J., Turgeman S., 1980,Finite elementmethod and limit ana- lysis theory for soilmechanics problems,ComputerMethods inAppliedMechanics andEngineering, 22, 1, 131-149 5. Chen J., Yin J.-H., Lee C.F., 2003, Upper bound limit analysis of slope stability using rigid finite elements and nonlinear programming,Canadian Geotechnical Journal, 40, 742-752 6. Ciria H., Peraire J., 2004, Computation of upper and lower bounds in limit analysis using second-order cone programming andmesh adaptivity, 9th ASCE Specialty Conference on Probabi- listic Mechanics and Structural Reliability 7. Ciria H., Peraire J., Bonet J., 2008,Mesh adaptive computation of upper and lower bounds in limit analysis, International Journal for Numerical Methods in Engineering, 75, 8, 899-944 8. Dorn M., 2012, Investigations on the Serviceability Limit State of Dowel-Type Timber Connec- tions, PhDThesis, Vienna University of Technology 9. Drucker D.C., Greenberg H.J., Prager W., 1951, The safety factor of an elastic-plastic body in plane strain, Journal of Applied Mechanics, 18, 371-378 432 M. Li et al. 10. Drucker D.C., Prager W., Greenberg H.J., 1952, Extended limit design theorems for con- tinuous media,Quarterly of Applied Mathematics, 9, 381-389 11. Füssl J., Lackner R., Eberhardsteiner J., Mang H.A., 2008, Failure modes and effective strength of two-phasematerials determined bymeans of numerical limit analysis,ActaMechanica, 1-4, 195, 185-202 12. Füssl J., Li M., Lukacevic M., Eberhardsteiner J., Martin C.M., 2017, Comparison of unit cell-based computationalmethods for predicting the strength ofwood,Engineering Structures, 141, 427-443 13. Hawksbee S., Smith C., Gilbert M., 2013, Application of discontinuity layout optimization to three-dimensionalplasticityproblems,Proceedings of theRoyal Society of LondonA:Mathematical, Physical and Engineering Sciences, 469, 2155 14. HochreinerG., Füssl J., Eberhardsteiner J., 2013,Cross-laminated timber plates subjected to concentrated loading – experimental identification of failure mechanisms, Strain, 50, 1, 68-81 15. Hochreiner G., Füssl J., Serrano E., Eberhardsteiner J., 2014, Influence of wooden board strength class on the performance of cross-laminated timber plates investigated bymeans of full-field deformationmeasurements, Strain, 50, 2, 161-173 16. Krabbenhøft K., Lyamin A.V., Hjiaj M., Sloan S.W., 2005, A new discontinuous upper bound limit analysis formulation, International Journal for Numerical Methods in Engineering, 63, 7, 1069-1088 17. Le C.V., Askes H., Gilbert M., 2019, Adaptive element-free Galerkin method applied to the limit analysis of plates, Computer Methods in Applied Mechanics and Engineering, 199, 37, 2487-2496 18. Li M., Füssl J., Lukacevic M., Eberhardsteiner J., Martin C.M., 2018, Strength predic- tions of clear wood at multiple scales using numerical limit analysis approaches, Computers and Structures, 196, 200-216 19. Liu F., Zhao J., 2013, Upper bound limit analysis using radial point interpolation meshless method and nonlinear programming, International Journal of Mechanical Sciences, 70, 26-38 20. LukacevicM., Füssl J., 2016,Application of amultisurface discrete crackmodel for clearwood taking into account the inherent microstructural characteristics of wood cells, Holzforschung, 70, 9, 845-853 21. Lukacevic M., Füssl J., Lampert R., 2014, Failure mechanisms of clear wood identified at wood cell level by an approach based on the extended finite elementmethod,Engineering Fracture Mechanics, 144, 158-175 22. Lyamin A.V., Sloan S.W., 2002a, Lower bound limit analysis using non-linear programming, International Journal for Numerical Methods in Engineering, 55, 5, 573-611 23. Lyamin A.V., Sloan S.W., 2002b, Upper bound limit analysis using linear finite elements and non-linear programming, International Journal for Numerical and Analytical Methods in Geome- chanics, 26, 2, 181-216 24. Lysmer J., 1970, Limit analysis of plane problems in soilmechanics,Journal of the Soil Mechanics and Foundations Division, 96, 4, 1311-1334 25. MaierG.,Zavelani-RossiA.,BenedettiD., 1972,Afinite elementapproachtooptimaldesign of plastic structures in plane stress, International Journal for Numerical Methods in Engineering, 4, 4, 455-473 26. Makrodimopoulos A., Martin C.M., 2006, Lower bound limit analysis of cohesive-frictional materials using second-order cone programming, International Journal for Numerical Methods in Engineering, 66, 4, 604-634 27. Makrodimopoulos A., Martin C.M., 2007, Upper bound limit analysis using simplex strain elements and second-order cone programming, International Journal for Numerical and Analytical Methods in Geomechanics, 31, 6, 835-865 A numerical upper bound formulation with sensibly-arranged velocity... 433 28. Makrodimopoulos A., Martin C.M., 2008, Upper bound limit analysis using discontinuous quadratic displacement fields, Communications in Numerical Methods in Engineering, 24, 11, 911-927 29. Martin C., MakrodimopoulosA., 2008, Finite-element limit analysis ofMohr-Coulombmate- rials in 3D using semidefinite programming, Journal of Engineering Mechanics, 134, 4, 339-347 30. Milani G., Lourenço P.B., 2009, A discontinuous quasi-upper bound limit analysis approach with sequential linearprogrammingmeshadaptation, International Journal ofMechanical Sciences, 51, 89-104 31. MOSEK ApS, 2014, The MOSEK optimization tools version 7.0, User’s Manual and Reference, http://www.mosek.com 32. Prandtl L., 1920, Über die Härte plastischerKörper,Nachrichten von der Gesellschaft derWis- senschaften zu Göttingen, Mathematisch-Physikalische Klasse, 12, 74-85 33. Sloan S.W., Kleeman P.W., 1995, Upper bound limit analysis using discontinuous velocity fields,Computer Methods in Applied Mechanics and Engineering, 127, 14, 293-314 34. Smith C., Gilbert M., 2007,Application of discontinuity layout optimization to plane plasticity problems,Proceedings of the Royal Society A, 463, 2461-2484 35. WuJ.-Y., Cervera M., 2014,On the equivalence between traction- and stress-based approaches for the modeling of localized failure in solids, Journal of the Mechanics and Physics of Solids, 82, 137-163 36. YuS., ZhangX., SloanS.W., 2016,A3Dupper bound limit analysis using radial point interpo- lationmeshless method and second-order cone programming, International Journal for Numerical Methods in Engineering, 108, 13, 1686-1704 37. ZouainN., Borges L., Silveira J.L., 2002,Analgorithm for shakedownanalysiswith nonlinear yield functions,Computer Methods in Applied Mechanics and Engineering, 191, 23, 2463-2481 Manuscript received January 11, 2017; accepted for print February 5, 2018