Displacement and equilibrium mesh-free formulation based on integrated radial basis functions for dual yield design

Engineering Analysis with Boundary Elements 71 (2016) 92–100

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Displacement and equilibrium mesh-free formulation based on integrated radial basis functions for dual yield design Phuc L.H. Ho a,b, Canh V. Le c,n, T. Tran-Cong d a

Faculty of Applied Mechanics and Civil Engineering, University of Technical Education - HCMC, Viet Nam Computational Engineering Laboratory, Institute for Computational Science and Technology, Ho Chi Minh City, Vietnam c Department of Civil Engineering, International University - VNU HCM, Viet Nam d Computational Engineering and Science Research Centre, University of Southern Queensland, Toowoomba, QLD 4350, Australia b

art ic l e i nf o

a b s t r a c t

Article history: Received 19 December 2015 Received in revised form 5 May 2016 Accepted 18 July 2016

This paper presents displacement and equilibrium mesh-free formulation based on integrated radial basis functions (iRBF) for upper and lower bound yield design problems. In these approaches, displacement and stress fields are approximated by the integrated radial basis functions, and the equilibrium equations and boundary conditions are imposed directly at the collocation points. In this paper it has been shown that direct nodal integration of the iRBF approximation can prevent volumetric locking in the kinematic formulation, and instability problems can also be avoided. Moreover, with the use of the collocation method in the static problem, equilibrium equations and yield conditions only need to be enforced at the nodes, leading to the reduction in computational cost. The mean value of the approximated upper and lower bound is found to be in excellent agreement with the available analytical solution, and can be considered as the actual collapse load multiplier for most practical engineering problems, for which exact solution is unknown. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Limit analysis Yield design Integrated radial basis function Mesh-free method Second order cone programming

1. Introduction The estimation of the plastic limit load of engineering structures has been of great interest to practical engineers. Elastic– plastic incremental analysis method can be employed to obtain such a limit load, but the plastic yield design (limit analysis) method has been viewed to be more effective [1,2], i.e. without performing a cumbersome series of incremental elastic–plastic analysis, and only requiring the knowledge of a local strength criterion. Based on bound theorems and numerical discretization techniques, the yield design approach can directly find the critical status of a structure or body under static loading and then the collapse load can be determined without intermediate steps. A lower bound on the actual limit load of a structure or body can be achieved by using the static theorem and approximated stress fields, while the upper bound is obtained as a result of combining displacement-based model and kinematic theorem [11]. In the static yield design formulation, the assumed stress fields are often expressed in terms of nodal stress values. In the framework of equilibrium finite elements, these approximated fields are also required to satisfy a priori equilibrium conditions n

Corresponding author. E-mail address: [email protected] (C.V. Le).

http://dx.doi.org/10.1016/j.enganabound.2016.07.010 0955-7997/& 2016 Elsevier Ltd. All rights reserved.

within elements and at their interfaces [3–6,11]. Due to these additional conditions, construction of such fields is often difficult. Compared with the equilibrium models, the displacement formulation is more popular. This may be because of the facts that the internal compatibility condition can be satisfied straightaway in the assembly scheme, and that essential (kinematic) boundary conditions can be enforced directly. In recent decades, the development and application of meshfree methods have attracted much attention due to their flexibility, e.g. requiring nodal data only and no need of nodal connectivity. The EFG method [7], one of the most widely used mesh-free methods, has been applied successfully to the framework of yield design problems [8–12], showing that the method is, in general, well suited for yield design problems and that accurate solutions can be obtained with a minimal computational cost. However, a typical limitation of the EFG method is that its shape functions do not hold Kronecker delta property, leading to difficulty in enforcing essential boundary conditions. Mesh-free methods using radial basis functions (RBFs) have also been developed in parallel [13–19]. These RBF-based methods have been applied successfully to a wide range of computational problems [20–26]. It has been shown that the RBF approximation is powerful to represent regular and smooth functions in arbitrary geometries and high-dimensional space, and enjoys spectral accuracy and exponential

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

convergence [17]. The aim of this paper is to study the performance of the integrated radial basis function-based mesh-free method in the framework of yield design problems. The indirect/integrated radial basis function (iRBF) approaches proposed in [28,29] are employed to approximate both displacement and stress fields. Multi-quadric iRBF method generally results in a high order approximation of the displacement fields (here, for convenience, ‘displacement rate or velocity’ is termed ‘displacement’), and hence volumetric locking phenomena in the kinematic yield design formulation can be prevented. Moreover, the stress fields constructed based on iRBF are smooth over the entire problem domain, and consequently there is no need to enforce continuity conditions at interfaces within the problem domain. With the use of iRBF-approximated stress fields the strong-form of equilibrium equations can be satisfied in a point-wise manner using a collocation method. In addition, the iRBF-based approximation possesses the Kronecker delta property. As a result, kinematic and static boundary conditions can be imposed as easily as in the finite element method (FEM). Finally, the kinematic and static formulations based on iRBF discretization are formulated as a conic optimization problem, ensuring that they can be solved using available efficient solvers.

