A polytree-based adaptive approach to limit analysis of cracked structures

Accepted Manuscript A polytree-based adaptive approach to limit analysis of cracked structures H. Nguyen-Xuan, S. Nguyen-Hoang, T. Rabczuk, K. Hackl PII: DOI: Reference:

S0045-7825(16)30785-X http://dx.doi.org/10.1016/j.cma.2016.09.016 CMA 11128

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 20 July 2016 Accepted date: 15 September 2016 Please cite this article as: H. Nguyen-Xuan, S. Nguyen-Hoang, T. Rabczuk, K. Hackl, A polytree-based adaptive approach to limit analysis of cracked structures, Comput. Methods Appl. Mech. Engrg. (2016), http://dx.doi.org/10.1016/j.cma.2016.09.016 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A polytree-based adaptive approach to limit analysis of cracked structures H. Nguyen-Xuan1,2*, S. Nguyen-Hoang2, T. Rabczuk3, K. Hackl1 1

Institute of Mechanics, Ruhr-University of Bochum, 44801 Bochum, Germany

2

Center for Interdisciplinary Research in Technology (CIRTech), Hutech University, Ho Chi Minh City, Vietnam

3

Institute of Structural Mechanics, Bauhaus of University, Marienstrae 15, 99423 Weimar, Germany

Abstract We in this paper present a novel adaptive finite element scheme for limit analysis of cracked structures. The key idea is to develop a general refinement algorithm based on a socalled polytree mesh structure. The method is well suited for arbitrary polygonal elements and, furthermore, traditional triangular and quadrilateral ones, which are considered as special cases. Also, polytree meshes are conforming and can be regarded as a generalization of quadtree meshes. For the aim of this paper, we restrict our main interest in plane-strain limit analysis to von Mises-type materials, yet its extension to a wide class of other solid mechanics problems and materials is completely possible. To avoid volumetric locking, we propose an approximate velocity field enriched with bubble functions using Wachspress coordinates on a primal-mesh and design carefully strain rates on a dual-mesh level. An adaptive mesh refinement process is guided by an L2 -norm-based indicator of strain rates. Through numerical validations, we show that the present method reaches high accuracy with low computational cost. This allows us to perform large-scale limit analysis problems favorably.

Keywords: Fracture; plasticity; incompressibility; polytree; adaptivity; limit analysis.

*

Corresponding author. Email address: [email protected] (H. Nguyen-Xuan).

1

1 Introduction The adaptive finite element approach to limit analysis was roughly discussed in the work of Christiansen [1]. Mesh-adaptive refinement strategies were thereafter investigated in [2-4]. The goal of any adaptive refinement scheme is to locally generate the mesh in zones of the problem domain where local singularities occurs due to high plastic deformation. To guide such adaptive meshes, a mesh refinement indicator needs to be devised. Basically, the errorbased indicator is needed in order to guide the adaptive mesh refinement. Nevertheless, since addressing a priori error estimate to the limit analysis problem is not really simple [2], the most suitable choice is to use a posteriori error scheme. Subsequently, Borges et al. [2, 3] proposed the directional error estimator (or the directional indicator) for an anisotropic adaptive strategy behind a mixed limit analysis formulation. Their method took the plastic multiplier field characterized for slip bands or localized plastic deformations as the control variable. As a result, localizing regions with the presence of high plastic deformation are adaptively refined and the improvement of the computational collapse loads is really desirable. Lyamin et al. [5] showed that extending the directional error estimator can be applied to lower bound limit analysis of soil mechanics. It was shown that a mesh-adaptive scheme, where the Lagrange multipliers are taken as the control variables, yields more accurate results for the collapse loads, especially with the presence of ‘fan’ zones due to stress singularities. In the other approach, Christiansen and Pedersen [4] used the 2-norm of strains or stresses (taken as the control variable) as the error indicator to measure locally plastic collapse zones. However, numerical observation showed that the adaptive mesh refinement is propagated far away from the undergoing plastic regions or is performed in large areas of the rigid parts of the solid. Other error indicators predominantly related to the dissipation gap or the difference between the lower and upper bounds were presented by Ciria et al. [6,7] and Munoz et al. [8]. Martin and co-workers [9-12] utilized dissipation-based adaptivity to obtain rigorous lower and upper bounds of the collapse multiplier. Further investigations [13] addressed this bounds gap to geotechnical stability analysis. Other application of the plastic dissipation-based indicator for plate structures was reported in [14]. Very recently, the plastic dissipation indicator was also involved in adaptive triangular finite element formulation enriched with bubble nodes for plane-strain limit analysis [15-17]. In addition, the error indicator based on yield stresses [18] was proposed for adaptive quasilower bound limit analysis of solids. In this paper, we exploit an L2 -norm-based indicator of strain rates to guide the adaptive mesh refinement process in high plastic deformation zones. 2

Note that such an indicator based on the L2 -norm of the deviatoric strain rate was previously used for solving shear bands in plane strain compression of a thermoviscoplastic solid [19]. Another key feature of any adaptive finite element technology is the automatically refinement mesh generator. In the context of limit analysis, the aforementioned approaches rely on triangular meshes. According to the authors’ best knowledge, there are no publications concerning adaptive polygonal meshes for limit analysis problems. From a computational point of view, n-node polygonal elements perform, broadly speaking, better than three-node triangular elements. Nevertheless, the design of a unified automatic mesh generator which is valid for triangles, quadrilaterals and polygons is really non-trivial. Even in cases using quadtree data structure [20], the generation of hanging nodes after each refinement step causes incompatibility in standard finite element approximations. As a consequence, this shortcoming required a special treatment with the enhanced complexity in both algorithms and interpolations [21-25]. Tabarraei and Sukumar [26] used natural neighbor (Laplace) basis functions over quadtree meshes to facilitate the C 0 -compatible requirement of approximations along edges being the presence of hanging nodes. The crucial idea is to delineate the interpolant over polygonal reference (or local) elements and then employ an affine map to establish conforming approximations on quadtree meshes. As a result, this approach enables an arbitrary number of hanging nodes per edge. It offers high flexibility in the automatically refinement mesh generator and adaptive numerical computations based on the quadtree data structure. We observe that the quadtree data structure cannot be applied directly to arbitrary polygon meshes, while various applications based on polygons are of increasing interest. In our work, it is desirable to investigate a polytree data structure and exploit the efficiency of the interpolant over polygonal reference elements, and also use an affine map to establish conforming approximations on polytree meshes for limit analysis problems. It was shown that polygonal elements possess greater flexibility in mesh design for complicated problems [27-29]. Nowadays, polygonal finite elements have been extensively applied to mechanics problems such as nonlinear constitutive modeling

of

polycrystalline

materials

[30,31,32],

topology

optimization

[33-35],

incompressible materials [36], fracture modeling [37,38,39] and so on. Limit analysis plays major roles in the context of elastic-plastic fracture and fracture-safe design.

Since

lower-order

displacement

finite

elements

perform poorly

in

the

incompressibility condition [40], several triangular or quadrilateral finite element

3

improvements have been revised in the literature. For example, by no means exhaustive, a mixed formulation based on an enhanced strain field was proposed in [41]. Several finite element formulations enriched with discontinuity velocity fields were reported in [42,43,44]. Linear strain finite elements (or quadratic velocity fields) where the strains vary as a simplex were investigated in [45]. Very recently, linear triangular finite element formulations enriched with bubble nodes for plane-strain limit analysis were formulated in [15-17]. However, a study on polygonal finite elements is also a necessary choice together with other existing elements in the context of limit analysis. This is our goal in this article. We in this paper propose an adaptive polygonal finite element method for limit analysis of cracked structures. We construct a general refinement algorithm based on a concept of a so-called polytree mesh structure. It is evident that polytree meshes are identical to quadtree meshes when meshes are triangles or quadrilaterals and the recursive partition of each element into four child elements is performed. Two meshing levels are used for adaptive meshing refinements and computations: a primal-mesh level is placed with the kinematic formulation relying on barycentric basis functions over polygonal elements, while a dualmesh level is used to compute numerical integrations and to enforce the flow constraint. We numerically show that the present formulation is well suited for the volumetric locking phenomenon. The discrete limit analysis problem can be further formulated in the compact form of the well-known second-order cone programming (SOCP) with the aim of efficiently utilizing the interior-point method. More importantly, an adaptive mesh procedure guided by the L2 -norm-based indicator of strain rates is performed over the primal-level mesh background to enhance the desired precision of the underlying solution. The method can be straightforwardly implemented into polygonal finite element packages. Finally, several benchmark problems are given to validate the robustness of the present method. The article is organized as follows: Section 2 briefly describes the kinematic theorem. Section 3 devotes a polygonal finite element limit analysis using Wachspress coordinates with an edge-based integration scheme. A polytree-based adaptive mesh refinement is given in Section 4. Numerical examples are presented in Section 5. Finally, we close this paper with some concluding remarks. 2 A brief on the kinematic theorem Consider a cracked structure made of rigid - perfectly plastic materials occupying the open domain Ω ⊂ R 2 bounded by continuous and discontinuous boundary parts such that 4

