An adaptive isogeometric-meshfree coupling approach for the limit analysis of cracked structures

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

An adaptive isogeometric-meshfree coupling approach for the limit analysis of cracked structures Jiazhao Huanga,b, Nhon Nguyen-Thanhb, Weidong Lia, Kun Zhoua, a b

T

⁎

School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore School of Electrical and Electronic Engineering, Nanyang Technological University, 50 Nanyang Avenue, 639798, Singapore

A R T I C LE I N FO

A B S T R A C T

Keywords: Plasticity Adaptive IGA Meshfree methods Fracture analysis Limit analysis

An adaptive isogeometric-meshfree (AIMF) coupling approach is proposed to predict the limit load factors of cracked structures. The concept of the present approach, which relies on forming an equivalence between the isogeometric analysis (IGA) and moving least-squares meshfree method, is developed based on the reproducing conditions, resulting in a unified formulation of the basis functions. Thereby, the refinement of IGA can be implemented in a straightforward meshfree manner for limit analysis problems. The adaptivity of refinement is realized by adopting an indicator of plastic dissipation to automatically identify the material regions associated with the dissipated work greater than a predefined threshold. Subsequently, the marked meshes are refined through the insertion of linear reproducing points. Enrichment functions are further introduced to the present approach for crack modelling. The resulting optimization formulation of limit analysis is re-expressed in the form of second-order cone programming which can be effectively tackled by the interior-point solvers. Through a series of numerical examples, the present approach has been proven to achieve both high convergence rates and accurate simulation results.

1. Introduction Limit analysis has been established as an effective tool to investigate the behaviour of structures under proportional loading in the plastic regime [1–3]. The application of limit analysis to the estimation of fracture failure of ductile cracked structures has been explored through analytical approaches such as the slip-line method [4]. Since it is not feasible to determine the exact limit loads for most engineering configurations, numerical limit analysis procedures based on the finite element method (FEM) have been developed as approximation methods [5–7]. Extended FEM (XFEM) has been further developed and proposed for modelling cracked structures without the tedious remeshing over discontinuity surfaces [8,9]. Another key feature of numerical solutions for limit analysis problems is mathematical programming. It has been shown that the yield criteria can be represented as an intersection of cones[10]. Limit analysis problems are thus cast in the form of a second-order cone programme (SOCP) which can be effectively solved by the primal-dual interior-point method [11,12]. As a potential alternative to FEM, isogeometric analysis (IGA) seamlessly integrates computer-aided designs and numerical simulations through non-uniform rational B-splines (NURBS) [13,14]. By combining isogeometric elements with XFEM, Ghorashi et al. [11]

⁎

further developed the extended isogeometric analysis (XIGA) for the fracture analysis of cracked structures. It has been proven that IGA preserves the exactness of the domains of conic sections and offers highorder continuity across elements [15–17]. Nguyen-Xuan et al. [18] integrated XIGA with SOCP for the limit analysis of cracked structures. Nevertheless, the XIGA-based limit analysis suffers from its significant computational cost due to the excessive overhead of control points with increasing levels of refinement. Moreover, localized plastic deformations could also lead to a slow convergence rate for simulations [19]. Therefore, an adaptive local refinement strategy should be developed for the IGA-based limit analysis. While the implementation of adaptive IGA mesh refinement is complicated due to the tensor product nature of its basis functions [13,20], it is much simpler to perform local refinement in a meshfree manner, especially for the fracture analysis [21–24]. Recently, Wang and Zhang [25] developed a consistently coupled isogeometric meshfree method to achieve a smoothing transition between IGA and reproducing kernel meshfree sub-domains based on the reproducing conditions. The geometry exactness can be well preserved by IGA basis functions, and local refinement can be easily implemented in the meshfree sub-domains. Zhang and Wang [26] further unified the isogeometric and meshfree shape functions across the entire domain, which allowed the local refinement of IGA to be implemented in a

Corresponding author. E-mail address: [email protected] (K. Zhou).

https://doi.org/10.1016/j.tafmec.2019.102426 Received 4 September 2019; Received in revised form 29 October 2019; Accepted 27 November 2019 Available online 16 December 2019 0167-8442/ © 2019 Elsevier Ltd. All rights reserved.

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

with respect of a(ξ ) to 0 as

direct meshfree fashion. The application of this approach to solving limit analysis problems should be further explored. In this paper, an adaptive isogeometric-meshfree (AIMF) approach is presented in tandem with an indicator of plastic dissipation to evaluate the limit load factors of cracked structures. These structures are assumed to be made of elastic-perfectly-plastic materials and subject to the von Mises yield criterion. The concept of the present approach relies on developing an equivalence between the IGA and the MLS meshfree basis functions, which preserves the geometric exactness of isogeometric discretizations and offers high flexibility for meshfree local refinement. The indicator that guides the adaptive local refinement is defined based on the amount of plastic dissipation across the entire domain. In this way, the regions containing dissipated work greater than a predefined threshold could be adaptively located and refined. Enrichment functions are applied to model the pre-existing cracks in structures. The resulting non-smooth optimization problem is then reformulated by minimizing the sum of Euclidean norms which can be effectively computed by SOCP [27–29]. The contents of this paper are arranged as follows: the development of an AIMF coupling approach is introduced in Section 2. The adaptive local refinement strategy is demonstrated in Section 3. Section 4 describes the solution procedure of the plastic collapse problem. A series of numerical examples are investigated in Section 5, and the conclusions are finally provided in Section 6.

a (ξ ) = A (ξ )−1B (ξ ) U,

(6)

}T