2. Fundamentals of yield design theory In this section, fundamentals of yield design theory are recalled, more details can be found in [27,45,51-54]. Consider a elastic–plastic body of area Ω ∈ 2 with fixed boundary Γu and free portion Γt, satisfying Γu ∪ Γt = Γ , Γu ∩ Γt = ⊘, and is subjected to body forces f and surface tractions t . Let Σ denotes a space of a statically admissible stress state, whereas @ is a space of a kinematically admissible displacement state. For smooth fields σ and u , the classical form of the equilibrium equation can always be transformed to a more precise variational form as

a (σ , u) = F (u),

∀u∈@

∫Ω σ T ϵ (u) dΩ

(2)

∫Ω f T u dΩ + ∫Γ

F (u) =

tT u dΓ

t

(4)

(5)

where the set ℓ is defined by ℓ = {u ∈ Y |F (u) = 1}, and the plastic dissipation D (u) is defined by

D (u) = maxa (σ , u), σ∈)

(7)

2 0⎤ ⎥ 4 0⎥ plane stress 0 1⎦ −1 0 ⎤ ⎥ 1 0⎥ plane strain 0 1⎦

(8)

The incompressibility condition, i.e. ΛT ϵ = 0, where Λ = [1 1 0]T , must be enforced for plane strain problem, ensuring that the plastic dissipation D (ϵ (u)) is finite. This kinematic constraint results in a reduction of the number of degrees of freedom, and hence the volumetric locking problem occurs in the kinematic formulation using low-order approximation of displacement fields.

3. Indirect/integrated radial basis function (iRBF) method In this approach, usually the highest-order (second-order in the present case) derivative of the original approximate function is computed first. The first-order partial derivative and the original function are then constructed by successive integration as follows [29]: 5

u,hij (x) =

∑ gI (x) aI = H2 (x) b

(9)

I=1

5

u,hi (x) =

5 + p1

∫ ∑ gI (x) aI dxj = ∑ I=1

5

uh (x) =

H1I (x) aI = H1(x) b

I=1

I=1

(10)

5 + p2

∫ ∫ ∑ gI (x) aI dxj dxi = ∑

H0I (x) aI = H 0 (x) b

I=1

(11)

where 5 denotes the number of nodes in the local support domain, b is a vector consisting of the vector of coefficients and integration constants, and

p2

H1(x) = [H11 (x), H12 (x), …, H15 + p1 (x),  0, …, 0 ] p1

H 0 (x) = [H01 (x), H02 (x), …, H05 + p2 (x)]

where the so-called yield criterion ψ (σ ) is convex. The kinematic principle of yield design, the dual form of (4), is u ∈ℓ

⎧ ⎡4 ⎪ 1⎢ ⎪ 3⎢2 ⎪ ⎣0 Θ=⎨ ⎪⎡ 1 ⎪ ⎢ −1 ⎪⎢ ⎩⎣ 0

(3)

⎪

λ = minD (u)

ϵT Θϵ

H2 (x) = [g1 (x), g2 (x), …, g5 (x),  0, …, 0 ]

where ϵ (u) = [ϵxx ϵ yy γxy ]T are strain rates. The static principle of yield design can now be expressed as

λ = max λ− ⎧ λ−F (u) = a (σ , u), ∀ u ∈ @ s. t ⎨ ⎪ ⎩ ∃ σ ∈ ), ) = { σ ∈ Σ |ψ (σ (x)) ≤ 0 ∀ x ∈ Ω}

∫Ω σp

where sp is the yield stress and

(1)

where the internal and external work rates can be respectively expressed as

a (σ , u) =

D (ϵ (u)) =

93

(6)

For the von Mises criterion, the power of dissipation can be formulated in terms of strain rates as

(12)

(13)

in which p1 and p2 are the number of integration constants, and p2 = 2p1; detailed calculation of H1I (x) and H0I (x) can be found in [28,29]. In this study, multi-quadrics (MQ) basis function, which is ranked the best in terms of accuracy among RBFs [13,28], will be employed. MQs are defined by

gI (x) =

rI2 + cI2

(14)

where rI = ∥ x − xI ∥, the shape parameter cI = βdI (β is a positive scalar), and dI is the minimum of distances from node I to its neighbors in the support domain. The evaluation of Eq. (11) at a set of collocation points yields the following equation in matrix form

94

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

u = HQ b

(15)

5

D=

where

⋯ ⋯ ⋯ ⎡ ⋯ ⎤ HQ = ⎢ H01 (x k ) H02 (x k ) ⋯ H05 + p2 (x k )⎥ ⎢⎣ ⋯ ⎥⎦ ⋯ ⋯ ⋯

(16)

By inversion, the vector b is expressed in terms of nodal function values u as

b=

H−Q1u

(17)

Substituting b into Eqs. (9)–(11), we obtain

∫Ω σ0