Γ = Γu ∪ Γt ∪ Γc , Γu ∩ Γt ∩ Γc = ∅, which Γc contains discontinuous boundary part. We

define a convex set S consisting of an admissible stress field σ such that S = {σ ∈ Σ ψ (σ ) ≤ κ}

(1)

where Σ contains symmetric stress tensors, ψ (σ ) is the yield function, and κ is the positive constant. Regarding plane strain problems, the von Mises yield function [55] reads

(σ 11 − σ 22 )

ψ (σ ) =

4

2

+ τ 122

(2)

Following the associative plasticity rule, plastic strain rates are amenable to the yield surface given by

ε= μ

∂ψ (σ ) ∂σ

(3)

where μ is the plastic multiplier. Also, we provide the L2 (Ω) -norm for strain rates which will be mentioned later as an alternative indicator in guiding polytree mesh refinements

η= ε

L2 ( Ω )

=

(∫ ε : ε)

1/ 2

Ω

(4)

Assumed that the structure is subjected to variable loads such as body forces λ f and external tractions λ g on the Neumann boundary Γt , a prescribed velocity vector u is given on the Dirichlet boundary Γu . The key of limit analysis is to seek a limit (or safety) load multiplier λ∗ at which the structure will collapse. According to Koiter’s kinematic (upper bound) theorem [46], the limit analysis of structures aims at finding a load factor that results in minimizing plastic dissipation (see Appendix for a detailed expression of von Mises plasticity) for any kinematically admissible velocity field, i.e,

λ * = min D(u)

(5)

u∈W

where D (u ) = ∫ κ Ω

{

( ε11 − ε 22 )

2

+ γ 122 dΩ

(6)

}

and W = u ∈ V (u) = ∫ f ⋅ udΩ + ∫ g ⋅ udΓ = 1 with a space of kinematically admissible Ω

Γt

velocity field V denoted by

5

V = {u ∈ ( H1 (Ω )) 2 , u = u on Γ u }

(7)

where κ = σ y / 3 , and σ y is the yield stress. It is known that plastic dissipation [40] is heavily constrained by the incompressibility condition, i.e, ∇ ⋅ u = 0 . From the computational point of view relying on lower-order finite element formulations, there remain challenging issues in the fulfillment of the incompressibility condition [40,42]. Contributing into the current study, we propose a polygonal finite element formulation with the enhancement of bubble functions for resolving this numerical difficulty. It is worth noting that our approach uses only a minimum number of velocity degrees of freedom and requires a minimum number of the incompressibility constraints. 3 On polygonal finite elements 3.1 A common scheme for lower-order polygonal finite elements enriched with bubble functions

It is assumed that an underlying domain Ω is discretized into a set T (at master level) of

ne polygons involving ns edges and nn nodes such that Ω ≈ Ω h = ∑e =1 Ω e . A finite element ne

approximation space V h ⊂ V of kinematically admissible velocity field is defined as

{

V h = u h ∈ ( H1 (Ω )) 2 , u h

Ωe

∈ P1 , u h

∂Ω e

}

∈C 0 , Ωe ∈ T

(8)

where P1 refers to so-called generalized barycentric coordinates associated with Ωe ∈ T . We state the discrete problem of (5) as: to find a collapse multiplier λ + :

λ + = min D h (u h )

(9)

⎧∇ ⋅ u h = 0 s.t. ⎨ h h ⎩ u ∈V

It is worth noting that, plane-strain problems, lower-order finite elements exhibit severe volumetric locking under plastic yield conditions [40,42-45]. Even for simple geometrical domains, they may not provide reasonable results when unstructured meshes are employed. Unfortunately, they alleviate partly volumetric locking, while the stability of the slip-line

6

pattern is not, in general, ensured. To remedy the shortcomings, we enhance a space V h with bubble functions which have been proved to be really necessary in the plane-strain limit analysis [15]. Ideally, bubble functions whose values are zeros on the element boundaries are used for any interior point of finite elements. For a polygonal element, such interior functions are formulated by the triangular simplex [58]. These are similar to bubble functions in the standard finite element literature. The supplementing space B h with a space of bubble functions is defined as A (x) ⎪⎧ ⎪⎫ , ∀Ωe ∈ T ⎬ B h = ⎨φbe = A ( x) + B ( x) ⎩⎪ ⎭⎪

(10)

where A (x ) vanishes on the boundary line, and B (x ) is minimum at the centroid of element Ω e and is smooth over the rest of the domain [58]. The construction of bubble shape

functions for polygonal elements with interior nodes will be presented in Section 3.2.2. The finite element space Vbh which is enriched with bubble functions is written as Vbh = V h ⊕ B h

(11)

Unlike the mixed displacement-pressure finite element formulation, e.g, the MINI element [59], our approach employs only a displacement-based formulation. As shown in [15], the inf-sup condition can be ensured. The displacement velocity u h ∈ Vbh on Ω e ∈ T is explicitly written as

uh (x) = ∑i =ne1 (φie (x)I 2 )die + (φbe (x)I 2 )dbe n

(12)

where nne denotes the number of vertices of polygonal element, I 2 is the unit matrix of 2nd rank, die is velocity degrees of freedom of u h (x) with respect to the ith vertex, φie is the shape function at the ith vertex, dbe is the unknown displacement associated with an extra centroid node ubh (x), and φbe indicates the bubble shape function at the centroid node of the element. 3.2 Shape functions on convex polygonal elements There have been numerous approaches to the construction of shape functions in polyhedral elements [60-71]. A comparison of different shape functions was depicted in [72,73]. In this 7

study, we focus on using Wachspress coordinates [60,62,70,71,73] with iso-parametric mapping on a reference element [74] to construct conforming approximations on polytree meshing structures. Wachspress coordinates have many nice properties [73]: •

Smooth C ∞ over polygonal domain and C 0 approximations along the boundary.

•

The partition of unity:

•

Kronecker-delta properties: φi (x j ) = δ ij

•

Non-negative: φi (x) ≥ 0

•

Linear precision property:

•

Wachspress coordinates are unchanged with any affine transformation [73]:

∑

n

φ ( x) = 1

i =1 i

∑

n

φ ( x) x i = x

i =1 i

φi (α x, α x1 , α x 2 ,… , α x n ) = φiW (x, x1 , x 2 ,… , x n ) ∀α ∈ R 2 . 3.2.1 Wachspress shape functions

Consider a polygon Ωe ∈ T with vertices x1 , x 2 ,… , x nne , nne ≥ 3 , in counter-clockwise ordering. For any v ∈ Ω e , Wachspress shape functions [60,62,70] can be defined by

φie =

∑

ϕi nne j =1

ϕj

=

A ( x i −1 , x i , x i +1 ) ϕi , with ϕ i = A ( v , x i −1 , x i ) A ( v , x i , x i +1 ) ψ

(13)

where A(x a , x b , x c ) denotes the signed area of the triangle [x a , x b , x c ] , see Fig. 1(a). An alternative approach to Wachspress coordinates is to use the perpendicular distances of v to the edges of Ω e [71]. Let hi (x ) denote the perpendicular distance of v to the edge e i , as shown in Fig. 1(b), and define p i ( x) = n i / hi (x) , where n i is the outward unit normal vector to the edge e i = [x i , x i +1 ] , with vertices indexed cyclically x n+1 = x1 . Then the Wachspress shape functions and their gradients are expressed as

ϕi , with ϕ i = det(p i −1 , p i ) ψ

(14)

∇ φie = φie [ϑi − ∑ jne=1 φ ej ϑ j ] where ϑ i = p i −1 + p i

(15)

φie =

ϕi = ψ

∑

ϕi nne j =1

ϕj

=

and n

Fig. 1(c) and Fig. 1(d) show 3D view and contour line of a shape function for an arbitrary node of a hexagonal element, respectively.

8

3.2.2 Bubble Wachspress shape functions Malsch and Dasgupta [58] constructed bubble shape functions over arbitrary convex polygons with interior nodes by using square of a distance function and a function which vanishes on the boundary. As presented in Section 3.1, a function A (x ) can be derived from the triangular simplex. As stated in Section 3.2.1, we can use the ratio of products of the Wachpress coordinate in Eq. (13) or Eq. (14) to establish A (x ) as A (x ) =