where U = {u1, u2, ⋯, unm is the vector set of nodal parameters, and matrices A and B can be obtained respectively by the equations: nm

⌢

∑ W (ξ − ξI ) p (ξI ) pT (ξI ),

A (ξ ) =

(7)

I=1

⌢

⌢

⌢

B (ξ ) = [W (ξ − ξ1 ) p (ξ1) W (ξ − ξ2 ) p (ξ2) ⋯ W (ξ − ξnm ) p (ξnm )]. (8) By substituting Eqs. (6)–(8) into Eq. (1), the displacement field given by

∑ ΨI (ξ ) uI (ξ ), I=1

m

∑ pJ (ξ )(A−1 (ξ ) B (ξ ))JI = pT (ξ )(A (ξ )−1B (ξ ))I .

ΨI (ξ ) =

J =1

nm I=1

The B-spline is developed from a knot set Ξ = {ξ1, ξ2, ⋯, ξn + p + 1} , where p and n are the polynomial degree and the total number of control points respectively. The B-spline basis functions Ni, p can thus be defined in a recursive form as shown in [13]

Since the meshfree formulation employed in this work is derived from the MLS approximation, the displacement field uh for the twodimensional (2D) domain ξ = [ξ , η] is given by [30,31]

{10

Ni, p (ξ ) =

ξi + p + 1 − ξ ξ − ξi Ni + 1, p − 1 (ξ ), p ⩾ 1. Ni, p − 1 (ξ ) + ξi + p + 1 − ξi + 1 ξi + p − ξi

(2)

∑ Ni (ξ ) Pi . i=1

n

S¯ (ξ ) =

where nm is the number of nodes within the support domain around ξ and uI denotes nodal parameter. The weight function Ẇ in Eq. (4) determines the degree of smoothness for meshfree basis functions, and its cubic splines form is given by 2

⎧ 3 − 4r 2 + 4r 3 ⎪ ̇ W (ξ − ξI ) = 4 − 4r + 4r 2 − ⎨3 ⎪0 ⎩

(15)

where the 2D NURBS basis functions Ri can be defined in terms of Ni and weights wi as

(4)

I=1

∑ Ri (ξ ) Pi , i=1

⌢

∑ W (ξ − ξI )[pT (ξI ) a (ξ ) − uI ]2 ,

(14)

A NURBS surface S¯ can be further defined in a similar form as [13]

(3)

In order to obtain the expression of a , the weighted L norm J of displacement fields should be first calculated by [32]

J=

(13)

n

S (ξ ) =

2

nm

(12)

The expression for a B-spline surface on the 2D domain ξ = [ξ , η] can be given by the control points Pi and B-spline basis functions by [13]

When the order p equals 2, the quadratic form of p is computed as

ξη ,

ξi ⩽ ξ < ξi + 1 otherwise

Ni,0 (ξ ) =

(1)

where a is the vector of coefficients and p is the pth order basis vector defined by

p(ξ ) = {1, ξ , η ,

(11)

2.2. Isogeometric approximation

2.1. MLS meshfree approximation

η2}T .

(10)

To reproduce the polynomial of complete order, the 2D consistency condition for meshfree basis functions is calculated by

In this section, a brief review of both MLS meshfree and isogeometric approximations is given first, followed by the development of AIMF coupling basis functions. It should be noted that IGA basis functions and meshfree shape functions are both denoted as the basis functions in the rest sections.

ξ 2,

(9)

where the MLS basis function ΨI is defined as

2. Formulation of AIMF coupling basis functions

p(ξ ) = {1, ξ , η , ξ 2, ξη , η2, ⋯, ξ p, ⋯, η p}T .

is

nm

uh (ξ ) = Ψ(ξ ) U =

∑ ΨI (ξ ) p (ξI ) = p (ξ ).

uh (ξ ) = pT (ξ ) a (ξ ),

uh

Ri (ξ ) =

Ni (ξ ) wi . n ∑^i = 1 N ^i (ξ ) w ^i

(16)

2.3. Reproducing conditions for isogeometric basis functions

r ⩽ 0.5 4r 3 3

0.5 < r ⩽ 1

with r =

|ξ − ξI | I d max

The 2D basis functions based on the reproducing conditions can be derived from the product of 1D B-spline basis functions as [25]

,

r⩾1

n

ξ a1 ηa2 =

(5)

m

∑ ∑ Ni,p (ξ ) Mj,q (η)(ξi[a1] )a1 (η[j a2 ] )a2 with a1 + a2 ⩽ l, i=1 j=1

I denotes the size of the rectangular or circular support where d max domain around the AIMF coupling node ξI [32]. Therefore, Eq. (4) can be reformulated by setting the derivatives of J

(17)

where N and M represent one-dimensional (1D) B-spline basis functions. 2

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 1. A 2D AIMF coupling shape function: (a) meshes and nodes and (b) basis functions.

By replacing Ni, p (ξ ) Mj, q (η) with NIpq (ξ ) , Eq. (17) can be further simplified into the following form as [25] nb

∑

NIpq (ξ ) p (ξI[l] )

3.1. Indicator of plastic dissipation Plane stress problems under the von Mises yield criterion are considered in this work, which further yields the plastic dissipation function D as [33]

= p (ξ ), (18)

I=1

D (ε )̇ = σ0

where nb is the total number of isogeometric basis functions and the 2D reproducing point vector is defined as shown in [25]

p(ξI[l] ) = {1, ξI[1] , ηI[1] , (ξI[2] )2 , ξI[1] ηI[1] , (ηI[2] )2, ⋯, (ξI[p] ) p, ⋯, (ηI[q] )q}T .

ε Ṫ Θε ̇ dΩ,

(19)

Θ=

1 ⎡4 3 0⎤ 3 4 0 . 3 ⎢0 0 1⎥ ⎣ ⎦

= ξi + 1 ξi + 2 , respectively. The normalized support size is set to 1.5 so that each node is influenced by three basis functions. The reproducing kernel formulation can further be derived as [26]

⌢

dicator to guide the refinement [19,29]. The total plastic dissipation De corresponding to the reproducing points within an element e can be obtained by

⌢

(20)

⌢

De =

where the moment matrix C can be expressed by [26]

ne

∑ Di , i=1

C (ξ ) =

⌢

∑ p (ξJ[l] ) pT (ξJ[l] )W(ξ − ξJ[1] ). J =1

ΦJ (ξ ) wJ . n ∑J¯ = 1 Φ J¯ (ξ ) w J¯

(25)

where nk is the number of nodes associated with the element k. The mesh density n¯ k = nk/Ak is defined to describe the number of nodes within a unit element area Ak . The averaged plastic dissipation D¯ k of the element k is calculated by

(21)

ξJ[1] in Eq. (21) is the linear reproducing points which can be regarded as additional meshfree field nodes. Zhang and Wang [26] proved that 2D B-spline basis functions NI are equivalent to ΦJ in Eq. (20) when the support size is (p + 1)/2. Reproducing conditions for a 2D NURBS surface S¯J can be further expressed in terms of ΦJ(ξ) and weights wJ as [26] S¯J (ξ ) =

(24)

Since plastic dissipation is concentrated within the regions of the material containing high strain rates, the plastic dissipation D¯ of each background element on the parametric domain is employed as an in-

ξi[2]

nm

(23)

in which σ0 is the yield stress and the matrix Θ is defined as

When quadratic basis functions are considered, the reproducing point vector is p(ξi[l] ) = {1, ξi[1] , (ξi[2] )2}T , where the linear and quadratic reproducing points are given as ξi[1] = (ξi + 1 + ξi + 2) 2 and

ΦJ (ξ ) = pT (ξJ[l] ) C−1 (ξ ) p (ξ )W(ξ − ξJ[1] ),

∫Ω

D D¯ k = k . n¯ k

(26)

A threshold γ* is further introduced herein to activate the refinement for a set of elements when their corresponding averaged mesh intensity D¯ k is larger than the threshold γ* .

(22)

3.2. Local mesh refinement Referring to the local mesh refinement presented by Zhang and Wang [26], Fig. 1(a) shows an initial quadratic mesh on the 2D domain with a set of knot vectors Ξ = {0, 0, 0, 1 4, 1 2, 3 4, 1, 1, 1} along each direction. By utilizing the quadratic basis functions derived from Eqs. (19)–(21), a set of linear reproducing points ξ [1] = {0, 1 8, 3 8, 5 8, 7 8, 1} could be generated. Its corresponding AIMF coupling basis functions which match the

3. Adaptive refinement procedure Adaptive local refinement techniques are preferred for limit analysis as substantial changes in plastic strain rates are found along the yield lines. An indicator of plastic dissipation is developed to guide the refinement. 3

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 2. Refinement of a 2D AIMF coupling surface: (a) meshes and nodes and (b) basis functions.

Fig. 3. (a) Full mode and (b) its quarter mode with symmetric boundary conditions of a square plate with a central hole.

Fig. 4. (a) Initial mesh and (b) mesh after 3 steps of refinement for a square plate with a central hole (r = 0.2 l).

4

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Table 2 Limit load factors λ of a square plate with various radius-to-length ratios r/l under uniaxial tension. r/l

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Methods Heitzer [44]

Tran et al. [42]

Tran [45]

Nguyen‐Xuan et al. [18]

Present method

0.8951 0.7879 0.6910 0.5720 0.4409 0.2556 0.1378 0.0565 0.0193

0.8932 0.7967 0.6930 0.5760 0.4011 0.2429 0.1277 0.0521 0.0133

0.9017 0.8015 0.7022 0.5914 0.4012 0.2425 0.1254 0.0523 0.0123

0.9005 0.8010 0.7016 0.5933 0.4168 0.2532 0.1315 0.0528 0.0119

0.8990 0.8001 0.7002 0.5872 0.4006 0.2410 0.1250 0.0509 0.0122

Fig. 2(b) further proves that a smooth transition of the basis functions between the refined and normal elements is achieved. 4. An AIMF coupling approach for the limit analysis of cracked structures

Fig. 5. The convergence of limit load factor solutions for a plate with a central hole obtained by different numerical methods.

The limit analysis is an effective tool of limit load factors which is of key importance in evaluating the load-carrying capacity of structures. Enrichment functions are incorporated into the upper bound limit analysis of cracked structures. The limit analysis problem is further transformed into a non-linear optimization problem which can be computed by a SOCP solver.

Table 1 Limit load factors λ of a square plate with a central hole (r = 0.2 l) under various loading conditions with σ1 = σ0. Methods

Gaydon and McCrum [37] Le et al. [5] Groβ-Weege [38] Chen et al. [39] Zouain et al. [40] Zhang et al. [41] Tran et al. [42] Nguyen‐Xuan et al. [43] Nguyen‐Xuan et al. [18] Present method

Loading conditions σ2 = σ0

σ2 = 0.5σ0

σ2 = 0

– – 0.882 0.8740 0.8940 0.8890 0.8960 0.8940 0.8955 0.8936

– – 0.891 0.899 0.9110 0.8980 0.9120 0.9110 0.9112 0.9088

0.8000 0.8006 0.7820 0.7980 0.8030 0.7840 0.8050 0.8020 0.8010 0.8001

4.1. Kinematic formulation A rigid-perfectly plastic body containing discontinuous parts is defined on a 2D problem domain Ω ∈ 2 with boundary conditions of Γ = Γu̇ ∪ Γt ∪ Γc and Γu̇ ∩ Γt ∩ Γc = ∅, where Γu̇ , Γt and Γc represent prescribed the boundary, free boundary and crack surface, respectively. The body is subject to body forces f on Γu̇ prescribed by the displacement velocity u̇ and surface traction g along Γt . u̇ = [u̇ v]̇ T is adopted to represent the plastic flow subject to a space V of kinematically admissible velocity fields. Consequently, the equilibrium equation can be given by [34]

Win (σ , u̇) = Wex (u̇) , ∀ u̇ ∈ V ,

(27)

where the internal work rate Win for σ and plastic velocity fields u̇ are given in a bilinear form as

Win (σ , u̇ ) =

∫Ω σ: ε ̇ (u̇ ) dΩ,

(28)

the external work rate Wex can be expressed in a linear form as

Wex (u̇ ) =

∫Ω f·u̇ dΩ + ∫Γ g·u̇ dΓ,

and space V can be defined by a Hilbert space

V = {u̇ ∈

(29)

t

(H1 (Ω))2 ,

u̇ = u̇ ∈ on Γu̇ }.

H1 (Ω)

as (30)

A convex set S consisting of the statically admissible stress field is further defined as

S = {σ ∈ E ; |ψ (σ ) ⩽ 0},

Fig. 6. Plastic dissipation distribution associated to σ2 = 0.

(31)

where E is a space of symmetric stress tensor and ψ represents the convex yield function. The plastic strain rates ε ̇ related to the associated flow rule can thus be determined by ψ as

isogeometric B-spline basis functions based on Eq. (21), are plotted in Fig. 1(b). As shown in Fig. 2(a), by inserting 3/8 and 5/8 into the initial knot vector, the new nodes set would become ξ [1] = {0, 1 8, 5 16, 7 16, 9 16, 11 16, 7 8, 1} through Eq. (19). It could be found that new linear reproducing points are added into the refined elements where the nodes within the rest domain remain constant.

ε̇ = μ

∂ψ (σ ) , ∂σ

(32)

where µ is a non-negative plastic multiplier. According to the virtual work theorem, the limit load factor λ for a 5

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 7. (a) Full mode and (b) its quarter mode with the symmetric boundary conditions of a grooved rectangular plate.

Fig. 8. (a) Initial mesh and (b) mesh after 3 steps of refinement for a grooved rectangular plate.

Table 3 Limit load factors of a grooved rectangular plate. Methods Prager and Hodge [48]

Yan [49]

Casciaro and Cascini [50]

Nguyen-Xuan et al. [18]

Present method

0.5000

0.5000–0.5770

0.5680

0.5593

0.5649

specific body is defined as the ratio between the external and internal virtual work. By further defining a space L = {u̇ ∈ V |Win (u̇) = 1} , λ is obtained by solving an optimization problem as shown in [35]:

λp = max{∃ σ ∈ S |Win (σ , u̇) = λWex (u̇), ∀ u̇ ∈ V } = min D (u̇), u̇ ∈ L

(33) where D (u̇) = max Win (σ , u̇) and σ are the admissible stresses on the σ∈S

Fig. 9. Plastic dissipation distribution for a grooved plate.

convex yield surface. Since the kinematic formulation is considered in this work, Eq. (33) will be used to evaluate the upper-bound limit load factors λp of structures.

n

n

uh (ξ ) = ∑I =S 1 ΦI (ξ ) q̇ I + ∑J H= 1 ΦJ (ξ )[H (ξ ) − H (ξ J )] cJ̇ n

{

4

β

}

+ ∑KK= 1 ΦK (ξ ) ∑β = 1 [Q (ξ ) − Q (ξ )] ḃ K ,

4.2. Enrichment functions for cracks

(34)

where q̇ I = [uİ , vİ ]T denotes the velocities of the nodal displaceβ ments associated with the standard AIMF coupling nodes; ċJ and ḃ K

Similar to the enrichment functions adopted in XFEM and XIGA, the velocity field of the cracked structures can be expressed as 6

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 10. (a) Top view and (b) front view of a rectangular thin plate subject to simply supported boundary conditions.

Fig. 11. (a) Initial mesh and (b) mesh after 4 steps of refinement for a square thin plate with l1/l2 = 1. Table 4 Limit load factors λp of a square thin plate. Methods Le et al. [52]

Le et al. [53]

Zhou et al. [54]

Le et al. [55]

Nguyen-Xuan et al. [56]

Present method

24.98

25.01

25.07

24.93

25.023

24.932

H (ξ ) =

{+− 11

if (ξ − ξ ∗)·n > 0 , otherwise

(35)

where ξ ∗ is the projection of node ξ on the crack and n represents the normal vector. The crack tip enrichment functions Q are defined in polar coordinates (r , θ) as

Q (r , θ) θ = ⎧ r sin ⎛ ⎞, ⎨ ⎝2⎠ ⎩

Fig. 12. Plastic dissipation distribution for a square thin plate with l1/l2 = 1.

denote those associated with the crack surfaces and tips, respectively. H in the above equation denotes the Heaviside functions that are defined by

θ r cos ⎛ ⎞, ⎝2⎠

θ r sin ⎛ ⎞ sin(θ), ⎝2⎠

θ r cos ⎛ ⎞ sin(θ)⎫ . ⎬ ⎝2⎠ ⎭ (36)

The compatible strain rates can be expressed by nodal velocities including both standard and enriched velocity unknows ḋ I as 7

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 13. (a) Initial mesh and (b) mesh after 4 steps of refinement for a rectangular thin plate with l1/l2 = 2.

BS and BT indicate the part of matrix B for the crack surface and tip which can be respectively calculated by ⎡ ΦI , x HI BSI = ⎢ ⎢ ⎣ ΦI , y HI ⎡ ΦI , x QI =⎢ ⎢ ⎣ ΦI , y QI

+ ΦI HI , x 0 ⎤ 0 ΦI , y HI + ΦI HI , y ⎥ and BTI ⎥ + ΦI HI , y ΦI , x HI + ΦI HI , x ⎦ + ΦI QI , x 0 ⎤ 0 ΦI , y QI + ΦI QI , y ⎥. ⎥ + ΦI QI , y ΦI , x QI + ΦI QI , x ⎦

(40)

The plastic dissipation of the body made of a rigid-perfectly material can be rewritten by

D h (u̇ h) = σ0

Fig. 14. Plastic dissipation distribution for a rectangular thin plate with l1/ l2 = 2.

ε̇ =

∑ BI ḋ I where the strain matrix B is given by

BI = [BNI

BSI

BTI ].

(38)

I=1

∫Ω

ε Ṫ Θε ̇ dΩ,

I

(41)

ng

D h (u̇ h) ≃ σ0

BN

in the Eq. (38) represents the standard part of matrix B and is defined by

⎡ ΦI ,1 0 ⎤ BNI = ⎢ 0 ΦI ,2 ⎥. ⎢Φ ⎥ ⎣ I ,2 ΦI ,1 ⎦

nt

ε Ṫ Θε ̇ dΩ = σ0 ∑

where nt denotes the total number of AIMF coupling nodes. Since the Gaussian integral is applied, the strain rate ε ̇ is evaluated across the Gauss points over the problem domain. Eq. (41) can be further reformulated as

(37)

I

∫Ω

∑ w¯ n |Jn|

εṅ T Θε ṅ ,

n=1

(42)

where ng is the total number of Gauss points; w¯ n and |Jn| is the weight and determinant of the Jacobian matrix associated with Gauss point n. The optimization problem shown in the Eq. (31) can be rewritten as

(39)

Fig. 15. (a) Top view and (b) front view of a circular thin plate subject to the simply supported boundary conditions. 8

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 16. (a) Initial mesh and (b) mesh after 3 steps of refinement for a circular thin plate. ng

D h (u̇ h) ≃ σ0

∑ w¯n

|Jn | ∥GTεṅ ∥,

n=1

(44)

where the matrix G is the Cholesky factor of Θ and can be defined as

G=

2 1 ⎡ ⎢1 3⎢ ⎣0

0 0⎤ 3 0 ⎥. 0 1⎥ ⎦

(45)

A vector of additional variables ρn is introduced to reduce the above equation and is given as

ρ ⎡ 1⎤ ρn = ⎢ ρ2 ⎥ = GTεṅ . ⎢ ⎣ ρ3 ⎥ ⎦

(46)

In this way, Eq. (43) can finally be reformulated into a problem of minimizing a sum of norms as

Fig. 17. Plastic dissipation distribution for a circular thin plate.

n

λ p = min {σ0 ∑i =g 1 w¯ i |Ji | vi} , Wex (u̇ ) = 1 ⎧ ⎪ u̇ = u̇ on Γu ⎪ s.t. ρ = CTε ̇ (i = 1, 2, ⋯n ) i g ⎨ i ⎪ T ⎪ ρi ρi ⩽ vi (i = 1, 2, ⋯ng ) ⎩

where λp is the limit load factor, and the first constraint in the Eq. (47) indicates the inequality conditions. The total number of variables for this optimization problem is nvar = ndofs + 4 ng, where ndofs is the total number of degree of freedoms for the entire domain. The Mosek optimization package is employed as an effective tool to solve Eq. (47) [18].

Fig. 18. A central cracked rectangular plate subject to uniaxial tension. n

λ p = min {σ0 ∑i =g 1 w¯ i |Ji |

5. Numerical results

εi̇ T Θε i̇ } .

W (u̇ ) = 1 on Γu s.t.⎧ ex ⎨ ⎩ u̇ = u̇

(47)

In this section, we shall evaluate the performance of the present method through a series of benchmark problems involving elastic-perfectly-plastic materials under the plane stress condition, where the von Mises yield criterion is adopted.

(43)

4.3. Solution procedure using second-order cone programming

5.1. Uncracked structures

Andersen et al. [36] proved that the optimization problem shown in the Eq. (43) could be further simplified into a problem of minimizing a sum of norms. Such a problem can be solved by a SOCP as

The limit analysis problems of uncracked structures are first studied to validate the present approach, where different material shapes and loading conditions are considered. 2 × 2 Gauss points are adopted for integration to minimize computational cost. 9

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 19. (a) Initial mesh and (b) mesh after 3 steps of refinement for a central cracked rectangular plate with lc/l = 0.5. Table 5 Limit load factors λp of a central cracked rectangular plate with various crack lengths lc. lc/l

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Methods Miller [58]

Tran et al. [59]

Nguyen-Xuan et al. [18]

Present method

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

0.9452 0.8426 0.7443 0.6398 0.5419 0.4356 0.3381 0.2292

0.9183 0.8239 0.7293 0.6303 0.5299 0.4295 0.3284 0.2266

0.8985 0.8105 0.7074 0.6098 0.5101 0.4099 0.3089 0.2080

Fig. 20. Plastic dissipation distribution for a central cracked plate with lc/ l = 0.5.

Fig. 22. The model of an edge cracked rectangular plate under uniaxial tension load.

with a uniformly coarse mesh as shown in Fig. 4(a). Fig. 4(b) displays the final mesh domain after 3 steps of refinement. It is evident that the elements associated with higher plastic deformations are automatically identified and refined. Gaydon and McCrum [37] provided the analytical solution to be λp = 0.8σ1/σ0 based on the von Mises yield criterion. Fig. 5 compares the convergence rates of limit load factors obtained by the present approach with the reference solutions: quadratic IGA and isogeometric-meshfree (IMF) coupling approach with global refinement. As expected, an AIMF coupling approach with global refinement achieves almost the same convergence rate of limit load factor as that obtained by quadratic IGA. After the adaptive refinement strategy is applied, the present approach significantly improves the convergence order to over 1.5. As listed in Table 1, the results obtained through the proposed approach are found to be in good agreement with both the analytical solutions and reference numerical results. It is obvious that the present method not only achieves the highest convergence rate but also offers the most accurate results with the least optimization variables. At the collaspe state of the plate, Fig. 6 demonstrates that the plastic dissipation distribution is captured with a high resolution. The geometric effect of the size of the central hole on the limit load factor is

Fig. 21. The convergence of limit load factor solutions for a central cracked plate with lc/l = 0.5 obtained by different numerical methods.

5.1.1. Square plate with a central hole As shown in Fig. 3(a), a square plate with a central circular hole under biaxial loads σ1 and σ2 is considered in this case. This is a classic benchmark problem which helps to validate the present approach. The length of the square is denoted by l while the radius of the central hole is indicated by r. Fig. 3(b) illustrates that this model can be further simplified into a quarter model based on the symmetry boundary conditions. By setting σ2 to 0, the adaptive refinement procedure starts 10

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 23. (a) Initial mesh and (b) mesh after 3 steps of refinement for an edge cracked rectangular plate with lc/l = 0.06.

Fig. 24. (a) Initial mesh and (b) mesh after 3 steps of refinement for an edge cracked rectangular plate with lc/l = 0.3.

investigated through various ratios of radius r to length l. The simulation results obtained by the present approach are summarized in Table 2. As expected, the present solutions show high reliability as the results are in good accordance with reference solutions.

∇2 = [∂2 ∂x 2 ∂2 ∂y 2 2∂2 ∂xy]T . By substituting the strain rates ε ̇ = zκ ̇ into Eq.(23), the plastic dissipation Dh for a thin plate can be given by

5.1.2. Grooved rectangular plate As shown in Fig. 7, a grooved rectangular plate under the uniaxial in-plane tensive load σ0 is considered. The dimensions of the plate are 2 l × l, while the groove on each side of the plate has a radius of r = 0.25 l. A half model with the symmetry boundary conditions is studied for this limit analysis problem. For this limit analysis problem, Prager and Hodge [46] determined the analytical solution of the limit load factor to be 0.5 based on the Tresca yield criterion, while Yan [47] obtained a range of 0.5 to 0.577 based on the von Mises yield criterion. Fig. 8 shows that the elements along the central line experience the highest strain rates at the collapse state. As presented in Table 3, it is observed that the result given by the present approach not only agrees well with reference solutions but also falls within the reliability interval of the analytical approach presented by Yan [47]. Fig. 9 plots plastic dissipation distribution across the grooved plate at the collapse state.

Consequently, the optimization problem shown in Eq. (47) can be reformulated as

D h (κ ̇) = σ0

σ0 h2 4

∫Ω

κ ̇TΘκ ̇ dΩ.

(50)

n

Wex (κ ̇) = 1 ⎧ ⎪ ̇ w = ẇ on Γw ⎪ s.t. ρ = CTB ẇ (i = 1, 2, ⋯n ) p g i i ⎨ i ⎪ T ⎪ ρi ρi ⩽ vi (i = 1, 2, ⋯ng ) ⎩

(51)

where and the matrix B for thin plate is defined as B = [Φi, xx Φi, yy 2Φi, xy ]T , and σp denotes the yield stress. As shown in Fig. 10, a thin plate with each of its side subject to simply supported boundary conditions is investigated. The plate is assumed to have a yield stress σ0 = 250 MPa and a thickness h = 0.1 l1. A transverse load σp = 1 MPa is uniformly applied on the plate along the z-direction. By setting mp to be 0.25σ0h2, the limit load factor λp for thin plate could be given by λp = σpl1l2/mp. Firstly, the plate is assumed to be a square with each side measuring 10 m. The adaptive refinement procedure starts with a coarse mesh as depicted in Fig. 11(a). The final refinement step depicted in Fig. 11(b) shows that the regions of the plate associated with high plastic dissipated power exist along the diagonal lines of the plate. Fig. 12 further demonstrates that the plate is most likely to yield diagonally at its plastic collapse state. Although analytical solutions have not been proposed for this limit analysis problem, the simulation results agree well with reference solutions, as presented in Table 4. Next, a rectangular thin plate with a length-width-ratio l1/l2 = 2 is investigated for its limit load factor. The adaptive refinement procedure of this plate is demonstrated in Fig. 13. It could be determined that the

(48)

The curvatures field κ can thus be expressed in terms of w as

κ = [ κ x κ y κ xy ]T = [∂2w ∂x 2 ∂2w ∂y 2 2∂2w ∂xy]T .

ε Ṫ Θε ̇ dΩdz =

λ p = min {0.25σ0 h2 ∑i =g 1 w¯ i |Ji | vi} ,

5.1.3. Thin plate The application of the proposed approach is extended to the limit analysis of a thin plate of which the mid-plane is denoted as Ω = 2 . By adopting the rotation-free plate elements [51], the transversal displacement w is regarded as the only independent variable. To satisfy the Kirchhoff’s thin plate conditions, the rest two displacement components are given by [51]

θx = ∂w ∂x and θy = ∂w ∂y.

h 2

∫−h 2 ∫Ω

(49)

ẇ is subsequently introduced to denote the velocity of w. The curvature rates κ̇ can be subsequently defined by ∇2 ẇ , where 11

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 26. A grooved cracked plate subject to uniaxial tension.

and Silva [57] to be 6.561. The simulation result obtained by the present method is found to be 6.569 which is in a good agreement with the reference solution.

5.2. Cracked structures After the robustness of the present approach has been verified for limit analysis of uncracked structures and plates, its application is extended to cracked structures under a uniaxial tension load σ0. The limit load factor λp for cracked structures is defined as σp/σ0, where σp is the plastic collapse load of the corresponding material. It should be noted that 2 × 2 Gauss points are applied for standard elements, 4 × 4 Gauss points for the elements along the crack surface, and 6 × 6 Gauss points for the elements around the crack-tip area for the following cases.

5.2.1. Central cracked plate As shown in Fig. 18, a central cracked 2 l × l rectangular plate is considered, with the crack length assumed to be lc. A pair of uniaxial tension loads σ0 = σp are uniformly applied along two short ends of the plate. Firstly, the plate with a crack length ratio lc/l of 0.5 is studied to examine the convergence of the AIMF coupling approach. The final steps of adaptive refinement are depicted in Fig. 19, which demonstrates that the elements within the regions of the material containing high plastic dissipation are automatically identified and refined along lines extended from the crack tips. Fig. 20 further demonstrates this phenomenon by plotting the plastic dissipation distribution over the problem domain at collapse state. Miller [58] developed an analytical solution of this limit analysis problem based on the von Mises yield criterion where λp = 1 − lc/l. By setting the crack length ratio lc/l = 0.5, the limit load factor can thus be calculated to be 0.5. The convergence rate achieved by the present method is compared with other available reference results in Fig. 21. Even if the proposed approach adopts the coarse mesh with fewer variables than those employed by the XFEM and quadratic XIGA, the present approach is still able to achieve the most accurate simulation results with the highest convergence rate. The effect of the crack geometry on the limit load factor was further investigated via various crack length ratios, which is tabulated in Table 5. Limit load factors achieved by the present approach are in the best agreement with analytical solutions compared with other numerical methods, which demonstrates the robustness of the proposed approach for limit analysis problems of cracked structures.

Fig. 25. Plastic dissipation distribution for edge cracked plates with (a) lc/ l = 0.06 and (b) lc/l = 0.3.

pattern of the area containing high plastic dissipation within the rectangular thin plate is similar to that possessed by the square plate. Fig. 14 further shows that the rectangular plate is also most likely to suffer plastic collapse along its diagonal lines. The limit load factor of a circular thin plate with a radius r and a thickness h = 0.02r is further studied to demonstrate the effectiveness of the present method on curved structures. As shown in Fig. 15, the outer edge of the circular plate is subject to simply supported boundary conditions while a transverse load σp = 1 MPa is uniformly applied across the plate. The initial mesh and final mesh after 3 refinement steps are depicted in Fig. 16. The plastic dissipation energy of the circular plate is further detailed in Fig. 17. As expected, both figures indicate that the plastic collapse is most likely to take place within the central area of the circular plate. By adopting the von Mises yield criterion, the reference solution for this problem was proposed by Capsoni

Table 6 Limit load factors λp of an edge cracked rectangular plate with various crack lengths. Methods

Nguyen-Xuan et al. [18] Present method

Crack length ratio lc/l 0.06

0.1

0.14

0.2

0.3

0.4

0.5

0.6

0.9364 0.9397

0.8816 0.8895

0.8327 0.8369

0.7197 0.7249

0.5298 0.5346

0.3662 0.3713

0.2339 0.2366

0.1352 0.1379

12

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

Fig. 27. (a) Initial mesh and (b) mesh after 3 steps of refinement for a grooved cracked rectangular plate with lc/l = 0.3.

5.2.3. Grooved cracked plate As depicted in Fig. 26, the last case investigates the plastic collapse behaviours of a grooved cracked plate under uniaxial tension σ. The length of the plate is twice its width L, and the radius R of the groove equals 0.25 l. The crack length ratio is defined by lc/(1 − 2r) for this cracked structure. The adaptive refinement process is plotted in Fig. 27, where high strain rates are found to exist around the crack tip. As depicted in Fig. 28, the plastic dissipation distribution of the grooved cracked plate at its collapse state is found to align along the direction of the crack path. As reported in Table 7, the limit load factors achieved by the present approach are in good agreement with the XIGA solutions. The crack lengths lc were also found to have a noticeable influence on the limit load factor of cracked structures as a larger discontinued domain would further weaken the elastic-perfectly-plastic materials. 6. Conclusions Fig. 28. Plastic dissipation distribution for a grooved cracked plate lc/l = 0.3.

In this paper, an AIMF coupling approach is combined with an indicator of plastic dissipation to perform limit analysis on cracked structures. In the present approach, IGA and meshfree basis functions are unified through an equivalence based on the reproducing conditions, thus the local refinement of IGA can be performed in a direct meshfree manner across the entire domain. In this way, the present approach inherits advantages from both IGA and meshfree approximation, including an exact geometric representation and flexibility of adaptive refinement. An effective indicator of the plastic dissipation has been further developed to guide the adaptive local refinement process for limit analysis problems. Only the region of a cracked structure containing sufficiently large dissipated work would be automatically captured and refined. Consequently, the computational cost incurred by the excessive overhead of IGA control points with increasing levels of refinement can be significantly reduced for limit analysis problems [18]. The robustness and effectiveness of the present method have been demonstrated through the investigations on a series of benchmark limit analysis problems. The application of the AIMF coupling approach on the limit and shakedown analysis of shell structures will be explored in future work.

Table 7 Limit load factors λp of a grooved cracked rectangular plate. Methods

Nguyen-Xuan et al. [18] Present method

Crack length ratio lc/(1 − 2r) 0.1

0.2

0.3

0.4

0.5

0.4787 0.4771

0.3746 0.3730

0.2865 0.2798

0.2006 0.1980

0.1395 0.1313

5.2.2. Edge cracked plate As shown in Fig. 22, the 2 l × l rectangular plate under uniaxial tension load σ0 now contains an edge crack with a length of lc. For this benchmark limit analysis problem, Ewing and Richards [60] proposed an analytical solution by employing the von Mises yield criterion. The plastic collapse behaviours of edge cracked plate would change once the crack length ratio lc/l exceeds 0.146. Therefore, the adaptive procedures for the edge cracked rectangular plates with the lc/l of 0.06 and 0.3 are plotted in Figs. 23 and 24, respectively. Obviously, the crack length affects the domains associated with large plastic dissipation. When the crack length ratio is below 0.146, the plastic collapse of the plate is most likely to occur at the crack tip and along the two directions symmetric to the crack surface. On the other hand, the edge cracked plate tends to yield along the crack surface when the crack length ratio lc/l exceeds 0.146. This phenomenon could be better explained by the plastic dissipation distribution over the problem domain as presented in Fig. 25. The limit load factors corresponding to various crack length ratios lc/l are presented in Table 6. It is found that the limit load factor is inversely proportional to the crack length lc. The simulation results match well with both the analytical and reference solutions as anticipated.

Acknowledgements This research work was conducted in the SMRT-NTU Smart Urban Rail Corporate Laboratory with funding support from the National Research Foundation (NRF), Singapore, SMRT, Singapore and Nanyang Technological University, Singapore. References [1] K. Krabbenhoft, L. Damkilde, A general non-linear optimization algorithm for lower bound limit analysis, Int. J. Numer. Meth. Eng. 56 (2) (2003) 165–184.

13

Theoretical and Applied Fracture Mechanics 106 (2020) 102426

J. Huang, et al.

[2] J.J. Muñoz, A. Huerta, J. Bonet, J. Peraire, A note on upper bound formulations in limit analysis, Int. J. Numer. Meth. Eng. 91 (8) (2012) 896–908. [3] M. Vicente da Silva, A.N. Antão, A non-linear programming method approach for upper bound limit analysis, Int. J. Numer. Meth. Eng. 72 (10) (2007) 1192–1218. [4] R. Hill, On discontinuous plastic states, with special reference to localized necking in thin sheets, J. Mech. Phys. Solids 1 (1) (1952) 19–30. [5] C.V. Le, H. Nguyen-Xuan, H. Askes, S.P.A. Bordas, T. Rabczuk, H. Nguyen-Vinh, A cell-based smoothed finite element method for kinematic limit analysis 83 (12) (2010) 1651–1674. [6] H. Li, Y. Liu, X. Feng, Z. Cen. Limit analysis of ductile composites based on homogenization theory, Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 2003. 459(2031): pp. 659-675. [7] H. Li, Y. Liu, X. Feng, Z. Cen, Upper bound analysis of plastic limit loads forcomposite materials and structures, Acta Metallurgica Sinica(English letters) 19 (1) (2002) 85–89. [8] M. Surendran, S. Natarajan, G.S. Palani, S.P.A. Bordas, Linear smoothed extended finite element method for fatigue crack growth simulations, Eng. Fract. Mech. 206 (2019) 551–564. [9] C. Jansari, S. Natarajan, L. Beex, K. Kannan, Adaptive smoothed stable extended finite element method for weak discontinuities for finite elasticity, Eur. J. Mech. A. Solids 78 (2019) 103824. [10] A. Makrodimopoulos, C.M. Martin, Lower bound limit analysis of cohesive-frictional materials using second-order cone programming, Int. J. Numer. Meth. Eng. 66 (4) (2006) 604–634. [11] E.D. Andersen, C. Roos, T. Terlaky, On implementing a primal-dual interior-point method for conic quadratic optimization, Math. Program. 95 (2) (2003) 249–277. [12] T.D. Tran, C.V. Le, Extended finite element method for plastic limit load computation of cracked structures 104 (1) (2015) 2–17. [13] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Eng. 194 (39) (2005) 4135–4195. [14] J.A. Cottrell, T.J.R. Hughes, A. Reali, Studies of refinement and continuity in isogeometric structural analysis, Comput. Methods Appl. Mech. Eng. 196 (41) (2007) 4160–4183. [15] S.S. Ghorashi, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, Int. J. Numer. Meth. Eng. 89 (9) (2011) 1069–1101. [16] G. Bhardwaj, S.K. Singh, I.V. Singh, B.K. Mishra, T. Rabczuk, Fatigue crack growth analysis of an interfacial crack in heterogeneous materials using homogenized xiga, Theor. Appl. Fract. Mech. 85 (2016) 294–319. [17] S.K. Singh, I.V. Singh, B.K. Mishra, G. Bhardwaj, T.Q. Bui, A simple, efficient and accurate bézier extraction based t-spline xiga for crack simulations, Theor. Appl. Fract. Mech. 88 (2017) 74–96. [18] H. Nguyen-Xuan, L.V. Tran, C.H. Thai, C.V. Le, Plastic collapse analysis of cracked structures using extended isogeometric elements and second-order cone programming, Theor. Appl. Fract. Mech. 72 (2014) 13–27. [19] H. Nguyen-Xuan, C.T. Wu, G.R. Liu, An adaptive selective es-fem for plastic collapse analysis, Eur. J. Mech. A. Solids 58 (2016) 278–290. [20] P. Hennig, M. Kästner, P. Morgenstern, D. Peterseim, Adaptive mesh refinement strategies in isogeometric analysis— a computational comparison, Comput. Methods Appl. Mech. Eng. 316 (2017) 424–448. [21] F. Liaghat, M.R. Hematiyan, A. Khosravifard, T. Rabczuk, A robust meshfree method for analysis of cohesive crack propagation problems, Theor. Appl. Fract. Mech. 104 (2019) 102328. [22] W. Ma, G. Liu, H. Ma, A smoothed enriched meshfree galerkin method with twolevel nesting triangular sub-domains for stress intensity factors at crack tips, Theor. Appl. Fract. Mech. 101 (2019) 279–293. [23] C. Liu, C. Ku, J. Xiao, W. Yeih, A novel spacetime collocation meshless method for solving two-dimensional backward heat conduction problems, Computer Modeling in Engineering and Sciences 118 (1) (2019) 229–252. [24] G. Kakuba, J.M. Mango, M.J.H. Anthonissen, Convergence properties of local defect correction algorithm for the boundary element method, Computer Modeling in Engineering and Sciences 119 (1) (2019) 207–225. [25] D. Wang, H. Zhang, A consistently coupled isogeometric–meshfree method, Comput. Methods Appl. Mech. Eng. 268 (2014) 843–870. [26] H. Zhang, D. Wang, Reproducing kernel formulation of b-spline and nurbs basis functions: A meshfree local refinement strategy for isogeometric analysis, Comput. Methods Appl. Mech. Eng. 320 (2017) 474–508. [27] C.V. Le, H. Nguyen-Xuan, H. Askes, S.P.A. Bordas, T. Rabczuk, H. Nguyen-Vinh, A cell-based smoothed finite element method for kinematic limit analysis, Int. J. Numer. Meth. Eng. 83 (12) (2010) 1651–1674. [28] A. Makrodimopoulos, C.M. Martin, Upper bound limit analysis using simplex strain elements and second-order cone programming, Int. J. Numer. Meth. Eng. 31 (6) (2007) 835–865. [29] H. Nguyen-Xuan, G.R. Liu, An edge-based finite element method (es-fem) with

[30]

[31]

[32] [33] [34] [35]

[36] [37]

[38] [39] [40]

[41] [42]

[43]

[44] [45] [46] [47]

[48] [49] [50] [51] [52] [53]

[54] [55]

[56]

[57]

[58] [59] [60]

14

adaptive scaled-bubble functions for plane strain limit analysis, Comput. Methods Appl. Mech. Eng. 285 (2015) 877–905. K.M. Liew, C. Feng, Y. Cheng, S. Kitipornchai, Complex variable moving leastsquares method: A meshless approximation technique, Int. J. Numer. Meth. Eng. 70 (1) (2007) 46–70. N. Sukumar, R.W. Wright, Overview and construction of meshfree basis functions: From moving least squares to entropy approximants, Int. J. Numer. Meth. Eng. 70 (2) (2007) 181–205. G.R. Liu, Y.T. Gu, An introduction to meshfree methods and their programming, Springer Publishing Company, Incorporated, 2010, p. 479. A. Capsoni, L. Corradi, A finite element formulation of the rigid–plastic limit analysis problem, Int. J. Numer. Meth. Eng. 40 (11) (1998) 2063–2086. E. Christiansen, K.D. Andersen, Computation of collapse states with von mises type yield condition, Int. J. Numer. Meth. Eng. 46 (8) (1999) 1185–1202. X.-L. Yang, Upper bound limit analysis of active earth pressure with different fracture surface and nonlinear yield criterion, Theor. Appl. Fract. Mech. 47 (1) (2007) 46–56. K. Andersen, E. Christiansen, M. Overton, Computing limit loads by minimizing a sum of norms, SIAM Journal on Scientific Computing 19 (3) (1998) 1046–1062. F.A. Gaydon, A.W. McCrum, A theoretical investigation of the yield point loading of a square plate with a central circular hole, J. Mech. Phys. Solids 2 (3) (1954) 156–169. J. Groβ-Weege, On the numerical assessment of the safety factor of elastic-plastic structures under variable loading, Int. J. Mech. Sci. 39 (4) (1997) 417–433. S. Chen, Y. Liu, Z. Cen, Lower-bound limit analysis by using the efg method and non-linear programming, Int. J. Numer. Meth. Eng. 74 (3) (2008) 391–415. N. Zouain, L. Borges, J. L. s. Silveira. An algorithm for shakedown analysis with nonlinear yield functions, Computer Methods in Applied Mechanics and Engineering, 2002. 191(23-24): pp. 2463-2481. X. Zhang, Y. Liu, Z. Cen, Boundary element methods for lower bound limit and shakedown analysis, Eng. Anal. Boundary Elem. 28 (8) (2004) 905–917. T.N. Tran, G.R. Liu, H. Nguyen-Xuan, T. Nguyen-Thoi, An edge-based smoothed finite element method for primal–dual shakedown analysis of structures, Int. J. Numer. Meth. Eng. 82 (7) (2010) 917–938. H. Nguyen-Xuan, T. Rabczuk, T. Nguyen-Thoi, T.N. Tran, N. Nguyen-Thanh, Computation of limit and shakedown loads using a node-based smoothed finite element method, Int. J. Numer. Meth. Eng. 90 (3) (2012) 287–310. M. Heitzer. Traglast-und einspielanalyse zur bewertung der sicherheit passiver komponenten, Berichte-forschungszentrum Julich Jul, 1999. 1(3704): pp. ALL-ALL. T.N. Trần, Limit and shakedown analysis of plates and shells including uncertainties (2008). W. Prager, P.G. Hodge, Theory of perfectly plastic solids, Wiley, New York, 1951. A.-M. Yan, L. Université de, a. Faculté des sciences. Contributions to the direct limit state analysis of plastified and cracked structures. Liege, Belgique: Universite de Liege, Faculte des Sciences appliquees; 1999. W. Prager, P.G. Hodge, Theory of perfectly plastic solids, Dover Publications, 1968. A.-M. Yan, Contributions to the direct limit state analysis of plastified and cracked structures (1999). R. Casciaro, L. Cascini, A mixed formulation and mixed finite elements for limit analysis, Int. J. Numer. Meth. Eng. 18 (2) (1982) 211–243. E. Oñate, F. Zárate, Rotation-free triangular plate and shell elements, Int. J. Numer. Meth. Eng. 47 (1–3) (2000) 557–603. C.V. Le, M. Gilbert, H. Askes, Limit analysis of plates and slabs using a meshless equilibrium formulation, Int. J. Numer. Meth. Eng. 83 (13) (2010) 1739–1758. C.V. Le, M. Gilbert, H. Askes, Limit analysis of plates using the efg method and second-order cone programming, Int. J. Numer. Meth. Eng. 78 (13) (2009) 1532–1552. S. Zhou, Y. Liu, S. Chen, Upper bound limit analysis of plates utilizing the c1 natural element method, Comput. Mech. 50 (5) (2012) 543–561. C.V. Le, H. Nguyen-Xuan, H. Nguyen-Dang, Upper and lower bound limit analysis of plates using fem and second-order cone programming, Comput. Struct. 88 (1) (2010) 65–73. H. Nguyen-Xuan, C.H. Thai, J. Bleyer, P.V. Nguyen, Upper bound limit analysis of plates using a rotation-free isogeometric approach, Asia Pacific Journal on Computational Engineering 1 (1) (2014) 12. A. Capsoni, M. Silva, A finite element formulation of mindlin plates for limit analysis, International Journal for Numerical Methods in Biomedical Engineering 27 (2011) 143–156. A.G. Miller, Review of limit loads of structures containing defects, Int. J. Press. Vessels Pip. 32 (1) (1988) 197–327. T.D. Tran, C.V. Le, Extended finite element method for plastic limit load computation of cracked structures, Int. J. Numer. Meth. Eng. 104 (1) (2015) 2–17. D.J.F. Ewing, C.E. Richards, The yield-point loads of singly-notched pin-loaded tensile strips, J. Mech. Phys. Solids 22 (1) (1974) 27–36.