(18)

u,hi (x) = H1(x) H−Q1u = Φ, i u

(19)

u,hij (x) = H2 (x) H−Q1u = Φ, ij u

(20)

in which the shape functions and its derivatives are defined by

∑ σ0 AI

(27)

⎧⎡ ρ ⎤ ⎡ 2 0 0⎤ ⎪ ⎢ 1⎥ ⎥ 1 ⎢ ⎪ ⎢ ρ2 ⎥ = ⎢ 1 3 0⎥ BI d plane stress 3⎢ ⎪ ⎢⎣ ρ3 ⎥⎦ ⎣ 0 0 1⎥⎦ ρI = ⎨ ⎪⎡ ⎤ ⎡ ⎤ ⎪ ⎢ ρ1⎥ = ⎢ B xxI d − B yyI d⎥ plane strain ⎪ ⎢⎣ ρ2 ⎥⎦ ⎢⎣ ⎥⎦ 2B xyI d ⎩

(28)

Hence, the optimization problem (5) associated with the iRBF method can now be rewritten as

⎪

(21)

I=1

where l χIk is the (I,k) element of the matrix H−Q1. It should be noted that the iRBF shape functions satisfy the Kronecker-delta property, and hence boundary conditions can be enforced in a way similar to one in the finite element method.

4. Dual yield design formulation using the iRBF method 4.1. iRBF discretization of kinematic formulation In a numerical upper-bound yield design problem, the displacement fields of the problem domain can be approximated by discretization methods. The approximate displacement fields can be expressed in terms of nodal displacement values as

⎡u ⎤

∑ ΦI (x) ⎢ vII ⎥ I=1

∥ ρI ∥

⎧ d = 0 on Γ u s. t ⎨ ⎪ ⎩ F (d) = 1

∑ gI (x)lχIk

5

∑ σ0 AI I=1

I=1 5

⎢⎣

⎥⎦

(22)

where ΦI (x) are the iRBF shape functions described above. The strain rates are then calculated by

⎡ϵ ⎤ ⎢ xx ⎥ ϵ (x) = ⎢ ϵ yy ⎥ = ⎢ γxy ⎥ ⎣ ⎦

∥ ρI ∥

I=1

5

∑ H1I (x)lχIk

⎡ ⎤ u uh (x) = ⎢ ⎥ = ⎢⎣ v ⎥⎦

(26)

5

D=

λ + = min

I=1 5

Φk, ij (x) =

(BI d)T ΘBI d

I=1

where AI is the area of the Ith nodal representative domain, i.e. Voronoi cell. In fact, a form containing a sum of norms can be used to compute the power of internal dissipation as follows:

5

∑ H0I (x)lχIk

Φk, i (x) =

∑ σ0 AI

where ρI are additional variables defined by

uh (x) = H 0 (x) H−Q1u = Φu

Φk (x) =

ϵT Θϵ dΩ =

⎡ ΦI, x (x) 0 ⎤⎡ ⎤ ⎥⎢ u ⎥ ⎢ ∑ ⎢ 0 ΦI, y (x)⎥ ⎢ vII ⎥ = B (x) d I = 1 ⎢ Φ (x) Φ (x) ⎥ ⎢ ⎥ ⎦⎣ ⎦ ⎣ I, y I, x

(29)

Introducing auxiliary variables t1, t2, … , t5 , the problem (29) can be formulated in the form of a standard conic programming as 5

λ + = min

∑ σ0 AI tI I=1

⎧ d = 0 on Γu ⎪ s. t ⎨ F (d) = 1 ⎪ ⎩ ρI ≤ tI I = 1, 2, …, 5

(30)

Note that for plane strain problems, incompressibility conditions, ΛT ϵ = 0, must be introduced. If low-order displacement approximations are used, volumetric locking phenomena in the kinematic formulations associated with the von Mises may occur due to these incompressibility conditions. However, here the iRBF method results in high-order displacement fields, and hence volumetric locking problem can be prevented. Moreover, it is evident that the size of optimization problem (30) depends on the number of integration points to be used. In this paper, the nodal integration technique is used, and hence the size of the resulting optimization problem is kept to be minimum.

4.2. iRBF discretization of static formulation

5

(23)

where d is the nodal displacement vector and B is the so-called strain-displacement matrix, and they are given by

dT = [u1, u2 , …, u5 , v1, v2, …, v5 ] ⎡ B xx ⎤ ⎡ Φ1, x Φ2, x ⋯ Φ5, x 0 0 ⋯ 0 ⎤ ⎥ ⎥ ⎢ ⎢ B = ⎢ B yy ⎥ = ⎢ 0 0 ⋯ 0 Φ1, y Φ2, y ⋯ Φ5, y ⎥ ⎥ ⎢B ⎥ ⎢Φ ⎣ xy ⎦ ⎣ 1, y Φ2, y ⋯ Φ5, y Φ1, x Φ2, x ⋯ Φ5, x ⎦ The internal dissipation power is computed as