1

ψ

1 = ~ ψ

(16)

In Eq. (10), B (x) can be defined in a simple way as B (x ) = x − x b

2

(17)

Ωe

where x b is the coordinates of the centroid of the element, Ω e . As a result, a bubble shape function φbe ∈ B h can be expressed as

φbe =

1 1 + x − xb

2 Ω

e

(18)

ψ

and its gradient is ∇ φbe = −

x − xb

2 Ωe

∇ψ + ∇ x − x b

( x − xb

2 Ωe

ψ + 1) 2

2 Ωe

ψ (19)

Fig. 2 shows the features of the bubble shape function φbe of triangle, quadrilateral and hexagonal elements. 3.3 The area-averaged edge-based strain rates projection As an alternative integration technique, we exploit an area-averaged edge-based projection technique, which follows the concept of edge-based smoothed finite element methods (ES-FEM) [56]. The present method is numerically shown to be insensitive to volumetric locking in incompressible limit. Regarding how to efficiently utilize the areaaveraged edge-based integration, a dual mesh based on edge-based mesh background T needs to be defined by connecting all vertices and centroids of elements of T . The domain Ω can be also partitioned into ns edge-based smoothing domains Ω ( k ) relied on edges of primal

9

elements such that Ω = ∪ks=1 Ω ( k ) and Ω ( i ) ∩ Ω ( j ) = ∅ for i ≠ j , which Ω (k ) is created by two n

sub-triangles having the common edge k as depicted in Fig. 3. Subsequently, we use a smoothing or constant (average) projector c h ( k ) to define an assumed strain rate on a dual mesh following the edge-based smoothing technique [56]: 1

ε ( ) = ∇ s u h = c h ( k )∇ s u h = k

where A(k ) = ∫

Ω( k )

A(

∫

k)

Ω

(k )

∇ s u h dΩ

(20)

dΩ denotes area of the smoothing domain Ω (k ) and ∇ s is the matrix of

differential operators. The volume expansion rate is ∇ ⋅ u h = c h ( k )∇ ⋅ u h =

1 A(

k)

∫

Ω

(k )

∇ ⋅ u h dΩ

(21)

which, under the incompressibility, leads to ∀Ω( k ) ∈ T , ∇ ⋅ u h = 0

(22)

Similar to the theoretical proof provided in [15], volumetric locking can be resolved through a careful enforcement of the incompressibility condition on the dual mesh via an alternative integration technique based on area-averaged edge-based projection.

3.4 Limit load approximation The discrete problem is to determine an approximated collapse multiplier λ +

λ + = min D h ( u h ) (23)

⎪⎧∇ ⋅ u h = 0 such that ⎨ h h ⎪⎩ u ∈ Vb where

D h (u h ) = ∫ κ Ω

(ε

− ε 22 ) + γ 122 dΩ = ∑ ks=1 κ A( k ) 2

11

n

(ε

− ε 22( k ) ) + γ 12( k )

(k ) 11

2

(24)

As observed from (24), the integration is numerically computed on the smoothing domains {Ω ( k ) }k =1, 2,…,ns and this is equivalent to one Gauss point used for each smoothing domain. Hence, only one incompressibility constraint is imposed on each entire smoothing 10

domain. Because the strain rates defined on the dual mesh via the edge-based strain smoothing are constant over the entire domain Ω (k ) , the flow condition is well-posed at any one point over Ω (k ) , and likewise everywhere in the whole problem domain. Hence, the present approach can be considered as a simply alternative form of constant strain elements, and, importantly, it attenuates volumetric locking. Last but not least, the limit analysis problem (23) forms a non-linear optimization problem with equality constraints. It turns out that this optimization model can be reformulated as

λ + = min ∑ k =1 κ A( k ) ρ( k ) ns

⎧c h ( k ) ∇ ⋅ u h = 0, k = 1,… , ns ⎪ s.t. ⎨u h = u on Γ u ⎪ h ⎩Wext (u ) = 1

(25)

⎡ε ( k ) − ε ( k ) ⎤ where ρ( k ) = ⎢ 11 ( k ) 22 ⎥ denotes auxiliary variables. ⎣ γ 12 ⎦

By introducing auxiliary variables t1 , t2 ,…, tns , the optimization problem (25) becomes a standard form of the following second-order cone programming (SOCP)

λ+ = min ∑k =1 κA ( k ) t k ns

⎧t k ≥ ρ ( k ) ⎪ ⎪c ( k )∇ ⋅ u h = 0, k = 1, 2,… , ns s.t. ⎨ h ⎪u h = u 0 on Γu ⎪ ⎩ (u) = 1

(26)

As a result, the conic interior-point optimizer of the academic MOSEK package [75] is an awesome choice for solving this problem. 4 A polytree-based approach to adaptive polygonal finite elements In plastic collapse analysis, the quality of the solution (the resolution and clearness of the yield lines) can be improved by using a finer mesh. However, this greatly increases the computational cost, which may grow rapidly and become very expensive in large problems. Since there are usually a large percentage of computational domain and plastic zones in the

11

plastic collapse solution, it is unnecessary and uneconomical to place a uniform-dense set of finite elements all over the computational domain. Thus, a more efficient strategy is to start from a relatively coarse analysis mesh and then adaptively refine it only within the certain regions as necessary. Therefore, a novel adaptive polytree approach is presented in this section to solve the plastic collapse problems. 4.1 L2 - norm-based indicator It is known that the presence of the localized plastic deformations leads to the slow convergence of finite element limit analysis solutions using uniform meshes. Ideally, the mesh is refined adaptively in the regions characterized by high strain rates. It means that the mesh refinements need to be done in these regions so that the approximate strain rate reaches a desired accuracy with the minimum computational cost. With that observation, we use an indicator based on the L2 - norm of plastic strain rates over each element regarded as an alternative indicator in adaptive refinement scheme and computed as

ηe =

1 nede

∑

e ned

η (k )

k =1

(27)

where nede is the number of edges of an element and η ( k ) is the L2 - norm of plastic strain rates on edge k of the element e and is computed by following Eq. (4) as

η (k ) = ε (k )

L2 ( Ω( k ) )

(28)

4.2 L2 - norm-based indicator and refinement strategy The corresponding L2 -norm defined in (28) is now utilized to guide the mesh refinements. The global indicator, η , is contributed from the local one, η e , for all the elements as follows

η = ∑ e=1η e ne

(29)

Adopting the Dorfler criterion [76], one uses the local refinement indicator to filter a minimal set M ⊆ T of the elements Ωe ∈ T for refinement such that

∑η

Ω ∈M e

e

≥ θ η , for some θ ∈ (0,1)

(30)

12

Subsequently, a new mesh T' is generated from T by refining (at least) the marked elements Ωe ∈ M by using polytree meshes. 4.3 Polytree meshes 4.3.1 Polytree data structure In the finite element method, the quadtree decomposition [77,78] plays an important key in the adaptive mesh generation of “triangular” and “quadrilateral" meshes. However, adaptation of quadtree meshes produce hanging nodes between two adjacent elements of different sizes after each refinement. In the literature, several approaches have been devised to overcome this incompatibility. Among them, a simple and effective method based on conforming polygon elements was presented by Tabarraei and Sukumar [79]. Herein, each quadtree element with the presence of hanging nodes can be considered as a polygon with nodes along the sides, called side-nodes. By employing C 0 approximations along edges of natural neighbor (Laplace) shape functions and iso-parametric map from polygonal reference elements, we can compute the polygonal element with n nodes and m side-nodes as the polygonal element with (n+m) nodes. Therefore, the restriction of the number of hanging nodes on each side is no need. Addressing advantages of polygonal meshes and adaptive technique of quadtree meshes, we introduce a novel spatial decomposition structure for polygonal mesh refinements, called a polytree. Unlike the quadtree decomposition with initial triangular or quadrilateral meshes, a general polytree decomposition follows an adaptation using initial polygonal meshes, which allows for the use of unstructured meshes for arbitrary domains. In the polytree structure, the domain is arbitrary with initial polygonal mesh as shown in Fig. 4(a). As shown in Fig. 4(b) and Fig. 4(c), a polygonal element (including triangular and quadrilateral elements) with n nodes (father) is decomposed into (n+1) new polygonal elements (children), and each of children is partitioned recursively until reaching a stopping criterion. A polytree mesh representative tree is presented in Fig. 4(d). It is clear that hanging nodes can be generated after each refinement if children and their adjacent elements are not the same levels of refinement. Here we adopt the technique of the conforming quadtree meshes [79] as presented above to treat incompatibility of polygonal elements with its hanging nodes (side-nodes). In this article, shape functions of the polygonal elements with side-nodes in polytree meshes are obtained from Wachspress coordinates and iso-parametric map J ξ = ∂x / ∂ξ from polygonal reference elements as shown in Fig. 5. Hence, there is no difference between side-nodes and vertices while using Wachspress shape functions with C 0 13

approximations along edges [60,62] and in fact all the nodes are considered as vertices. It is worth noting that the construction of shape functions and all computations of triangular or quadrilateral elements with hanging nodes referring to the quadtree structure are efficiently integrated into the common framework of the polytree structure. Therefore, the quadtree algorithm by Tabarraei and Sukumar [79] can be considered as a special case of the polytree one. 4.3.2 Subdivision of polygonal elements on polytree mesh structure The key of the polytree is that a father at the previous step, with an arbitrary shape and n nodes, is further subdivided into (n+1) children. There are several cases where polygonal meshes are described in matching and non-matching new meshes [81]. Besides, by connecting the centroid of all the triangular elements circumventing a particular node of father, we can construct children from a triangular mesh background. A sequence of steps required to construct children from an underlying triangular mesh is shown in Fig. 6. As shown in Fig. 7, this method remains its limitations in a general case with many iterations. Recently, a simple and efficient polygonal mesh generator, PolyMesher, based on the context of centroidal Voronoi tessellation (CVT) written in Matlab, was developed by Talischi et al. [82]. By using centroidal Voronoi tessellation algorithm, PolyMesher can construct high quality polygonal meshes of an arbitrary domain shape, see Fig. 7. Note that the quality of the polygonal mesh generated by using CVT depends on a set of initial random points called seeds, which are firstly inserted within the domain of polygonal element. Here, the (n + 1) seeds are inserted easily at centroid of father and mid-point of lines that connects the centroid of father to n vertices of father. After that (n + 1) children are constructed within domain of father by using CVTs as shown in Fig. 8. Note that new nodes of children are located at the arbitrary places of the edges of father where the non-matching meshes can be created as shown in Fig. 9(a) (left). We can enhance the compatibility and matching for polytree meshes, see Fig. 9(a) (right), by dividing the edges of father into four equal segments with triangular nodes as shown in Fig. 9(b), and then moving new nodes to the nearest triangular node. Fig. 10 shows a sequence of subdivisions of arbitrary polygonal elements using CVTs. 5 Numerical validations

In this section, we demonstrate the robustness of the present method through several benchmark problems. We abbreviate the polygonal elements as follows

14

ƒ

W-Poly: n-node polygonal element,

ƒ

bW-Poly: n-node polygonal element enriched with a bubble node.

ƒ

MINI: mixed displacement/pressure finite element enriched with cubic bubble functions [59].

Unlike the classical FEM, we use an alternative quadrature scheme relying on the areaaveraged edge-based projection for polygonal elements. It is noted that using the standard quadrature rule, polygonal finite elements are sensitive to volumetric locking and do not perform well for plastic collapse plane-strain solids. A flowchart for solving limit analysis problems using the polytree meshing structure is also summarized in Fig. 11. The program is performed on ASUS X550C (Intel Core I5, 1.8GHz CPU, 4G RAM). 5.1 Cook’s membrane problem: test for volumetric locking The aim of this example is to examine the performance of the present formulation in the incompressible limit. The well-known Cook’s membrane problem under plane strain conditions is considered. This example has been known as a standard classical benchmark for incompressible problems. Fig. 12 describes the geometry and boundary condition. The leftside boundary is fixed, and a vertical load is applied to the right-side boundary. The material properties are taken to be E = 1 and ν = 0.4999999 . Four discretizations using (n+1)-node polygonal elements with enriched nodes are shown in Fig. 13. The reference values [84] of the strain energy and the vertical displacement at the center tip section are 9.27769 and 18.5002, respectively. In Fig. 14, the numerical pressure field with respect to mesh #4 (950 nodes and 316 elements) using bW-Poly and others are compared. It is indicated that strong oscillations of the pressure distribution are found using the W-Poly element, while MINI and bW-Poly elements work well. MINI element exhibits small oscillation of the pressure distribution which is more pronounced near the right-side boundary, while the present element produces stable results. The relative errors in log-log scale of the vertical displacement and the strain energy are plotted in Fig. 15(a) and Fig. 15 (b), respectively. Herein, the average length h = AΩ / ne (where AΩ stands for the area of the entire problem domain) is used. It is indicated that the convergence rates of the bW-Poly element (r = 2.68 and 2.4) are much higher than the MINI element (r = 1.47 and 1.45). It is observed that the present method is volumetric locking-free and shows good performance for nearly incompressible solids.

15

5.2 Notched tensile specimen This is the benchmark problem occured in [40] and has become popular in many researches. Fig. 16 (a) shows a notched tensile specimen with the presence of two external thin symmetric cuts subjected to in-plane tensile stresses q. Due to the unknown exact solution, the best reference upper bound value of the collapse multiplier or limit load factor using the Richardson's extrapolation was found to be λ * = 0.9241, 1.1316, and 1.3833 with respect to a = 1/3, 1/2, and 2/3. An n-node polygonal meshes enriched with centroid nodes are plotted in Fig. 17 (b) for modelling the upper-right quarter due to its symmetry. Again, we compare the results between W-Poly and bW-Poly for the crack length a = 1/ 2 . Defining N var as the total number of variables of the optimization problem, one has

N var = ndof + 3ns . Fig. 17 shows the convergence of the collapse multipliers with respect to Nvar. As mentioned in subsection 3.1, the W-Poly element without the enrichment of bubble functions gives reasonable results in this mesh case. However, as shown in Fig. 18(a), this element does not properly reproduce the distribution of plastic dissipation and collapse mechanism. As shown in Fig. 17 and Fig. 18 (b), the bW-Poly element works well. Without loss of generality, we consider the computational efficiency of the bW-Poly element. Table 1 shows the computational cost evaluated through steps, N var , and optimizer times. Regarding the maximum number of variables ( N var = 69547 ), the computational time for solving the optimization problem based on the interior-point algorithm only takes around 102 seconds. This proves that the present approach is totally suitable for solving large-size sparse SOCP problems. Next, we analyze the specimen for various crack lengths using a polytree meshing structure presented in Section 4 . The mesh refinement process is initiated at a uniformly coarse mesh as shown in Fig. 19(a). Several mesh steps are then illustrated in Fig. 19(b-e). Fig. 20(b) exhibits a super-convergence rate when polytree mesh refinement is used. This proves the high efficiency of polytree refinement in comparison with the uniform one. Similarly, Fig. 21 shows the final polytree meshes for two crack lengths a = 1/ 3 and 2/3, respectively. The relative error of the present solutions versus number of variables N var for two crack lengths a = 1/ 3 and 2/3 are shown Fig. 20(a) and Fig. 20(c). As can be seen, the present method provides reliable solutions for uniform and polytree meshes. It is evident that the mesh is refined adaptively in collapse plastic zones characterized by high strain rates. As

16

shown in Table 2, the relative error percentage in the collapse multiplier demonstrates the high performance of the present method based on polytree meshes. For comparison, Table 3 shows that the present results are in good agreement with other published ones.

5.3 Central cracked plate subjected to tension Consider a plate subjected to tension as shown in Fig. 22(a). The plate has dimension b × H (with b = 2, H = 4 ) and a center crack length a. The analytical collapse multiplier is

given by [50]

λ* =

2 (1 − a / b) 3

(31)

Due to its symmetry, the upper-right quarter of the plate is modeled. We verify for three cases x = a / b . Final polytree meshes for x = 1/ 4 , x = 1/ 2 and x = 3 / 4 are displayed in Fig. 23, respectively. The collapse multipliers are plotted in Fig. 24. These results are also listed in Table 4. Present solutions converge well with the analytical ones. It is seen again that the present method using the polytree meshes retains a low number of variables while producing numerical solutions which match very well with the analytical values. For example, in the x = 1/ 2 case, using a polytree final mesh yields λ + = 0.5803 (with the error of 0.5%) with a low number of 15932 variables, while the collapse multiplier value λ + = 0.5965 (with the error of 3.32%) requires a uniformly fine mesh of 22507 variables.

5.4 Single-edge cracked plate subjected to tension Consider a plate of length H width b and a single edge crack of length a subjected to tension. The system and its symmetric model are depicted in Fig. 25. The analytical solution with respect to x = a / b was reported in [48] as

λ=

{

}