(24)

(25)

While in the upper bound formulation the displacement fields are approximated, here the stress fields need to be approximated. With the use of the iRBF method, approximations of these stress fields can be presented as

⎡ σh ⎤ ⎢ xx ⎥ h ⎥ = σ h (x) = ⎢ σyy ⎢ ⎥ h ⎢⎣ σxy ⎥⎦

⎤ ⎡ ⎢ σxxI ⎥ ∑ ΦI (x) ⎢⎢ σyyI ⎥⎥ = Cs σ I=1 ⎢⎣ xyI ⎥⎦ 5

(31)

where

sT = [σxx1, σxx2, …, σxx5 , σyy1, σyy2, …, σyy5 , σxy1, σxy2, …, σxy5 ]

(32)

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

⎡ C xx ⎤ ⎡ Φ Φ ⋯ Φ 0 0 ⋯ 0 0 0 ⋯ 0⎤ 5 ⎥ ⎥ ⎢ 1 2 ⎢ C C = ⎢ yy ⎥ = ⎢ 0 0 ⋯ 0 Φ1 Φ2 ... Φ5 0 0 ⋯ 0 ⎥ ⎢C ⎥ ⎢ 0 0 ⋯ 0 0 0 ⋯ 0 Φ Φ ⋯ Φ ⎥ ⎣ xy ⎦ ⎣ 1 2 5⎦

(33)

These approximated stress fields must be ensured to be statically admissible, meaning that equilibrium and continuity conditions within elements and on their boundary must be satisfied. While the strong form of equilibrium equations can be treated using collocation method, its equivalent weak form (involving integrals) is often handled using the weighted residual method. The strong-form method is simple and fast, and hence the collocation method using the iRBF will be considered in this study. The equilibrium equations can be imposed at 5 nodes, and are expressed as

A1σ 1 + A2σ 3 = 0

(34)

A1σ 3 + A2σ 2 = 0

(35)

where

⋯ ⋯ ⋯ ⎤ ⎡ ⋯ A1 = ⎢ Φ1, x (x k ) Φ2, x (x k ) ⋯ Φ5, x (x k )⎥ ⎢⎣ ⋯ ⋯ ⋯ ⋯ ⎥⎦5 × 5

(36)

⋯ ⋯ ⋯ ⎤ ⎡ ⋯ A2 = ⎢ Φ1, y (x k ) Φ2, y (x k ) ⋯ Φ5, y (x k )⎥ ⎢⎣ ⋯ ⋯ ⋯ ⋯ ⎥⎦5 × 5

(37)

T σ 1 = ⎡⎣ σxx1 ⋯ σxx5 ⎤⎦

(38)

T σ 2 = ⎡⎣ σyy1 ⋯ σyy5 ⎤⎦

(39)

T σ 3 = ⎡⎣ σxy1 ⋯ σxy5 ⎤⎦

(40)

{ ρ ∈ R |ρ ≥ ∥ ρ { ρ ∈ R |ρ ≥ ∥ ρ 3

3

1

1

2 → 4 ∥L2 =

ρ 22 + ρ32 + ρ42

2 → 3 ∥ L2 =

ρ 22 + ρ32

}

}

plane stress

plane strain

(41)

where

ρ1 = σ 0

ρ2 → 4

⎡ρ ⎤ ⎡2 − 1 0 ⎤ ⎢ 2⎥ 1 ⎢ ⎥ ρ = ⎢ 3⎥ = ⎢ 0 3 0 ⎥ Cs plane stress 2 ⎢⎣ ρ4 ⎥⎦ ⎢⎣ 0 0 2 3 ⎥⎦

⎤ ⎡ ⎤ ⎡1 ⎢ ρ2 ⎥ ⎢ (C xx s − C yy s)⎥ ρ2 → 3 = ⎢ ρ ⎥ = ⎢ 2 ⎥ plane strain 3 ⎥⎦ ⎢⎣ ⎥⎦ ⎢⎣ C xy s

equilibrium equations and yield criterion are enforced at nodes only, and therefore the strict property of the lower bound λ
is not guaranteed. However, using a fine nodal distribution one can hope to achieve a reliable approximated lower bound on the actual limit load multiplier. Moreover, by enforcing the equilibrium equations and yield criterion at nodes only the number of constraints in optimization problem (45) is kept to be minimum, and hence the presented static method is computationally inexpensive.

5. Numerical examples The described procedures are tested by their application to solve various problems for which, in most cases, exact and numerical solutions are available. Upper and lower bound solutions based on direct radial basis function (dRBF) are also carried out for comparison purpose. Optimization problems (30) and (45) are implemented in the Matlab environment. Mosek optimization solver version 6.0 is used to solve the conic optimization problem obtained (using a 2.8 GHz Intel Core i5 PC running Window 7). 5.1. Prandtl problem