3.404 ⎡ 2 2 1/2 2.06 − x ) + 0.5876 (1 − x ) ⎤ + ( 2.06 − x ) , ∀x ≥ 0.545 ( ⎦ 3 ⎣

(32)

and in case x < 0.545 , the solution is bounded by [50]

(

)

(

2 2 2 1 − x − 1.232 x 2 + x3 ≤ λ ≤ 1 − x − 1.232 x 2 + x3 + 22 x3 ( 0.545 − x ) 3 3

)

(33)

17

For x = 0.5 , the actual value of the collapse multiplier falls into a small range between 0.366 and 0.373. Final polytree meshes for the edge-cracked tensile specimen with x = 1/ 4, x = 1/ 2 and x = 3 / 4 are displayed in Fig. 26. As expected, a high density of polygon

elements in collapse plastic regions is shown. Fig. 27 plots the collapse multiplier versus Nvar for various crack lengths. In addition, Table 5 summarizes present solutions and their corresponding errors to analytical values. A very good convergence is confirmed. Plastic dissipation profile plotted in Fig. 28 shows a good agreement with the one reported in [53]. 5.5 Cracked cylinders under internal pressure Finally, we deal with a cylinder with longitudinal crack a subjected to internal pressure, cf. Fig. 29. The cylinder has inner radius Ri and outer radius R0 . Analytical solutions and numerical solutions are available in the literature. For comparison, several following investigations are summarized as below •

Approximate bound of Chell [49]:

h−a 2 σy p Ri + a h−a 3 λa = L = = p0 p0 (Ri + a ) ln R0 Ri

(34)

where h = R0 − Ri and p0 stands for the limit pressure of the un-cracked cylinder given by p0 =

•

2 R σ y ln 0 Ri 3

Lower bound solution of Miller [50]: ⎛ R ⎞ ⎛ R ⎞ 2 σ y ln ⎜ 0 ⎟ ln ⎜ 0 ⎟ 3 p ⎝ Ri + a ⎠ = ⎝ Ri + a ⎠ λ− = L = R p0 p0 ln 0 Ri

•

(35)

(36)

Approximate analytic solution of Yan and Nguyen-Dang [52]:

18

R 2 − ( Ri + a ) 1 σy 0 2 R02 − ( Ri + a ) R0 Ri pL 3 a = = λ = R p0 p0 2 R0 Ri ln 0 Ri 2

(37)

In addition, finite element solutions are also used for comparison: •

Yan [52] a a2 λ = 1 − 0.7716 − 0.2267 2 h h

(38)

a

•

Kim et al. [83]

λa =

pL = p0

2 2h σy R0 + Ri 3

⎛ a a2 a3 ⎞ 1 0.356 1.6882 1.0442 − − + ⎜ ⎟ h h2 h3 ⎠ ⎝ p0

(39)

⎛ a a2 a3 ⎞ 2 h ⎜ 1 − 0.356 − 1.6882 2 + 1.0442 3 ⎟ h h h ⎠ ⎝ = R ( R0 + Ri ) ln 0 Ri

Our computations are performed for three cases of cracked length and final polytree meshes are provided in Fig. 30.

It is, once again, found that the mesh refinement

concentrates primarily in plastic zones. The profiles of the plastic dissipation match well with those of singular 8-node finite element reported in [52]. Finally, Fig. 32 closes numerical tests with the high reliability of the present results. 6 Conclusions

Plastic collapse analysis of cracked plane-strain solids has been investigated in this paper. The numerical results obtained showed the high efficiency of the proposed method. Our approach was to establish a polygonal finite element formulation behind Wachspress shape functions with the presence of bubble basis functions defined on the polytree primal-mesh structure. The plastic dissipation and incompressibility constraints were performed on the polytree dual-meshes through the edge-based mesh repartitioning scheme. The area-averaged edge-based projection technique was used with the aim to redistribute the plastic strain rates in order to guarantee a minimum number of incompressibility constraints, and ensure that the flow rule constraint held everywhere in the entire problem domain. In addition, the limit

19

analysis problem was solved efficiently by interior-point solvers, which are robust for largescale optimization problems. More importantly, we showed that the automatic local refinement procedure guided by the L2 -norm-based indicator of strain rates in combination with the polytree mesh structure is really useful for plastic collapse problems. Finally, the present method can be straightforwardly extended into various mechanics problems without more requirements of implementation and computational cost. Acknowledgments

The first author would like to thank the Alexander von Humboldt Foundation for granting the Georg Forster Research Award. This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2014.24. These supports are gratefully acknowledged. Appendix

We express the von Mises criterion for perfectly rigid-plastic materials as

ψ = J2 = κ

A1

1 1 J 1 = σ ii , sij = σ ij − J 1δ ij , J 2 = s ij s ij 2 3

A2

where

and δ ij is the Kronecker delta. The relation between strain rates and the yield function is expressed through the associated yield rule:

εij = μ

s ∂ψ = μ ij ∂σij 2 J2

A3

Hence, the plastic multiplier is computed by

μ = ( 2ε ij ε ij )

1/ 2

A4

Otherwise, the power of the plastic dissipation is written as

 = κμ D = σij εij = μψ

A5

D = κ (2ε ij ε ij )1/2

A6

Therefore, one has

For plane strain case ( ε 33 , γ 13 , γ 23 = 0 ), s13 , s23 are zero and

20

σ 33 =

σ 11 + σ 22

A7

2

The invariants reduce to:

J2 =

1 sij sij 2

2 2 2 ⎛ σ 11 − σ 22 ⎞ ⎛ σ 11 − σ 22 ⎞ ⎤ 2 1⎡ = ⎢(σ 11 − σ 22 ) + ⎜ ⎟ +⎜ ⎟ ⎥ + τ 12 6⎢ 2 2 ⎝ ⎠ ⎝ ⎠ ⎥⎦ ⎣ 2 ⎡1 ⎤ = ⎢ (σ 11 − σ 22 ) + τ 122 ⎥ ⎣4 ⎦

A8

The yield function becomes 2 ⎡1 ⎤ ψ = ⎢ (σ11 − σ 22 ) + τ122 ⎥ ⎢⎣ 4 ⎥⎦

1/ 2

=κ

A9

Substituting A9 into A3 and after chains of proof, we obtain

μ=

( ε 11 − ε 22 )

2

+ γ 122

A10

+ γ 122

A11

and the power of plastic dissipation becomes D =κ

(ε11 − ε 22 )

2

In addition, one obtains the incompressibility condition ∇ ⋅ u h = ε 11 + ε 22 = 0

A12

21

References 1. E. Christiansen, K.D. Andersen, Computation of collapse states with von Mises type yield condition, International Journal for Numerical Methods in Engineering 46 (1999) 1185–1202. 2. L.A. Borges, R.A. Feijoo, N. Zouain, A directional error estimator for adaptive limit analysis, Mechanics Research Communications 26 (5) (1999) 555 – 563. 3. L.A. Borges, N. Zouain, C. Costa, R. Feijoo, An adaptive approach to limit analysis, International Journal of Solids and Structures 38 (2001) 1707–1720. 4. E. Christiansen, O.S. Pedersen, Adaptive mesh refinement in limit analysis, International Journal for Numerical Methods in Engineering 50 (2001)1331–1346. 5. A.V. Lyamin, S.W. Sloan, K. Krabbenhoft, M. Hjiaj, Lower bound limit analysis with adaptive remeshing, International Journal for Numerical Methods in Engineering 63(2005) 1961–1974. 6. H. Ciria, J. Peraire, Computation of upper and lower bounds in limit analysis using secondorder cone programming and mesh adaptivity, Proceedings of 9th ASCE Specialty Conference on Probabilistic Mechanics and Structural Reliability, 2004. 7. H. Ciria, J. Peraire, J. Bonet, Mesh adaptive computation of upper and lower bounds in limit analysis, International Journal for Numerical Methods in Engineering 75 (2008) 899–944. 8. J. Munoz, J. Bonet, A. Huerta, J. Peraire, Upper and lower bounds in limit analysis: adaptive meshing strategies and discontinuous loading, International Journal for Numerical Methods in Engineering 77 (2009) 471–501. 9. C.M. Martin, The use of adaptive finite-element limit analysis to reveal slip-line fields, Géotechnique Letters 1 (2011) 23–29. 10. C.M. Martin, Undrained collapse of a shallow plane-strain trapdoor, Géotechnique 59 (10) (2009) 855–863. 11. C.M. Martin, The use of adaptive finite-element limit analysis to reveal slip-line fields, G´eotechnique Lett. 1 (2011) 23–29. 12. D.S.K. Mana, S. Gourvenec, C.M. Martin, Critical skirt spacing for shallow foundations under general loading, J. Geotech. Geoenviron. Eng. 139 (9) (2013) 1554–1566. 13. S. W. Sloan, Geotechnical stability analysis, Géotechnique, 63 (7) (2013) 531–572. 14. V.C. Le, A stabilized discrete shear gap finite element for adaptive limit analysis of Mindlin-Reissner plates, International Journal for Numerical Methods in Engineering 96 (2013) 231–246.

22

15. H. Nguyen-Xuan, G.R. Liu, An edge-based finite element method (ES-FEM) with adaptive scaled-bubble functions for plane strain limit analysis, Computer Methods in Applied Mechanics and Engineering 285 (2015) 877–905. 16. H. Nguyen-Xuan and T. Rabczuk, Adaptive selective ES-FEM limit analysis of cracked plane-strain structures. Frontiers in Civil Engineering 9 (4) (2015) 478–490. 17. H. Nguyen-Xuan, C.T. Wu, G.R. Liu (2016). An adaptive selective ES-FEM for plastic collapse analysis, European Journal of Mechanics A/Solids 58 (2016) 278-290. 18. Canh V. Le. Yield-stress based error indicator for adaptive quasi-static yield design of structures, Computers and Structures 171 (2016) 1-8. 19. R.C. Batra, K.I. Ko. An adaptive mesh refinement technique for the analysis of shear bands in plane strain compression of a thermoviscoplastic solid, Computational Mechanics 10 (1992) 369 - 379. 20. H. Samet, The quadtree and related hierarchical data structures, ACM Computing Surveys; 16 (1984) 187–260. 21. P. Hansbo, C. Lovadina, I. Perugia, G. Sangalli, A Lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes, Numerische Mathematik 100 (2005) 91–115. 22. N. Palle, JA. Dantzig, An adaptive mesh refinement scheme for solidification problems, Metallurgical and Materials Transactions A 27 (1996) 707–718. 23. M. Ainsworth, B. Senior, Aspects of an adaptive hp-finite element method: adaptive strategy conforming, approximation and efficient solvers, Computer Methods in Applied Mechanics and Engineering 150 (1997) 65–87. 24. P. Krysl, P. Grinspun, E. Schroder, Natural hierarchical refinement for finite element methods. International Journal for Numerical Methods in Engineering 56 (8) (2003) 1109– 1124. 25. P. Kagan, A. Fischer, P.Z. Bar-Yoseph, Mechanically based models: adaptive refinement for B-spline finite element, International Journal for Numerical Methods in Engineering 57 (2003) 1145–1175. 26. A. Tabarraei, N. Sukumar, Adaptive computations on conforming quadtree meshes, Finite Elements in Analysis and Design 41 (2005) 686–702. 27. D.W. Spring, S.E. Leon, G.H. Paulino, Unstructured polygonal meshes with adaptive refinement for the numerical simulation of dynamic cohesive fracture, International Journal of Fracture 189 (2014) 33–57.

23

28. S. E. Leon, D. W. Spring, G. H. Paulino, Reduction in mesh bias for dynamic fracture using adaptive splitting of polygonal finite elements, International Journal for Numerical Methods in Engineering 100 (2014) 555–576. 29. H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, Polygonal finite elements for finite elasticity, International Journal for Numerical Methods in Engineering 101 (2015) 305–328. 30. K.Y. Sze, N. Sheng, Polygonal finite element method for nonlinear constitutive modeling of polycrystalline ferroelectrics, Finite Elements in Analysis and Design 42 (2) (2005) 107– 29. 31. A. Simone, C.A. Duarte, E. Van der Giessen, A generalized finite element method for polycrystals with discontinuous grain boundaries, International Journal for Numerical Methods in Engineering 67 (8) (2006) 1122–45. 32. A. Menk, S. Bordas, Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals, International Journal for Numerical Methods in Engineering 83 (7) (2010) 805–28. 33. C. Talischi, G.H. Paulino, A. Pereira, I.F.M Menezes, Polygonal finite element for topology optimization: a unifying paradigm, International Journal for Numerical Methods in Engineering 82 (6) (2007) 671–98. 34. C. Talischi, G.H. Paulino, A. Pereira, I.F.M. Menezes, PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes, Journal of Structural and Multidisciplinary Optimization 45 (3) (2012) 329-357. 35. A.L. Gain, G.H. Paulino, L.S. Duarte, I.F.M. Menezes, Topology optimization using polytopes, Computer Methods in Applied Mechanics and Engineering 293 (2015) 411-430. 2015. 36. C. Talischi, A. Pereira, G.H. Paulino, I.F.M Menezes, M.S. Carvalho, Polygonal finite elements for incompressible fluid flow, International Journal for Numerical Methods in Fluids 74 (2014) 134–151. 37. A. Tabarraei, N. Sukumar, Extended finite element method on polygonal and quadtree meshes, Computer Methods in Applied Mechanics and Engineering 197 (5) (2008) 425–438. 38. J.E. Bishop. Simulating the pervasive fracture of materials and structures using randomly closed packed Voronoi tessellations, Computational Mechanics 44 (4) (2009) 455–471. 39. E.T. Ooi, C.M. Song, F. Tin-Loi, Z.J. Yang, Polygon scaled boundary finite elements for crack propagation modeling, International Journal for Numerical Methods in Engineering 91 (3) (2012) 319–342.

24

40. J.C. Nagtegaal, D.M. Parks, J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Computer Methods in Applied Mechanics and Engineering 4 (1974) 153–177. 41. A. Capsoni, L. Corradi, A finite element formulation of the rigid-plastic limit analysis problem, International Journal for Numerical Methods in Engineering 40 (1997) 2063–2086. 42. S.W. Sloan, P.W. Kleeman, Upper bound limit analysis using discontinuous velocity fields, Computer Methods in Applied Mechanics and Engineering 127 (1995) 293–314. 43. K. Krabbenhøft, A.V. Lyamin, M. Hjiaj, S.W. Sloan, A new discontinuous upper bound limit analysis formulation, International Journal for Numerical Methods in Engineering 63 (2005) 1069–1088. 44. A.V. Lyamin, S.W. Sloan, Upper bound limit analysis using linear finite elements and nonlinear programming, International Journal for Numerical and Analytical Methods in Geomechanics 26 (2002) 181–216. 45. A. Makrodimopoulos, C.M. Martin, Upper bound limit analysis using simplex strain elements and second-order cone programming, International Journal for Numerical and Analytical Methods in Geomechanics 31 (2007) 835–865. 46. W.T. Koiter, General theorems for elastic plastic solids, Proceedings of Solid Mechanics (edited by Sneddon I. N. and Hill R.), Nord-Holland, Amsterdam, 1960. 47. R. Hill, On discontinuous plastic states, with special reference to localized necking in thin sheets, Journal of the Mechanics and Physics of Solids 1 (1952) 19−30. 48. D.J.F. Ewing, C.E. Richards, The yield-point loading of singly-notched pin loaded tensile strips, Journal of the Mechanics and Physics of Solids 22 (1974) 27−36. 49. G.G. Chell, Elastic-plastic fracture mechanics. In Developments in Fracture Mechanics-1, Applied Science Publishers: London (1979) 67–105. 50. A.G. Miller, Review of limit loading of structures containing defects, International Journal of Pressure Vessels and Piping 32 (1988) 197−327. 51. E. Melan, Theorie statisch unbestimmter Systeme aus ideal plastischem Baustoff. Sitzber. Akad. Wiss. Wien IIa 145 (1936) 195–218. 52. A.M. Yan, H. Nguyen-Dang, Limit analysis of cracked structures by mathematical programming and finite element technique, Computational Mechanics 24 (1999) 319–333. 53. I.A. Khan, A.K. Ghosh, A modified upper bound approach to limit analysis for plane strain deeply cracked specimens, International Journal of Solids and Structures 44 (2007) 3114– 3135.

25

54. I.A. Khan, V. Bhasin, J. Chattopadhyay, R.K. Singh, K.K. Vaze, A.K. Ghosh, An insight of the structure of stress fields for stationary crack in strength mismatch weld under plane strain mode–I loading–Part II: Compact tension and middle tension specimens, International Journal of Mechanical Sciences 87 (2014) 281–296. 55. Drucker DC, Prager W. Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics 1952; 10:157–165. 56. G.R. Liu, T. Nguyen-Thoi, K.Y. Lam, An edge-based smoothed finite element method (ESFEM) for static, free and forced vibration analyses of solids, Journal of Sound and Vibration 320 (2009) 1100–1130. 57. C. Talischi, A. Pereira, G.H. Paulino, I.F.M. Menezes, M.S. Carvalho, Polygonal finite elements for incompressible fluid flow, International Journal for Numerical Methods in Fluids 74 (2014) 134–151. 58. E.A. Malsch, G. Dasgupta, Shape functions for polygonal domains with interior nodes, International Journal for Numerical Methods in Engineering 61 (8) (2004) 1153–1172. 59. D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337-344. 60. E.L. Wachspress, A rational finite element basis. Academic Press, New York, USA, 1975. 61. R. Sibson, A vector identity for the Dirichlet tesselation, Mathematical Proceedings of the Cambridge Philosophical Society 87 (1980) 151–155. 62. J. Warren, Barycentric coordinates for convex polytopes, Advances in Computational Mathematics 6 (1) (1996) 97–108. 63. H. Hiyoshi, K. Sugihara, Two generalizations of an interpolant based on Voronoi diagrams, International Journal of Shape Modeling 5 (2) (1999) 219–231. 64. M. Floater, Mean value coordinates, Computer Aided Geometric Design 20 (2003) 19–27. 65. E.A. Malsch, J.J. Lin, G. Dasgupta, Smooth two dimensional interpolants: a recipe for all polygons, Journal of Graphics Tools 10 (2) (2005) 27–39. 66. K. Hormann, N. Sukumar, Maximum entropy coordinates for arbitrary polytopes, Computer Graphics Forum 27 (5) (2008) 1513–1520. 67. J. Manson, S. Schaefer, Moving least squares coordinates, Computer Graphics Forum 29 (5) (2010) 1517– 1524. 68. X.Y. Li, S.M. Hu, Poisson coordinates, IEEE Transactions on Visualization and Computer Graphics 19 (2013) 344–352. 69. X.Y. Li, T. Ju, S.M. Hu, Cubic mean value coordinates, ACM Transactions on Graphics 32 (4) (2013) 126. 26

70. M. Meyer, A. Barr, H. Lee, M. Desbrun, Generalized barycentric coordinates on irregular polygons, Journal of Graph Theory 7 (2002) 13–22. 71. M.S. Floater, A. Gillette, N. Sukumar, Gradient bounds for Wachspress coordinates on polytopes, SIAM Journal of Numerical Analysis 52(1) (2014) 515–532. 72. N. Sukumar, E. Malsch, Recent advances in the construction of polygonal finite element interpolants, Archives of Computational Methods in Engineering 13 (1) (2006) 129–163. 73. M. S. Floater, K. Hormann, and G. Kós. A general construction of barycentric coordinates over convex polygons. Advances in Computational Mathematics, 24 (1–4) (2006) 311–331. 74. N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, International Journal for Numerical Methods in Engineering 61 (12) (2004) 2045–2066. 75. Mosek, The MOSEK optimization toolbox for MATLAB manual, http://www.mosek.com, Mosek ApS, Version 5.0 Edition, 2009. 76. W. Dorfler, A convergent adaptive algorithm for Poisson’s equation, SIAM Journal on Numerical Analysis 33 (1996) 1106–1124. 77. M. Yerry, M. Shepard, A modified quadtree approach to finite element mesh generation, IEEE Computer Graphics and Applications 3 (1983) 39–46. 78. H. Samet, The quadtree and related hierarchical data structures, ACM Computing Surveys; 16 (1984) 187–260. 79. A. Tabarraei, N. Sukumar, Adaptive computations on conforming quadtree meshes, Finite Elements in Analysis and Design 41 (2005) 686–702. 80. C. Baiocchi, F. Brezzi, L. Franca, Virtual bubbles and the Galerkin least-squares method, Computer Methods in Applied Mechanics and Engineering 105 (1) (1993) 125–141. 81. E.T. Filipov, J. Chun, G.H. Paulino, J. Song, Polygonal multiresolution topology optimization

(PolyMTOP)

for

structural

dynamics,

Journal

of

Structural

and

Multidisciplinary Optimization, (2015) DOI 10.1007/s00158-014-1211-y 82. C. Talischi, G.H. Paulino, A. Pereira, I.F.M. Menezes, PolyMesher: A general-purpose mesh generator for polygonal elements written in MATLAB, Structural and Multidisciplinary Optimization 45 (3) (2011) 309–328. 83. Y.-J. Kim, D.-J. Shim, N.-S. Huh, Y.-J. Kim, Plastic limit pressures for cracked pipes using finite element limit analyses, International Journal of Pressure Vessels and Piping 79 (2002) 321–30. 84. Liu GR, Nguyen-Xuan H, Nguyen-Thoi T. A variationally consistent FEM (VCFEM) for solution bounds and nearly exact solution to mechanics problems using quadrilateral elements. Int J Numer Methods Eng 2011;85(4):403–536. 27

85. C.V. Le, H. Askes, M. Gilbert, A locking-free stabilized kinematic EFG model for plane strain limit analysis, Computers and Structures 106–107 (2012) 1– 8.

28

List of Figures

a) Triangles

b) Perpendicular distances

c) 3D view

d) Contour line

Fig. 1. Shape function for a boundary node of the hexagonal element using the Wachspress coordinates.

29

a)

b)

c) Fig. 2. 3D view and contour bubble shape functions: (a) triangle, (b) quadrilateral and (c) hexagon.

30

Fig. 3. Edge-based sub-domain associated with edge k.

31

a)

b)

32

c)

d)

Fig. 4. Polytree meshes (a)-(c) and its representative tree (d).

33

Fig. 5. Mapping from a regular hexagonal element (right) to an arbitrary pentagonal element with one side-node (left).

Fig. 6. Sequence of steps for constructing children polygonal elements from an underlying triangular mesh.

34

a) 1st iteration

b) 2nd iteration

c) 3rd iteration

Fig. 7. Polytree mesh under three iterations by using (left) connecting the centroid of all the triangular elements circumventing a particular node to the centroid of the element, and (right) centroid Voronoi tessellation (CVT).

35

a)

b)

c)

Fig. 8. Subdivisions of three arbitrary pentagonal elements using the Centroidal Voronoi tessellations (CVTs): a) initial seeds within domain of element, b) polygonal mesh at first iteration, and c) polygonal mesh after few iterations.

36

a) Non-matching

b) Subdivision for six nodes polygon element and its modified positions.