Additionally, the approximated stress fields must belong to a convex domain, ) . In other words, these stress fields must satisfy the following second-order cone constraints obtaining from the von Mises criterion ⎧ ⎪ 3 PS = σ h ( x) ∈ ) ≡ ⎨ ⎪3 = ⎩ PD

(42)

(43)

(44)

The first example is the classical punch problem presented in [30], as shown in Fig. 1. Due to symmetry, a rectangular region of dimensions B ¼5 and H¼ 2 is considered. Appropriate displacement and stress boundary conditions are enforced as shown in Fig. 2. For a load of 2τ0, the analytical limit multiplier is λ = 2 + π = 5.142. Approximations of upper and lower bounds on the actual limit load for both dRBF and iRBF methods with various nodal discretizations are reported in Table 1. From these results, it can be seen that for both kinematic and static formulations the iRBFbased method can provide more accurate solutions than the dRBFbased method. Convergence analysis and relative errors in collapse multipliers versus number of variables are also shown in Figs. 3 and 4, respectively. It should be stressed that the mean values of upper and lower approximations obtained using the iRBF-based numerical procedures are in excellent agreement with the analytical solutions for all nodal discretizations, as shown in Fig. 3, with less than 0.4% even for coarse nodal distribution. Furthermore, as mentioned, the present procedure cannot theoretically provide strict lower bound solutions, it is evident that all approximated lower bound results are below the exact value. Note that in the kinematic formulation, volumetric (or isochoric) locking often occurs when adding the incompressibility condition to the low-order displacement based yield design problem. The volumetric locking behavior of the Prandtl yield design problem has been studied in [35,12]. In these papers, it has been demonstrated that when smoothed strains were used, the volumetric locking problem can be eliminated. Here, we have shown that the iRBF method used in combination with direct nodal integration can remove such the volumetric locking behavior and also result in stable and accurate solutions. In Table 2 the solutions obtained using the present methods

Hence the static yield design formulation (4) can now be expressed as

λ− = max λ ⎧ A1σ 1 + A2σ 3 = 0 ⎪ s. t ⎨ A1σ 3 + A2σ 2 = 0 ⎪ k ⎩ ρ ∈ 3k, k = 1, 2, …, 5

95

(45)

and accompanied by appropriate boundary conditions. It should be emphasized that in the present static formulation

Fig. 1. Prandtl problem: geometry and loading.

96

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

7 6.5

Collapse multiplier

6 5.5 5 4.5

Analytical solution dRBF−Upper bound iRBF−Upper bound dRBF−Lower bound iRBF−Lower bound dRBF−Mean value iRBF−Mean value

4 3.5 3 2.5

0

200

400 600 Number of nodes

800

1000

Fig. 3. Bounds on the collapse multiplier versus the number of nodes.

50

Fig. 2. Prandtl problem: appropriate displacement and stress boundary conditions.

Table 1 Prandtl problem: upper and lower bound approximation of collapse multiplier.

Relative error in collapse multiplier

45 dRBF−Upper bound dRBF−Lower bound iRBF−Upper bound iRBF−Lower bound

40 35 30 25 20 15 10 5

Nodes

dRBF

iRBF

0

40 160 360 640 1000

Upper-bound

Lower-bound

Upper-bound

Lower-bound

λþ

e (%)

λ


e (%)

λþ

e (%)

λ


e (%)

6.857 5.548 5.289 5.211 5.189

33.35 7.90 2.86 1.34 0.91

2.761 4.244 4.800 5.042 5.125

46.31 17.46 6.65 1.95 0.33

6.036 5.279 5.209 5.191 5.180

17.39 2.66 1.30 0.95 0.74

4.218 4.960 5.061 5.108 5.134

17.97 3.54 1.58 0.66 0.16

with 2560 nodes are compared with those obtained previously by different yield design approaches using FEM, smoothed finite element (SFEM) and EFG simulations. In general, the present solutions are close to results in the literature. Considering upper solutions, the result obtained using the iRBF method is slightly lower than the one obtained using the EFG mesh-free method with the same nodal discretization [12]. 5.2. Thin square plates with cutouts subjected to tension load Next, two thin square plates with a central square cutout and a thin crack subjected to a uniform tension load, as shown in Fig. 5, are considered. These problems have been investigated numerically by finite elements [37,39], symmetric Galerkin boundary elements [38], and mesh-free methods [40,41]. Owing to the symmetry, only the topright quarter of plates is modeled, as shown in Fig. 6. Uniform nodal distribution is used to discretize the computational domain, see Fig. 7.

0

1000

2000 3000 4000 number of variables

5000

6000

Fig. 4. Relative error in collapse multipliers versus the number of variables.

Table 2 The punch of problem: comparison with previous solutions. Authors

Present method (dRBF) Present method (iRBF) Vicente da Silva and Antao [31] Sloan and Kleeman [32] Makrodimopoulos et al. [33,34] Le et al. [12], EFG Le et al. [35], SFEM Capsoni and Corradi [36]

Approach

Kinematic and static Kinematic and static Kinematic Kinematic Kinematic and static Kinematic Kinematic Mixed formulation

Collapse multiplier Upperbound

Lowerbound

5.159

5.133

5.146

5.140

5.264

–

5.210 5.148

– 5.141

5.147 5.143 5.240

– –

Limit load multipliers obtained using uniform nodal distributions are reported in Tables 3 and 4. Collapse load multiplier versus the number of nodes is also shown in Fig. 8. Again, it can be observed that the iRBF-based method can provide more accurate solutions than the dRBF-based method, particularly for the static

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

97

Fig. 5. Thin square plates.

Fig. 6. The upper-right quarter of plates.

Fig. 7. Uniform nodal discretization.

approach. Table 5 shows that the results obtained by using RBF methods are in good agreement with previously reported numerical solutions. Considering upper bound limit factor, the present results are close to Zhou and Liu's solutions, with the maximum error of only 2.79%. It is important to note that the estimated lower bounds reported in [38,40] are higher than the present lower bound solutions, and surpasses the upper bound of the present iRBF method for the plate with square cutout. This can be explained by the fact that in [38,40] the strong form of the equilibrium equations was transformed into the so-called weak form, and to be satisfied locally in an average sense using approximated virtual

Table 3 Collapse multipliers for the square plate with a central square cutout. Nodes

273 589 851

dRBF

iRBF

Upper-bound

Lower-bound

Upper-bound

Lower-bound

0.750 0.750 0.749

0.734 0.732 0.732

0.702 0.715 0.715

0.712 0.726 0.729

98

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

Table 4 Collapse multipliers for the square plate with a central thin crack. Nodes

441 625 841

dRBF

Table 5 Plates with cutouts problem: comparison with previous solutions.

iRBF

Authors

Upper-bound

Lower-bound

Upper-bound

Lower-bound

0.520 0.516 0.514

0.482 0.484 0.487

0.532 0.523 0.516

0.499 0.502 0.504

Approach

Present method (dRBF) Present method (iRBF) Belytschko and Hodge [37] Zhang et al. [38] Zhang et al. [39] Chen et al. [40] Zhou and Liu (Sibson) [41] Zhou and Liu (Laplace) [41]

Square cutout

Thin crack

UB

LB

UB

LB

Kinematic and static Kinematic and static Static

0.749

0.715

0.514

0.487

0.732

0.729

0.516

0.504

–

0.693

–

0.498

Static Kinematic Static Kinematic

– 0.764 – 0.753

0.747 – 0.736 –

– 0.534 – 0.515

0.514 – 0.513 –

Kinematic

0.752

–

0.513

–

including higher-order displacement-based finite element method [43], mixed finite elements [36,44,45] and discontinuous elements [46,32,47], mesh-free methods [12], smoothed finite elements [35,48]. Owing to the symmetry, only the upper-right quarter of the double notched problem is discretized. Several uniform nodal distributions are employed. Computed solutions and convergence analysis are presented in Table 6 and Fig. 10. Table 7 compares the present solutions with those obtained previously. The mean values of the dRBF and iRBF results are 1.1343 and 1.1342, respectively. It can be observed that these mean values are very close to the benchmark solution obtained using mixed formulation by Christiansen and Andersen [45].

Fig. 8. Convergence of limit load factor for the plates.

displacement fields. Therefore, the static method in [38,40] may result in a higher value than the actual limit multiplier. In contrast, it is clear that all the present lower bound solutions obtained are below the upper bounds reported in Table 5. 5.3. Notched tensile specimen Finally, a double notched specimen consists of a rectangular specimen with two thin cracks under in-plane tensile stresses τ0, as shown in Fig. 9, is also considered. This problem exhibits volumetric locking phenomena [42] and became a popular benchmark test for plastic yield design procedures. The locking problem was handled using various techniques proposed in the literature,

Fig. 9. Double notched specimen.

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

and lower bounds are accurate, and can be considered as the actual collapse load multiplier for most practical engineering problems, for which exact solution is unknown.

Table 6 Double notched specimen: limit load multiplier. Nodes

289 441 625 841

dRBF

99

iRBF

Upper-bound

Lower-bound

Upper-bound

Lower-bound

λþ

e (%)

λ


e (%)

λþ

e (%)

λ


e (%)

1.156 1.152 1.149 1.146

0.21 0.33 0.29 0.20

1.092 1.116 1.120 1.122

5.56 2.24 0.37 0.20

1.150 1.146 1.143 1.141

0.56 0.34 0.23 0.15

1.100 1.114 1.122 1.127

2.79 1.31 0.66 0.48

Acknowledgments This research has been supported by the Institute for Computational Science and Technology (ICST), Ho Chi Minh City.

References

Fig. 10. Convergence study for the double notched specimen problem.

Table 7 The double notched specimen: comparison with previous solutions. Authors