Fig. 9. Non-matching mesh and modified mesh.

37

Fig. 10. A sequence of subdivisions of arbitrary polygonal elements using the Centroidal Voronoi tessellations (CVTs).

38

Fig. 11. Flowchart of solving limit analysis problems using polytree meshes.

39

Fig. 12. Cook’s membrane problem.

a) Mesh 1st (137 nodes and 45 elements)

b) Mesh 2nd (242 nodes and 80 elements)

c) Mesh 3rd (488 nodes and 162 elements)

d) Mesh 4th (950 nodes and 316 elements)

Fig. 13. Domain discretization of the Cook’s membrane using four polygonal meshes.

40

a) W-Poly

b) Present (bW-Poly)

c) MINI

Fig. 14. An illustration of the level line of pressure distribution for Cook’s membrane in nearly incompressible limit.

41

log10 (Relative error of vertical displacement)

0 ‐0.5 ‐1 ‐1.5 ‐2 ‐2.5 ‐3 ‐3.5 0.2

MINI (r=1.47) Present (bW‐Poly) (r=2.68) 0.3

0.4

0.5

0.6

0.7

0.8

0.9

log10 (h)

a)

log10 (Relative error of strain energy [%])

2 1.5 1 0.5 0 ‐0.5 ‐1 ‐1.5 0.2

MINI (r=1.45) Present (bW‐Poly) (r=2.4) 0.3

0.4

0.5 0.6 log10 (h)

0.7

0.8

0.9

b) Fig. 15. The relative errors of the vertical displacement (a) and strain energy (b) for Cook’s membrane in nearly incompressible limit.

42

Fig. 16. A notched tensile specimen: fully and its quarter model.

W‐Poly Present (bW‐Poly) Reference

1.35

Load factor

1.3

1.25

1.2

1.15

1.1

0

5000

10000 Nvar

15000

Fig. 17. Notched tensile specimen (a = 1/ 2) : The convergence of collapse multiplier

versus number of variables.

43

a) W-Poly

b) bW-Poly

Fig. 18. Collapse mechanism of the notched tensile specimen (a = 1/ 2) .

a) Initial mesh

b) 1st step

44

d) 3rd step

c) 2nd step

e) final step Fig. 19. Mesh refinements for the notched tensile specimen with a = 1/ 2 .

45

Relative error [%]

10

Uniform meshes Polytree meshes

1

10

0

10

3

10

4

Nvar

a) a = 1/ 3

Relative error [%]

10

Uniform meshes Polytree meshes

1

10

0

10

3

10

4

Nvar

b) a = 1/ 2

46

Uniform meshes Polytree meshes 1

Relative error [%]

10

10

0

10

3

10

4

Nvar