Present method (dRBF) Present method (iRBF) Ciria et al. (Uniform) [49] Ciria et al. (Adaptive) [49] Le et al. [35] Le et al. [12] Krabbenhoft et al. [50] Tin-Loi and Ngo [43] Christiansen et al. [45]

Approach

Kinematic and static Kinematic and static Kinematic and static Kinematic and static Kinematic Kinematic Static Static Mixed formulation

Collapse multiplier Upperbound

Lowerbound

1.146 1.141 1.149 1.139 1.137 1.154 – – 1.136

1.122 1.127 1.131 1.132 – – 1.132 1.166

6. Conclusions The present contribution has presented displacement and equilibrium mesh-free formulation based on integrated radial basis functions (iRBF) for dual yield design problems. In the kinematic formulation, the high-order approximation of the displacement fields using the integrated radial basis functions can prevent volumetric locking. Moreover, direct nodal integration of the iRBF approximation not only results in inexpensive computational cost, but also overcomes the instability problems. In the static formulation, with the use of iRBF approximation of the stress fields in combination with the collocation method, equilibrium equations and yield conditions only need to be enforced at the nodes, leading to the reduction in computational effort. It has been shown in several examples that the mean values of the iRBF upper

[1] Save MA, Massonnet CE, de Saxce G. Plastic analysis and design of plates, shells and disks.North-Holland series in applied mathematics and mechanics, 43. Amsterdam: Elsevier; 1997. [2] Salencon J. Yield design. Wiley.com; 2013. [3] Sloan SW. Lower bound limit analysis using finite elements and linear programming. Int J Numer Anal Methods Geomech 1988;12:61–77. [4] Krenk S, Damkilde L, Hoyer O. Limit analysis and optimal design of plates with equilibrium elements. J Eng Mech 1994;120:1237–54. [5] Poulsen PN, Damkilde L. Limit state analysis of reinforced concrete plates subjected to in-plane forces. Int J Solids Struct 2000;37:6011–29. [6] Krabbenhoft K, Damkilde L. Lower bound limit analysis of slabs with nonlinear yield criteria. Comput Struct 2002;80:2043–57. [7] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [8] Chen SS, Liu YH, Cen ZZ. Lower-bound limit analysis by using the EFG method and non-linear programming. Int J Numer Methods Eng 2008;74:391–415. [9] Le CV, Askes H, Gilbert M. Adaptive element-free Galerkin method applied to the limit analysis of plates. Comput Methods Appl Mech Eng 2010;199:2487– 96. [10] Le CV, Gilbert M, Askes H. Limit analysis of plates using EFG method and second order cone programming. Int J Numer Methods Eng 2009;78:1532–52. [11] Le CV, Gilbert M, Askes H. Limit analysis of plates and slabs using a meshless equilibrium formulation. Int J Numer Methods Eng 2010;83:1739–58. [12] Le CV, Askes H, Gilbert M. A locking-free stabilized kinematic EFG model for plane strain limit analysis. Comput Struct 2012;106–107:1–8. [13] Franke R. Scattered data interpolation: tests of some methods. Math Comput 1982;38:181. [14] Franke R, Schaback R. Solving partial differential equations by collocation using radial functions. Appl Math Comput 1998;93:73–82. [15] Fasshauer GE. Solving differential equations with radial basis functions: multilevel methods and smoothing. Adv Comput Math 1999;11:139–59. [16] Mai-Duy N, Tran-Cong T. Numerical solution of differential equations using multiquadric radial basis function networks. Neural Netw 2001;14(2):185–99. [17] Cheng AHD, Golberg MA, Kansa EJ, Zammito G. Exponential convergence and H-c multiquadrics collocation method for partial differential equations. Numer Methods Partial Diff Equ 2003;19:571–94. [18] Hu HY, Li ZC, Cheng AHD. Radial basis collocation method for elliptic equations. Comput Math Appl 2005;50:289–320. [19] Hu HY, Chen JS, Hu W. Weighted radial basis collocation method for boundary value problems. Int J Numer Methods Eng 2007;69:2736–57. [20] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246(1):29–46. [21] Xiao JR, McCarthy MA. A local heaviside weighted meshless method for twodimensional solids using radial basis functions. Comput Mech 2003;31:301– 15. [22] Liu GR, Zhang GY, Gu YT, Wang YY. A meshfree radial point interpolation method (RPIM) for three-dimensional solids. Comput Mech 2005;36:421–30. [23] Chen JS, Wang L, Hu HY, Chi SW. Subdomain radial basis collocation method for heterogeneous media. Int J Numer Methods Eng 2009;80:163–90. [24] Xia P, Long SY, Cui HX, Li GY. The static and free vibration analysis of a nonhomogeneous moderately thick plate using the meshless local radial point interpolation method. Eng Anal Bound Elem 2009;33:770–7. [25] Wang L, Chen JS, Hu HY. Subdomain radial basis collocation method for fracture mechanics. Int J Numer Methods Eng 2010;83:851–76. [26] Khosravifard A, Hematiyan MR, Marin L. Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method. Appl Math Model 2011;35:4157–74. [27] Christiansen E. Limit analysis of collapse states. In: Handbook of numerical analysis, vol. IV. Amsterdam: North-Holland; 1996. p. 193312 [chapter II]. [28] Mai-Duy N, Tran-Cong T. Approximation of function and its derivatives using radial basis function networks. Appl Math Model 2003;27:197–220. [29] Mai-Duy N, Tran-Cong T. An efficient indirect RBFN-based method for numerical solution of PDEs. Numer Methods Partial Diff Equ 2005;21:770–90. [30] Prandtl L. Ueber die Haerte plastischer Koerper. Nachr Akad Wiss Gott II: Math-Phys Kl II 1920;12:74–85. [31] Vicente da Silva M, Antao AN. A non-linear programming method approach for upper bound limit analysis. Int J Numer Methods Eng 2007;72:1192–218. [32] Sloan SW, Kleeman PW. Upper bound limit analysis using discontinuous

100

P.L.H. Ho et al. / Engineering Analysis with Boundary Elements 71 (2016) 92–100

velocity fields. Comput Methods Appl Mech Eng 1995;127:293–314. [33] Makrodimopoulos A, Martin CM. Upper bound limit analysis using simplex strain elements and second-order cone programming. Int J Numer Anal Methods Geomech 2006;31:835–65. [34] Makrodimopoulos A, Martin CM. Lower bound limit analysis of cohesivefrictional materials using second-order cone programming. Int J Numer Methods Eng 2006;66:604–34. [35] Le CV, Nguyen-Xuan H, Askes H, Bordas S, Rabczuk T, Nguyen-Vinh H. A cellbased smoothed finite element method for kinematic limit analysis. Int J Numer Methods Eng 2010;83:1651–74. [36] Capsoni A, Corradi L. A finite element formulation of the rigid-plastic limit analysis problem. Int J Numer Methods Eng 1997;40:2063–86. [37] Belytschko T, Hodge PG. Plane stress limit analysis by finite element. J Eng Mech Div 1970;96:931–44. [38] Zhang XF, Liu YH, Zhao YN, Cen ZZ. Lower bound limit analysis by the symmetric Galerkin boundary element method and the complex method. Comput Methods Appl Mech Eng 2002;191:1967–82. [39] Zhang F, Mingwan L, Kehchih H. A mathematical programming algorithm for limit analysis. Acta Mech Sin 1991;7:267–74. [40] Chen S, Liu Y, Cen Z. Lower-bound limit analysis by using the EFG method and non-linear programming. Int J Numer Methods Eng 2008;74:391–415. [41] Zhou ST, Liu YH. Upper-bound limit analysis based on the natural element method. Acta Mech Sin 2012;28:1398–415. [42] Nagtegaal JC, Parks DM, Rice JC. On numerically accurate finite element solutions in the fully plastic range. Comput Methods Appl Mech Eng 1974;4:153–77. [43] Tin-Loi F, Ngo NS. Performance of the p-version finite element method for limit analysis. Int J Mech Sci 2003;45:1149–66.

[44] Andersen KD, Christiansen E, Overton ML. Computing limit loads by minimizing a sum of norms. SIAM J Sci Comput 1998;19:1046–162. [45] Christiansen E, Andersen KD. Computation of collapse states with von Mises type yield condition. Int J Numer Methods Eng 1999;46:1185–202. [46] Bottero A, Negre R, Pastor J, Turgeman S. Finite element method and limit analysis theory for soil mechanics problems. Comput Methods Appl Mech Eng 1980;22:131–49. [47] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and nonlinear programming. Int J Numer Anal Methods Geomech 2002;26:181–216. [48] Le CV, Nguyen-Xuan H, Askes H, Rabczuk T, Nguyen-Thoi T. Computation of limit load using edge-based smoothed finite element method and secondorder cone programming. Int J Comput Methods 2013;10:1340004. [49] Ciria H, Peraire J, Bonet J. Mesh adaptive computation of upper and lower bounds in limit analysis. Int J Numer Methods Eng 2008;75:899–944. [50] Krabbenhoft K, Damkilde L. A general nonlinear optimization algorithm for lower bound limit analysis. Int J Numer Methods Eng 2003;56:165–84. [51] Le CV, Nguyen-Xuan H, Nguyen-Dang H. Upper and lower bound limit analysis of plates using FEM and second-order cone programming. Comput Struct 2010;88:65–73. [52] Tran DT, Le CV. Extended finite element method for plastic limit load computation of cracked structures. Int J Numer Methods Eng 2015;104:2–17. [53] Le CV, Nguyen PH, Chu TQ. A curvature smoothing Hsieh–Clough–Tocher element for yield design of reinforced concrete slabs. Comput Struct 2015;152:59–65. [54] Le CV, Ho PLH, Nguyen PH, Chu TQ. Yield design of reinforced concrete slabs using a rotation-free meshfree method. Eng Anal Bound Elem 2015;50:231–8.