c) a = 2 / 3

Fig. 20. The relative error of collapse multiplier of the notched tensile specimen.

47

a)

b) Fig. 21. Fine polytree meshes of the notched tensile specimen with (a) a = 1/ 3 and (b) a = 2 / 3.

48

Fig. 22. Central cracked plate: system and its quarter model.

49

a) x = 1/ 4

b) x = 1/ 2

50

c) x = 3 / 4 Fig. 23. A fine polytree mesh of the central cracked plate with x = 1/ 4 , x = 1/ 2 and x = 3/ 4.

51

1.02 Analytical solution Uniform meshes Polytree meshes

1 0.98

Load factor

0.96 0.94 0.92 0.9 0.88 0.86

0

0.5

1

1.5

2

Nvar

x 10

4

a)

Analytical solution Uniform meshes Polytree meshes

0.72 0.7

Load factor

0.68 0.66 0.64 0.62 0.6 0.58 0

0.5

1

1.5 Nvar

2 x 10

4

b)

52

0.45 Analytical solution Uniform meshes Polytree meshes

Load factor

0.4

0.35

0.3

0

2000

4000

6000 Nvar

8000

10000

12000

c)

Fig. 24. The convegence rate of collapse multiplier of the central cracked plate: (a) x = 1/ 4 ; (b) x = 1/ 2 ; and (c) x = 3 / 4 .

53

Fig. 25. Single-edge cracked plate: fully and its quarter model.

a) x = 1/ 4

54

b) x = 1/ 2

c) x = 3 / 4 Fig. 26. A fine adaptive mesh of the edge-cracked tensile specimen with x = 1/ 4, x = 1 / 2 and x = 3 / 4 .

55

0.92 Analytical solution (UB) Analytical solution (LB) Uniform meshes Polytree meshes

0.9

Load factor

0.88 0.86

0.84

0.82 0.8 0

0.5

1 Nvar

1.5

2 x 10

4

a) 0.54 Analytical solution (UB) Analytical solution (LB) Uniform meshes Polytree meshes

0.52 0.5

Load factor

0.48 0.46 0.44 0.42 0.4 0.38 0.36 0

0.5

1 Nvar

1.5

2 x 10

4

b)

56

0.12 Analytical solution Uniform meshes Polytree meshes

0.11

Load factor

0.1

0.09

0.08

0.07

0.06

0

0.5

1 Nvar

1.5

2 x 10

4

c)

Fig. 27. The convegence rate of collapse multiplier of the single-edge cracked plate: a) x = 1/ 4 ; b) x = 1/ 2 ; and c) x = 3 / 4 .

57

a)

b)

c)

Fig. 28. Plastic dissipation of the single-edge cracked plate: (a) x = 1 / 4 ; (b) x = 1 / 2 ; and (c) x = 3 / 4 .

Fig. 29. Cracked cylinders under internal pressure: fully and its half model.

58

a) a / h = 1/ 4

59

b) a / h = 1/ 2

b) a / h = 3 / 4

Fig. 30. Fine polytree meshes of the cracked cylinders under internal pressure.

60

a) a / h = 1/ 4

b) a / h = 1/ 2

61

c) a / h = 3 / 4

Fig. 31. Plastic dissipation of the cracked cylinders under internal pressure: (a) a / h = 1/ 4 ; (b) a / h = 1/ 2 ; and (c) a / h = 3 / 4 .

62

Chell Miller Yan & Nguyen‐Dang  Present

1

Load factor

0.8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

a/h

a) Comparison with analytically approximated solutions

1

Yan ‐ Q8 Kim et al. Present

0.9 0.8

Load factor

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

a/h

b) Comparison with numerical solutions

Fig. 32. Limit load factor of the cylinder with a longitudinal crack under internal pressure.

63

List of Tables

Table 1: The computational efficiency of the bW-Poly element for the notched tensile specimen (a = 1/2) λ+

1.2460

1.1922

1.1668

1.1478

1.1423

Variables

592

2257

7327

30517

69547

Iteration

12

15

16

17

18

Optimizer time (s)

1.65

2.42

6.22

30.55

102.04

Table 2: The relative error in the collapse multiplier of the notched tensile specimen. Uniform meshes N var

λ+

Polytree meshes

Relative error (%)

N var

λ+

Relative error (%)

a = 1/ 3 (Reference value λ * = 0.9241)

592

1.0143

9.76

592

1.0143

9.76

2257

0.9654

4.47

1667

0.9497

2.78

7327

0.9505

2.86

4092

0.9367

1.37

30517

0.9350

1.18

11462

0.9304

0.68

a = 1 / 2 (Reference value λ * = 1.1316)

592

1.2469

10.19

592

1.2469

10.19

2257

1.1926

5.40

1742

1.1638

2.85

7327

1.1669

3.12

4247

1.1474

1.40

30517

1.1480

1.45

11822

1.1390

0.65

a = 2 / 3 (Reference value λ * = 1.3833)

592

1.5616

12.89

592

1.5616

12.89

2257

1.4703

6.29

1972

1.4313

3.47

7327

1.4197

2.63

4732

1.4078

1.77

17497

1.4041

1.50

13072

1.3898

0.47

64

Table 3: Notched tensile specimen: Comparison with literature solutions Approach Kinematic

Authors Le et al. [85]

a = 1/ 3

a = 1/ 2

a = 2/3

0.9412

1.1535

1.4090

Ciria et al. [7] •

Uniform meshes

--

1.1495

--

•

Adaptive meshes

--

1.1390

--

Present element •

Uniform meshes

0.9434

1.1480

1.4041

•

Polytree meshes

0.9304

1.1390

1.3898

Mixed

Christiansen and Andersen [1]

0.9276

1.1358

1.3884

Static

Ciria et al. [7]

Reference

•

Uniform meshes

--

1.131

--

•

Adaptive meshes

--

1.132

--

0.9241

1.1316

1.3833

Christiansen and Andersen [1]

Table 4: The relative error in the collapse multiplier of the central cracked plate. Uniform meshes N var

λ+

Polytree meshes

Relative error (%)

N var

λ+

Relative error (%)

x = 1/ 4 (Analytical value λ * = 0.8660)

487

1.0117

16.82

487

1.0117

16.82

2257

0.9429

8.88

1307

0.9088

4.94

4507

0.9159

5.76

3142

0.8863

2.34

22507

0.8896

2.72

24592

0.8717

0.65

x = 1/ 2 (Analytical value λ * = 0.5774)

487

0.7199

24.69

487

0.7199

24.69

2257

0.6469

12.04

1022

0.6171

6.88

4507

0.6267

8.54

2202

0.5967

3.36

12002

0.6080

5.31

5597

0.5874

1.73

22507

0.5965

3.32

15932

0.5803

0.50

x = 3 / 4 (Analytical value λ * = 0.2887)

65

487

0.4392

52.14

487

0.4392

52.14

2257

0.3610

25.06

927

0.3312

14.72

4507

0.3331

15.39

1537

0.3084

6.85

7507

0.3228

11.82

3037

0.2979

3.21

12002

0.3160

9.48

6737

0.2936

1.69

Table 5: The relative error in the collapse multiplier of the single-edge cracked plate. Uniform meshes N var

λ+

Polytree meshes

Relative error (%)

N var

λ+

Relative error (%)

x = 1/ 4 (Analytical value: λ + = 0.8297, λ − = 0.7952)

592

0.9142

10.19

592

0.9142

10.19

2257

0.8806

6.13

2012

0.8543

2.98

6007

0.8608

3.75

5942

0.8369

0.87

19502

0.8450

1.85

17452

0.8306

0.11

x = 1/ 2 (Analytical value: λ + = 0.3725, λ − = 0.3660)

592

0.5279

41.73

592

0.5279

41.73

2257

0.4302

15.50

1662

0.4158

11.63

6007

0.4047

8.64

3882

0.3824

2.66

12007

0.3961

6.35

11102

0.3779

1.46

x = 3 / 4 (Analytical value: λ * = 0.0644)

592

0.1123

74.42

592

0.1123

74.42

2257

0.0834

29.58

1132

0.0786

22.02

6007

0.0775

20.30

2187

0.0718

11.56

12007

0.0732

13.64

5067

0.0666

3.35

19502

0.0716

11.15

13897

0.0659

2.30

66

Highlights •

We propose a novel adaptive finite element scheme for limit analysis of cracked planestrain structures.

•

The method relies on a volumetric locking-free polygonal finite element formulation.

•

An adaptive primal-mesh strategy is guided by the L2 -norm-based indicator of strain rates.

•

A novel polytree date structure over polygons is addressed.

•

The present approach reaches high accuracy with low computational cost.

67