An adaptive strategy based on conforming quadtree meshes for kinematic limit analysis

TL;DR

This paper introduces an adaptive finite element method using conforming quadtree meshes for efficient kinematic limit analysis of structures with cohesive and frictional materials, addressing computational challenges like hanging nodes.

Contribution

It presents a novel adaptive mesh strategy with conforming quadtree meshes and enhanced shape functions to improve efficiency and reliability in collapse plastic analysis.

Findings

Adaptive mesh strategy improves computational efficiency.

Conforming quadtree meshes resolve hanging node issues.

Numerical results validate the approach's reliability.

Abstract

We propose a simple and efficient scheme based on adaptive finite elements over conforming quadtree meshes for collapse plastic analysis of structures. Our main interest in kinematic limit analysis is concerned with both purely cohesive-frictional and cohesive materials. It is shown that the most computational efficiency for collapse plastic problems is to employ an adaptive mesh strategy on quadtree meshes. However, a major difficulty in finite element formulations is the appearance of hanging nodes during adaptive process. This can be resolved by a definition of conforming quadtree meshes in the context of polygonal elements. Piecewise-linear shape functions enhanced with generalized bubble functions in barycentric coordinates are used to approximate the velocity field. Numerical results prove the reliability and benefit of the present approach.

Tables

Table 1. Table 1 : A smooth strip footing: Comparison with the exact solution using adaptive mesh.

$\phi =0^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	137	807	2107	5737	16542
		Mosek time (s)	0.05	0.06	0.16	0.45	2.70
		${\alpha }^{+}$	5.437	5.203	5.169	5.154	5.149
		Relative error (%)	5.75	1.19	0.53	0.24	0.14
	Uniform mesh (bL-Q4)	$N_{var}$	137	1742	6682	26162	103522
		Mosek time (s)	0.05	0.11	0.58	2.38	18.67
		${\alpha }^{+}$	5.437	5.202	5.171	5.156	5.149
		Relative error (%)	5.75	1.17	0.57	0.28	0.14
$\phi ={35}^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	1417	2467	5332	14647	44017
		Mosek time (s)	0.11	0.16	0.39	2	9.5
		${\alpha }^{+}$	51.922	48.044	46.939	46.488	46.292
		Relative error (%)	12.57	4.16	1.76	0.79	0.36
	Uniform mesh (bL-Q4)	$N_{var}$	1417	5432	21262	84122
		Mosek time (s)	0.19	0.39	2.05	11.03
		${\alpha }^{+}$	51.922	48.065	46.971	46.519
		Relative error (%)	12.57	4.21	1.84	0.86

Table 2. Table 2 : A smooth strip footing: Comparison with other methods.

Approach	Authors	$\phi =0^{\circ }$		$\phi ={35}^{\circ }$
		${\alpha }^{+}$	(error (%))	${\alpha }^{+}$	(error %))
Kinematic	Sloan & Kleeman [10]	5.21	(-1.33)	–	–
	Makrodimopoulos & Martin [17]	5.148	(-0.12)	46.37	(-0.52)
	Present (Adaptive bL-Q4)	5.148	(-0.14)	46.29	(-0.36)
Mixed	Capsoni and Corradi [9]	5.24	(-1.91)	–	–
Static	Makrodimopoulos & Martin [16]	5.141	(-0.02)	46.07	(-0.12)
Analytical	Prandtl	2 + $\pi$	0	46.124	0

Table 3. Table 3 : A smooth strip footing: limit load factor for various internal frictional angles.

${\phi }^{\circ }$	${\alpha }^{+}$ Present Method	${\alpha }^{+}$ Prandt [70]	Error ($\%$)
$0$	$5.145$	$5.142$	$0.060$
$5$	$6.493$	$6.489$	$0.066$
$10$	$8.350$	$8.345$	$0.065$
$15$	$10.986$	$10.977$	$0.084$
$20$	$14.845$	$14.835$	$0.072$
$25$	$20.735$	$20.721$	$0.072$
$30$	$30.166$	$30.140$	$0.086$
$35$	$46.160$	$46.124$	$0.079$
$40$	$75.384$	$75.313$	$0.095$
$45$	$134.098$	$133.874$	$0.168$

Table 4. Table 4 : Block with two symmetrical circular holes: The convergence of the present solutions using adaptive meshes.

$\phi =0^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	1821	6028	25418	56818
		${\alpha }^{+}$	1.831	1.828	1.816	1.814
		Relative error (%)	1.05	0.89	0.23	0.12
	Uniform mesh (bL-Q4)	$N_{var}$	1821	6668	25818	101558
		${\alpha }^{+}$	1.831	1.832	1.825	1.818
		Relative error (%)	1.05	1.11	0.72	0.34
$\phi ={30}^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	1821	5713	28478	67458
		${\alpha }^{+}$	1.073	1.063	1.06	1.059
		Relative error (%)	1.41	0.46	0.18	0.09
	Uniform mesh (bL-Q4)	$N_{var}$	1821	6668	25818	101558
		${\alpha }^{+}$	1.073	1.068	1.062	1.060
		Relative error (%)	1.41	0.94	0.37	0.18

Table 5. Table 5 : Block with two symmetric holes: Comparison with literature solutions.

Approach	Authors	$\phi =0^{\circ }$	$\phi ={30}^{\circ }$
Kinematic	Makrodimopoulos &Martin [17]	1.825	1.063
	Munoz et al. [28]	1.8351	1.0652
	Present (Adaptive bL-Q4)	1.814	1.059
Static	Makrodimopoulos &Martin [16]	1.8089	1.0562
	Munoz et al. [28]	1.8119	1.0581
Reference	Zouain et al [18]	1.8131	–

Table 6. Table 6 : Slope stability: Comparison with the exact solution using adaptive mesh.

$\phi ={20}^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	1492	1837	2707	5332	13212
		${\alpha }^{+}$	9.496	8.791	8.449	8.313	8.266
		Relative error (%)	15.66	7.08	2.91	1.25	0.68
	Uniform mesh (bL-Q4)	$N_{var}$	1492	5702	22282	88082
		${\alpha }^{+}$	9.496	8.83	8.543	8.393
		Relative error (%)	15.66	7.55	4.06	2.23
$\phi ={35}^{\circ }$	Present approach (Adaptive bL-Q4)	$N_{var}$	1492	1772	2507	4287	9027
		${\alpha }^{+}$	20.075	16.38	14.842	14.195	13.984
		Relative error (%)	46.00	19.13	7.94	3.24	1.70
	Uniform mesh (bL-Q4)	$N_{var}$	1492	5702	22282	88082
		${\alpha }^{+}$	20.075	16.457	15.028	14.407
		Relative error (%)	46.00	19.69	9.29	4.78

Table 7. Table 7 : Slope stability: Comparison with other methods in the literature.

Approach	Authors	$\phi ={20}^{\circ }$	$\phi ={35}^{\circ }$
Kinematic	Krabbenhoft et al.[13]	8.440	–
	Lyamin & Sloan [15]	8.440	–
	Makrodimopoulos & Martin [17]	8.366	14.19
	Chen [72]	8.300	13.86
	Present (Adaptive bL-Q4)	8.266	13.984
Static	Makrodimopoulos & Martin [16]	8.210	13.75

Equations

$$W_{int}(\mathbf{\sigma },\dot{\mathbf{u}})=W_{ext}(\dot{\mathbf{u}}),\forall \dot{\mathbf{u}}\in \mathcal{V}$$

$$W_{int}(\mathbf{\sigma },\dot{\mathbf{u}})={\int }_{\Omega }\mathbf{\sigma }:\mathbf{\epsilon }(\dot{\mathbf{u}})\mathrm{d}\Omega$$

$$W_{ext}(\dot{\mathbf{u}})={\int }_{\Omega }\mathbf{f}\cdot \dot{\mathbf{u}}\mathrm{d}\Omega +{\int }_{{\Gamma }_t}\mathbf{g}\cdot \dot{\mathbf{u}}\mathrm{d}\Gamma$$

$$\mathcal{V}=\{\dot{\mathbf{u}}\in {\mathcal{H}}^1(\Omega {)}^2,\dot{\mathbf{u}}=\dot{\bar{\mathbf{u}}}\text{on}{\Gamma }_{\dot{u}}\}$$

$$\eta=\|\dot{\bm{\varepsilon}}\|_{\mathcal{L}^{2}(\Omega)}=\Big{(}\int_{\Omega}\mathit{\dot{\varepsilon}}_{ij}\mathit{\dot{\varepsilon}}_{ij}\text{d}\Omega\Big{)}^{1/2}$$

$$\dot{\mathbf{\epsilon }}=\dot{\mu }\frac{\partial \psi (\mathbf{\sigma })}{\partial \mathbf{\sigma }}$$

$$\displaystyle\mathcal{B}=\{\mathbf{\bm{\sigma}\in\sum\big{|}\psi(\bm{\sigma})\leq 0}\}$$

$$\displaystyle{\color[rgb]{0,0,1}\lambda^{*}}$$

$$=\underset{\mathbf{\sigma }\in \mathcal{B}}{\max }\underset{\dot{\mathbf{u}}\in \mathcal{C}}{\min }{W}_{int}(\mathbf{\sigma },\dot{\mathbf{u}})$$

$$=\underset{\dot{\mathbf{u}}\in \mathcal{C}}{\min }\underset{\mathbf{\sigma }\in \mathcal{B}}{\max }{W}_{int}(\mathbf{\sigma },\dot{\mathbf{u}})$$

$$=\underset{\dot{\mathbf{u}}\in \mathcal{C}}{\min }D(\dot{\mathbf{u}})$$

$$D(\dot{\mathbf{u}})=\underset{\mathbf{\sigma }\in \mathcal{B}}{\max }{W}_{int}(\mathbf{\sigma },\dot{\mathbf{u}})$$

$$\psi (\mathbf{\sigma })$$

$${\sigma }_m$$

$$D(\dot{\mathbf{u}})$$

$$\lambda \ge 2\sqrt{J_2(\dot{\mathbf{e}})}\text{and}{\dot{e}}_m=\varsigma \lambda ,\forall \varsigma \ge 0$$

$${\alpha }^{+}$$

$$s.t.$$

$$\displaystyle{\color[rgb]{0,0,1}\dot{e}^{h}_{m}=0}\quad\text{or}\quad\bm{\nabla}\cdot\dot{\mathbf{u}}=0$$

$$\displaystyle\mathcal{V}^{h}=\Big{\{}\dot{\mathbf{u}}^{h}\in\mathcal{H}^{1}(\Omega)^{2},\dot{\mathbf{u}}^{h}\Big{|}_{\Omega^{e}}\in[\mathcal{Q}(\Omega^{e})]^{2},\dot{\mathbf{u}}^{h}\Big{|}_{\Omega^{e}}\in\mathcal{C}^{\infty},\dot{\mathbf{u}}^{h}\Big{|}_{\partial\Omega^{e}}\in\mathcal{C}^{0},\Omega^{e}\in\mathcal{T}\Big{\}}$$

$$\displaystyle\mathcal{V}^{h}_{b}=\Big{\{}\dot{\mathbf{u}}^{h}_{b}\in\mathcal{H}^{1}(\Omega)^{2},\dot{\mathbf{u}}^{h}_{b}\Big{|}_{\partial\Omega^{e}}=0\ \text{and}\ \dot{\mathbf{u}}^{h}_{b}|_{\Omega^{e}}\in\mathcal{C}^{0}\Big{\}}$$

$${\mathcal{W}}^h={\mathcal{V}}^h\oplus {\mathcal{V}}_b^h$$

$$\dot{\mathbf{u}}^{h}(\mathbf{x})\Big{|}_{\Omega^{e}}=\sum_{i=1}^{n_{nod}}(N^{e}_{i}(\mathbf{x})\mathbf{I}_{2})\dot{\mathbf{d}}^{e}_{i}+(N^{e}_{b}(\mathbf{x})\mathbf{I}_{2})\dot{\mathbf{d}}^{e}_{b}$$

$${\varphi }_i^e({\mathbf{x}}_j)={\delta }_{ij}=\{\begin{array}{l}1{\mathbf{x}}_i={\mathbf{x}}_j\\ 0{\mathbf{x}}_i\ne {\mathbf{x}}_j\end{array}$$

$${\varphi }_i^e({\mathbf{x}}_c)=1/n_{nod}$$

$$N_i^e(\mathbf{x})$$

$$\nabla N_i^e(\mathbf{x})$$

$${\dot{\mathbf{\epsilon }}}^{(k)}={\mathcal{P}}^{(k)}({\nabla }_s{\dot{\mathbf{u}}}^h)=\frac{1}{A^{(k)}}{\int }_{{\Omega }^{(k)}}{\nabla }_s{\dot{\mathbf{u}}}^h\mathrm{d}\Omega ,\text{with}{\dot{\mathbf{u}}}^h\in {\mathcal{W}}^h$$

$${\nabla }_s={\left[\begin{array}{ccc}{\partial }_{,x}& 0& {\partial }_{,y}\\ 0& {\partial }_{,y}& {\partial }_{,x}\end{array}\right]}^{\mathrm{T}}$$

$${\dot{\bar{\epsilon }}}_m^h={\mathcal{P}}^{(k)}({\dot{\epsilon }}_m^h)={\mathcal{P}}^{(k)}(\nabla \cdot {\dot{\mathbf{u}}}^h)=\frac{1}{A^{(k)}}{\int }_{{\Omega }^{(k)}}\nabla \cdot {\dot{\mathbf{u}}}^h\text{d}\Omega ,\text{with}{\dot{\mathbf{u}}}^h\in {\mathcal{W}}^h$$

$$\forall {\Omega }^{(k)}\in \bar{\mathcal{T}},{\dot{\bar{\epsilon }}}_m^h=0$$

$$\displaystyle{\color[rgb]{0,0,1}\dot{\bar{\bm{\varepsilon}}}^{h}=\sum_{I=1}^{n_{nk}+n_{b}}\bar{\mathbf{B}}^{(k)}_{I}\dot{\bar{\mathbf{d}}}^{(k)}_{I}=\bar{\mathbf{B}}^{(k)}\dot{\bar{\mathbf{d}}}^{(k)}}$$

$${\bar{\mathbf{B}}}_I^{(k)}=\frac{1}{A^{(k)}}\sum _{i=1}^{\text{nc}}{\mathbf{B}}_i^{\ast },\text{nc is the number of sub-triangles in the edge k}$$

Full text

1 Introduction

Limit analysis has been known as a power tool to directly obtain the ultimate load bearing capacity and plastic collapse path of structures without any requirement of iterative or incremental analysis. The iteratively elastic-plastic analysis that might be undertaken with numerical methods is capable of providing a collapse load factor. However, the computational cost and convergence of the nonlinear solution is still questionable for large-scale structures. As an alternative methodology, limit analysis allows to show the most important features of the limit state of structures in the plastic regime. Until now, the slip-line field (SLF) theory is an analytical approach for evaluating the load bearing capacity of structures [1, 2, 08Miller]. Practically, SLF is well suited to the given problem with simple geometry and loading conditions [4]. The extensive development of different numerical methods has become more attractive beyond analytical approaches. Various numerical methods basically rely on the Koiter’s kinematic (upper bound) theorem [5] or the Melan’s static (or lower bound) theorem [6]. The upper bound limit analysis is to find a load factor that results in minimizing the plastic dissipation problem (without the work of any additional loads) for any kinematically admissible displacement field, while the lower bound limit analysis determines statically and plastically admissible stress field that maximizes the plastic load factor. Among them, finite element methods have become popular in limit analysis [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. Moreover, lower-order finite elements based on triangular or quadrilateral types are largely used due to its simplicity and efficiency, and specially, the flexibility in adaptive mesh refinements. Nevertheless, they do not perform well with unstructured meshes [17] and are sensitive to volumetric locking in the incompressibility limit. Several advanced technologies have therefore been devised in the literature [16, 17, 18, 19].

In attempts to enhance in the accuracy of the limit analysis solution, adaptive mesh refinement becomes rather important. We know that localized plastic deformations cause the slow convergence of the numerical approaches [20]. Therewith, the mesh should be automatically refined along plastic zones. Theoretically, the error-based indicator has to be known to conduct adaptive mesh refinement[21, 22, 23, 24]. However, finding a priori error estimate for the plastic problem is not easy. Hence, a posteriori error scheme is the most suitable choice. Initially, a posteriori indicator based on the recovery technique of the Hessian matrix was studied in [20, 25]. Then, several alternative indicators related to the plastic dissipation and both the static and kinematic bound problems were studied in [26, 27, 28, 29]. The plastic dissipation indicator in association with adaptive bubble-enhanced triangular finite element formulation was presented in [19]. On the other hand, we observe that plastic strain rates can be used to measure localized plastic deformation at limit state. In addition, there exist some regions in the structure represented by highly plastic strain rates. It requires that these zones must be refined to obtain the solution that converges quickly to the actual value with the lowest computational cost.

Regarding advances in mesh generation, quadtree meshes [30] have been well known and applied largely to a wide range of various engineering multi-disciplines, e.g, computer graphics and image processing [31], moving fluid interface[32], large-scale earthquake ground motion simulation [33]. They have been successfully incorporated into finite elements in the framework of adaptive mesh strategy to enhance the accuracy of approximate solution, especially for singular problems with high demand of computational resource or solution formulation [34, 35]. Recently, quadtree meshes in combination with the scaled boundary finite element method (SBFEM) [36, 37] showed a robust computational tool for mechanics problems. However, in the context of the standard FEM, after each refinement, hanging nodes are generated unless new elements and their neighbours are at same levels of refinement. For instance, Fig. 1 (a) illustrates three different levels of refinements. We see that hanging nodes are vertices of the child (smaller) element that is located on the side of the father (larger) element, i.e, diamond-shape nodes displayed in Fig. 1 (a). Consequently, hanging nodes leads to the incompatibilities in finite element approximations. To overcome these difficulties, several special treatments have been studied such as Lagrange multipliers or penalty method [38], adding temporary elements [39], constraining hanging nodes to corner nodes [40], hierarchical enrichment [41], B-splines [42] [43], and natural neighbor basis functions [44]. In [44], natural neighbor (Laplace) basis functions were utilized to obtain $\mathcal{C}^{0}$– continuous approximations along sides involving hanging nodes. The shape functions are constructed from the polygonal reference elements through the affine map. By this approach, we do not require the number of hanging nodes lied on each side and the computations are, moreover, performed over the framework of polygonal elements. We recall that polygonal finite elements have been applied to an extensive class of mechanics problems such as nonlinear constitutive modeling of polycrystalline materials [45, 46, 47], solid mechanics problems [48, 49, 50, 51, 52, 53], topology optimization [54, 55, 56], incompressible materials [57], fracture modeling [58, 59, 60] and so on. However, we emphasize that, in the context of polygonal finite elements, the numerical solutions of limit analysis are less accurate for cohesive–-frictional materials, especially for purely cohesive material or incompressible material due to the presence of the flow rule constraints.

Last but not least, another important issue is how to efficiently solve a minimization problem existing in limit analysis models. Traditionally, linear and non-linear programming is totally possible, but linear programming problems lead to a large number of additional variables. It is interesting that most of the plastic yield criteria can be represented as an intersection of cones where the primal–dual interior point method [62, 63] implemented in optimization packages, e.g. MOSEK software, is really suitable to evaluate the limit analysis problem efficiently. For example, this optimization tool used for the limit analysis problems has been implemented in [16, 17].

The main contribution of this study focuses on some following aspects:

A new bubble-enhanced FE formulation with piecewise-linear shape functions is presented to eliminate the volumetric locking in fully plastic regime.

• 2.

Quadrilateral elements with the presence of hanging nodes are resolved as conforming polygonal elements.

• 3.

Second-order cone programming (SOCP) is exploited to solve the large-scale optimization problems.

• 4.

A quadtree-based adaptive mesh generator guided by the $\mathcal{L}^{2}$ - norm of plastic strain rates is presented.

The paper is outlined as follows: A brief review on the upper bound theorem at hand is described in the next section. Section 3 presents a polygonal finite element formulation implemented with bubble functions through barycentric coordinates and a quadrature scheme based on an area-averaged edge-based projection technique. Section 4 states a solution procedure of the discrete problem. An adaptive quadtree mesh procedure is presented in Section 5. Section 6 summarizes a numerical implementation. Numerical validations are given in Section 7. Section 8 closes the paper with the concluding remarks.

2 A short review of the kinematic theorem

Let us consider a problem domain $\Omega$ bounded by a continuous boundary $\Gamma=\Gamma_{\dot{u}}\cup\Gamma_{t},\Gamma_{\dot{u}}\cap\Gamma_{t}=\emptyset$. The rigid-perfectly plastic body is acted by body forces f and external tractions g on $\Gamma_{t}$, and the boundary $\Gamma_{\dot{u}}$ is prescribed by the displacement velocity vector $\mathbf{\dot{u}}$. The equilibrium equation can be delineated as follows

[TABLE]

in which the internal work rate for stress tensor $\bm{\sigma}$ and velocity vector $\dot{\mathbf{u}}$ is described by

[TABLE]

and the external work rate is defined by

[TABLE]

in which $\mathcal{V}$ is a space of kinematically admissible velocity field denoted by

[TABLE]

We also define $\mathcal{L}^{2}(\Omega)$ -norm for strain rates which is regarded as a mesh-control indicator in adaptive mesh refinements as

[TABLE]

where plastic strain rates $\dot{\bm{\varepsilon}}$ adhere to the associated flow rule defined by

[TABLE]

and $\dot{\mu}$ is a non-negative plastic multiplier. In addition, we define a convex set $\mathcal{B}$ which contains statically admissible stresses

[TABLE]

where $\sum$ is a space of symmetric stress tensors, $\bm{\sigma}$ satisfies the yield condition for assumed material and $\psi(\bm{\sigma})$ denotes the convex yield function.

Defining $\mathcal{C}=\{\mathbf{\dot{u}}\in\mathcal{V}|W_{int}(\mathbf{\dot{u}})=1\}$, the limit analysis problem is generally speaking to seek the limit load factor $\lambda^{*}$ yielding the following optimization problem [17]

[TABLE]

in which

[TABLE]

For plane strain problems, von Mises and Mohr–Coulomb models 111Our approach is also available for the Drucker–Prager model. Under plane strain conditions, the Drucker–Prager yield criterion can reduce to the Mohr–Coulomb yield criterion. can be coined in the following compact form [17]

[TABLE]

where

[TABLE]

and $(\varsigma=\sin\varphi,~{}\kappa=c\cos\varphi)$ is available for the Mohr–Coulomb criteria, c denotes the cohesion coefficient, I is the identity matrix, and $\varphi$ is the internal friction angle. The Mohr-Coulomb criteria reduces to the von Mises criteria as $\varsigma=0,\kappa=\sigma_{y}/\sqrt{3},$ in which $\sigma_{y}$ is the yield stress.

Under in the associated flow rule, the plastic dissipation can be expressed as a function of plastically admissible strains

[TABLE]

where

[TABLE]

and $\dot{e}_{m}=tr(\dot{\bm{\varepsilon}})=\bm{\nabla}\cdot\dot{\mathbf{u}}$ denotes the volume expansion rate with $\dot{e}_{ij}=\dot{\varepsilon}_{ij}-\dfrac{1}{2}\dot{e}_{m}\delta_{ij}$, and $J_{2}(\dot{\mathbf{e}})=\dfrac{1}{2}\dot{\mathbf{e}}\colon\dot{\mathbf{e}}$. As $\varsigma=0$, the plastic dissipation [11] is subjected to the incompressibility condition $\bm{\nabla}\cdot\dot{\mathbf{u}}=0$. Although lower-order finite element formulations are extremely convenient for computation, they are sensitive to the volumetric locking [7, 10, 12]. In the current study, we introduce a simple formulation based on quadrilateral elements with piecewise-linear shape functions defined over barycentric coordinates [64] and employ an area-averaged edge-based projection technique [65] to achieve the numerical stability in the plastic regime.

3 A finite element limit analysis formulation

3.1 A finite element space enriched with bubble piecewise-linear basis functions

We consider a bounded domain $\Omega$ discretized into a set $\mathcal{T}$(at primal mesh) of $\mathit{n_{e}}$ non-overlapping polygonal elements 222We use the concept of consider polygonal elements in order to emphasize that quadrilateral elements with presence of hanging nodes are defined as polygonal elements. involving a set $\partial\mathcal{T}$ of edges $\mathit{n_{ed}}$ and $\mathit{n_{n}}$ nodes such that $\Omega\approx\Omega^{h}=\bigcup^{n_{e}}_{e=1}\mathbf{\Omega^{e}}$. Let $\mathcal{V}^{h}\subset\mathcal{V}$ be a finite element approximation space of kinematically admissible velocity fields. The discrete form of limit analysis is to seek an approximately collapse load factor $\alpha^{+}$ as

[TABLE]

It turns out that lower-order finite elements are not stable for purely cohesive materials. They perform poorly under the incompressibility condition. Mathematically, the numerical solution must satisfy the following condition

[TABLE]

It is clear that such a constraint requires a special treatment to obtain a stable element formulation in plastic regime. Roughly speaking, the total number of degrees of freedom is not great enough to dominate the number of incompressible constraints in the FE discretization [7, 8, 9, 10]. Several advanced remedies have been devised to provide locking-free finite elements for limit analysis such as mixed finite elements [9, 12], discontinuity velocity fields [10, 13, 15], and linear strain elements [17]. We now describe the finite dimension space $\mathcal{V}^{h}$ as

[TABLE]

where $\mathcal{Q}(\Omega^{e})$ contains barycentric basis functions over the convex polygon $\Omega^{e}$. As shown in [7], the space $\mathcal{V}$ needs be enriched to overcome volumetric locking in the fully plasticity regime. To address this, we introduce a piecewise-linear bubble space $\mathcal{V}^{h}_{b}$ such that

[TABLE]

The finite element space for the velocity field which is enriched with such bubble functions is defined as

[TABLE]

It will be shown in numerical examples of plane strain problems that a following displacement-based FE formulation in combination with bubble functions and dual-mesh integration scheme can avoid volumetric locking in the plastic collapse analysis. We write explicitly the velocity $\mathbf{\dot{u}}^{h}\in\mathcal{W}^{h}$ over $\Omega^{e}\in\mathcal{T}$ as follows

[TABLE]

where $n_{nod}$ is the number of vertices of element, $\mathbf{I}_{2}$ denotes the unit matrix of $2^{nd}$ rank, $\dot{\mathbf{d}}^{e}_{i}$ is the vector of nodal degrees of freedom of $\dot{\mathbf{u}}^{h}(\mathbf{x})$ associated to the {\color[rgb]{0,0,1}i^{th}} vertex of the element, $N^{e}_{i}$ is the standard nodal basis shape function at the {\color[rgb]{0,0,1}i^{th}} vertex of element, $\dot{\mathbf{d}}^{e}_{b}$ is unknowns associated with the centroid node, and $N^{e}_{b}$ is the shape function at a centroid node of the element. In next section, we construct piecewise-linear shape functions for polygonal elements, which are used to define quadrilateral elements with the presence of hanging nodes.

3.2 Basis functions on arbitrary polygonal elements

It is recalled that construction of shape functions {\color[rgb]{0,0,1}N_{i}(\mathbf{x})} over arbitrary polygonal elements cannot be defined in a similar way as those for the standard triangular or quadrilateral elements. Some approaches to basis functions have been developed such as Wachspress [66, 67], Mean-value, Laplace [68] and piecewise-linear shape functions. A comparison of different shape functions of four regular polygonal elements is shown in Fig. 2. More simply, Floater et al. [64] proposed sharp upper and lower bound piecewise-linear functions (cf. Fig. 3(a) and (b)), which can be used for basis functions. These shape functions have the essential feature of interpolants such as non-negative, partition of unity, Kronecker delta, and linear precision. On the other hand, the Gaussian integration rule is valid for triangular and quadrilateral finite elements. For arbitrary polygonal elements, the Gaussian integration is often performed over sub-triangles. However, this integration technique leads to the huge increase of incompatibility constraints resulting in the poor accuracy of limit analysis. Moreover, we learn in limit analysis that no significant difference in solution using Wachspress, mean-value, Laplace and piecewise-linear shape functions, while rational basis functions of Wachspress, mean-value and Laplace coordinates raise high computational cost. Therefore, piecewise-linear shape functions are recommended. They are defined over the sub-triangles which are created by connecting the centroid of the polygonal element to the vertices of the element as illustrated in Fig.4.

We recall that the shape functions at the vertices of polygonal elements fulfill the Kronecker delta property:

[TABLE]

The shape functions at the centroid $\mathbf{x}_{c}$ are then evaluated by

[TABLE]

Now shape functions and their derivatives are accomplished by using linear shape functions over sub-triangles:

[TABLE]

where $N_{\triangle}^{l}(\mathbf{x})$ and $\bm{\nabla}N_{\triangle}^{l}(\mathbf{x})$ are linear shape functions and their derivatives on each sub-triangle $\Omega_{\triangle}$, $\phi_{i}^{e}(\mathbf{x}_{l})$ denotes the shape function of node i at the node l of $\Omega_{\triangle}$. In Eq. 22, shape functions are the piecewise-linear continuous functions over sub-triangles and have the following basic properties: (1) Kronecker delta $N_{i}^{e}(\mathbf{x}_{j})=\delta_{ij}$; (2) partition of unity $\sum_{i=1}^{n}N_{i}^{e}(\mathbf{x})=1$; (3) non-negative $N_{i}^{e}(\bar{\mathbf{x}})\geq 0$ and linear compatibility $\sum_{i=1}^{n}N_{i}^{e}(\mathbf{x})\mathbf{x}_{i}=\mathbf{x}$.

In addition, we introduce piecewise-linear bubble shape functions which help to enrich a space of the velocity field for solving incompressible constraints. On each sub-triangle $\Omega_{\triangle}$, we now construct a linear shape function using barycentric coordinates with vertices including $\bm{x}_{c}$ and two endpoints of the edges. In other words, $N^{e}_{b}$ is a linear function on $\Omega_{\triangle}$. A 3D view of the bubble functions is illustrated in Fig. 5. It is evident that the bubble shape function is zero along the elemental boundary and unity at the centroid. Furthermore, to facilitate computation of underlying quantities in elements with hanging nodes, we adopt the following concept of conforming polygonal element[44].

3.3 On polygon with side-nodes

Consider “polygonal” element $\Omega^{e}$ displayed in Fig. 6 (a), which has one side-node along the edge. To construct shape functions for a convex polygon with side nodes, we construct the shape functions on a polygonal reference element $\Omega^{\xi}$ with the reference coordinate $\mathbf{\xi}(\xi_{1},\xi_{2})\in\Omega^{\xi}$ [44]. As shown in the Fig. 6 (b), the nodes of the polygonal reference element are located at $\xi_{i}=\{\cos(2\pi i/n),\ \sin(2\pi i/n)\}$. The shape functions of polygon with side-nodes are then obtained by mapping from this polygonal reference element to the polygonal physical element. It should be emphasized that there is no distinction in using piecewise-linear interpolation between side-nodes and corner nodes. This means that all the nodes are regarded as corner nodes and the shape functions of all the nodes are obtained in a similar manner.

3.4 Quadrature scheme

It is known that the lower-order quadrilateral element (Q4) does not work well for nearly-incompressible and incompressible materials, and even cohesive-frictional materials. This shortcoming also happens with Q4 having presence of hanging nodes. As an alternative way mentioned in Section 3.2, piecewise-linear shape functions are defined over the sub-triangles of quadrilateral elements. To reduce incompressibility constraints, we employ an alternative integration technique of a so-called area-averaged edge-based projection technique, which follows the concept of stabilized conforming nodal integration (SCNI)[69] and edge-based smoothed finite element method (ES-FEM) [65]. As a result, the present method suppresses volumetric locking in purely cohesive materials and performs well for cohesive-frictional materials. To use efficiently this alternative integration scheme, a dual mesh is defined over edge-based mesh background $\bar{\mathcal{T}}$ by connecting vertices and centroids of elements of $\mathcal{T}$. Then, $\Omega$ is subdivided into $n_{s}$ edge-based domains $\Omega^{(k)}$ based on edges of elements such that $\Omega=\bigcup\!^{n_{s}}_{k=1}\Omega^{(k)}$ and $\Omega^{(i)}\cap\Omega^{(j)}=\varnothing$ for $i\neq j$. The edge-based integration domain $\Omega^{(k)}$ or smoothing domain associated with the edge k is obtained by connecting two end-nodes of the edge k to the centroids of adjacent elements, as depicted in Fig. 8.

It is now noticeable that in our study assumed strain rates need to be defined over an average projector $\mathcal{P}^{(k)}$ of compatible strain rates on a dual mesh [65]

[TABLE]

where $A^{(k)}$ is the area of the edge-based integration domain $\Omega^{(k)}$ and $\bm{\nabla}_{s}$ denotes a matrix of differential operators defined as

[TABLE]

We write

[TABLE]

In addition, the method must satisfy the incompressibility condition for the purely cohesive material:

[TABLE]

Substituting Eq. 19 into Eq. 26, the approximation of the constant strain rates on can be expressed as

[TABLE]

where $n_{nk}$ is the number of the neighboring nodes of edge $k$, $n_{b}$ is the number of central nodes or bubble nodes in $\Omega^{(k)}$, $\dot{\bar{\mathbf{d}}}^{(k)}_{I}$ is the nodal degrees of freedom at the $I^{th}$ node of domain $\Omega^{(k)}$, and $\bar{\mathbf{B}}^{(k)}_{I}$ is the area-averaged edge-based strain-displacement matrix on $\Omega^{(k)}$ computed by

[TABLE]

where

[TABLE]

Eq. 29 shows that the computations are performed on sub-domain $\Omega^{(k)}_{i}$ of $\Omega^{(k)}$. To achieve numerical integration, the polygonal reference element $\Omega^{\xi}$ introduced in Section 3.3 is divided into sub-triangles and the Gaussian integration is performed over sub-triangle $\Omega^{\xi}_{\triangle}$ of the reference element $\Omega^{\xi}$ corresponding to sub-domain $\Omega^{(k)}_{i}$ of the physical element $\Omega^{e}$, see Fig. 9. The piecewise-linear shape functions and their derivatives at all vertices evaluated at the Gaussian point $\bm{\xi}_{\triangle m}$ in the sub-triangle $\Omega^{\xi}_{\triangle}$ are computed as follows

[TABLE]

The numerical integration in 30 becomes:

[TABLE]

where $n_{g}$ indicates the number of Gauss points, $\nabla_{s}N^{e}_{I}(\mathbf{x}_{\triangle Gm})$ denotes a matrix of shape function derivatives at Gauss point $\mathbf{x}_{\triangle Gm}$ in the physical coordinates and $w_{\triangle m}^{\xi},w_{\triangle m}^{\eta}$ represents the corresponding Gaussian quadrature weights defined on polygonal reference element and on three-node triangle natural element, respectively. From the Jacobian matrix $\mathbf{J}_{\xi}$, the shape function derivatives $\bm{\nabla}N^{e}_{I}(\mathbf{x}_{\triangle Gm})$ in the physical coordinates are determined via the relation: $\bm{\nabla}N^{e}_{I}(\mathbf{x}_{\triangle Gm})=\mathbf{J}^{-1}_{\xi}\bm{\nabla}N_{I}^{e}(\bm{\xi}_{\triangle m})$, in which $\bm{\nabla}N_{I}^{e}(\bm{\xi}_{\triangle m})$ are shape function derivatives which are computed at Gauss points $\bm{\xi}_{\triangle m}$ in the polygonal reference element $\Omega^{\xi}$.

3.5 Bubble-enriched finite element limit analysis

The discrete problem 14 is rewritten in the following form

[TABLE]

where

[TABLE]

It is evident that the flow rule can be satisfied at every point over $\Omega^{(k)}$, and likewise everywhere in a set $\bar{\mathcal{T}}$ of the problem domain. On the other hand, the edge-based averaging strains rates defined on a dual mesh somewhat debilitates the compatibility due to the smoothed strain rate which is constant over the smoothing domain. Hence, the present method cannot produce properly a strict upper bound, but the collapse load factor found is reliable and applicable.

4 Solution procedure of the discrete problem

The aforementioned limit analysis problem can be solved by a general non-linear optimization solver [9]. Furthermore, it can be rewritten in a sum of norms [11] in second-order cone programming (SOCP), which is solved by the interior point method [62, 63]. To address this, the limit analysis problem with $n_{s}$ constraints can be formed as follows

[TABLE]

where $t_{k}\in\mathcal{R}^{+},k=1,...,n_{s}$ or $\mathbf{t}\in\mathbf{R}^{n_{s}}$ denotes optimization variables and the coefficients are $c_{k}\in R,\mathbf{R}_{k}\in R^{m_{dim}\times n_{s}},\mathbf{v}_{k}\in R^{m_{dim}},\mathbf{y}_{k}\in R^{n_{s}}$ and $z_{k}\in R$. The SOCP problem is used for $m_{dim}$ = 2 or $m_{dim}$ = 3. When $m_{dim}$ = 1, the SOCP problem reduces to a linear programming problem. Now we represent the limit analysis problem in the form of the second-order quadratic cone as

[TABLE]

The limit analysis problem given in 34 can be reduced to the problem of minimizing a sum of norms [11]. Now the optimization problem can be written in a common way of minimizing a sum of norms as follows

[TABLE]

where $\mathbf{\rho}^{k}$ denotes auxiliary variables. The optimization problem in 39 is classified into a second-order cone programming (SOCP) with the constraint of quadratic cones in Eq. 34. As a result, we can use the fast solver of the Mosek optimization package for this problem. The robust feature of Mosek is to provide high computational efficiency and accuracy for solving the large-scale optimization problems in practice.

5 On an adaptive quadtree mesh refinement

Since uniform meshes are computationally expensive to solve plastic collapse problems, adaptive mesh refinements are very useful for adjusting mesh resolution in order to improve the solution with the desired accuracy. In the limit analysis, the solution accuracy depends on the fast detection of the yield lines or plastic collapse zones over the problem domain. This means that a mesh should be refined along the yield lines while it is gradually coarse away from the yield lines. Regarding computational efficient, the mesh size should be adjusted and refined automatically or adaptively [27]. Theoretically, an effective error estimator or indicator to guide the mesh refinement process is required yet it is available for elasticity. Nevertheless, in plastic limit analysis, finding such an error indicator is challenging. Therefore, we in this study introduce an alternative indicator based on the $\mathcal{L}^{2}$ - norm measure of quantities in plastic collapse zones based on plastic strain rates. One of our important contributions in the present method is highly capable of generating the automatic local refinement with a quadtree mesh. Ideally, the quadtree mesh is a recursive spatial decomposition of the quadrilateral element into four new child elements (children). It is worth emphasizing that such a mesh procedure is performed easily together with the proposed indicator.

5.1 $\mathcal{L}^{2}$ - norm-based indicator

We know that the presence of localized plastic deformations or high strain rates at limit state very significantly affects the accuracy of numerical solutions. Using uniform meshes leads to huge computational cost to gain the limit load factor of the desired accuracy due to refined meshes at unnecessary regions. Therefore, adaptive meshes should be done so that the approximated velocity/strains rates converge the actual ones with the lowest computational cost as possible. To achieve this, we introduce an alternative indicator in adaptive refinement procedure based on the $\mathcal{L}^{2}$ - norm of plastic strain rates, which can be defined over each element as

[TABLE]

where $n^{e}_{edge}$ indicates the number of elemental sides and {\color[rgb]{0,0,1}\Theta^{(k)}} is the $\mathcal{L}^{2}$ - norm of plastic strain rates on side k of the element $e$ and is calculated by Eq. 5 as

[TABLE]

5.2 Refinement strategy

Once plastic strain rates are found, it is used to evaluate the corresponding $\mathcal{L}^{2}$ - norm and to guide the mesh refinements. The global indicator, $\Theta$, can be described as the contribution of the local refinement indicator, $\Theta^{e}$, for all the individual elements

[TABLE]

We furthermore mark a set $\mathcal{M}\subseteq\mathcal{T}$ of $\{\Omega^{e}\}\in\mathcal{T}$ such that

[TABLE]

A new mesh $\mathcal{T}^{{}^{\prime}}$ is then derived from $\mathcal{T}$ by refining the marked elements $\Omega^{e}\in\mathcal{M}$ by using quadtree meshes.

5.3 Description of quadtree meshes

The quadtree meshes were introduced in [30], where a spatial decomposition structure is used for quad mesh refinements. In the quadtree mesh, each father element is partitioned recursively into four child quadrants. In Fig. 1, a quadtree mesh and its representative tree are presented. After each refinement, hanging nodes are created if new elements and their neighbours are not at same levels of refinement. Hanging nodes consist of the vertices of the child element that is placed on the side of the father element (e.g., diamond-shape nodes in Fig. 1 (a)). The presence of these hanging nodes leads to the incompatibilities in finite element approximations. Several methods have been devised to treat these issues such as Lagrange multipliers or penalty method [38], adding temporary elements [39], constraining hanging nodes to corner nodes [40], hierarchical enrichment [41], B-splines [42], natural neighbor basis functions [44]. In [44], authors constructed natural neighbor (Laplace) basis functions to obtain $C^{0}$ continuous approximations along sides containing hanging nodes. These shape functions were defined over the polygonal reference elements with an affine map. Hence, the number of hanging nodes on each side is arbitrary as shown in several available techniques [38, 39, 40, 41, 42]. In this study, we utilize $C^{0}$ admissible approximations along sides by using a mapping from the polygonal reference element but we employ the barycentric coordinates to construct piecewise-linear shape functions of quadrilaterals with hanging nodes. Note that the classical Wachspress coordinates are not easy to construct the shape functions of polygon with side-nodes (including hanging nodes) on the quadtree meshes. As described in the Section 3.3, we can construct shape function for quadrilaterals with side-nodes using a mapping $\mathbf{J}_{\xi}=\partial\mathbf{x}/\partial{\xi}$. Note that hanging nodes and vertices are now regarded as vertices. By this way, conforming shape functions for quadtree meshes with the number of hanging nodes can be well established.

6 Numerical implementation

In this section, we summarize algorithms for the implementation of the present method. Three main algorithms are given as

6.1 Algorithm 1: Main program

Define the problem domain (coordinate, load and boundary conditions): $\Omega,\Gamma_{\dot{\mathbf{u}}},\Gamma_{\mathbf{t}}$.

• -

Give number of iterations, $n_{iter}$, and refinement factor $\theta$.

• -

Given plastic parameters: $c$, $\varphi$.

Create $\mathbf{\Lambda}=\displaystyle\sum_{e=1}^{n_{e}}\{\Lambda\}^{e}$ % element connectivity matrix

Create $\mathbf{\Lambda}_{p}=\mathbf{\Lambda}$ % element connectivity matrix on quadtree

for i = 1 $\rightarrow n_{iter}$

  Add boundary conditions $\Gamma_{\dot{\mathbf{u}}},\Gamma_{\mathbf{t}}$.

  $\mathbf{\Lambda}\rightarrow$ add one node at the centroid of each element $\rightarrow\mathbf{\Lambda}_{c}$

  Create edge-based integration domains $\Omega=\bigcup\!^{n_{s}}_{k=1}\Omega^{(k)}$

  for k = 1 $\rightarrow n_{s}$ % $n_{s}$ *is the number of edge-domains

 *   Call algorithm 2: compute strain matrix using $\mathbf{\Lambda}_{c}$

  **end for

 ** Compute auxiliary variables and define SOCP

  Call MOSEK: optimization solution

  Evaluate element-based strain-norm indicator $\bm{\eta}=\displaystyle\sum_{e=1}^{n_{e}}\eta^{e}$

  Call algorithm 3: adaptation on quadtree meshes

end for

6.2 Algorithm 2: Compute strain matrix

Input: $\mathbf{\Lambda}_{c}$, information of $\Omega^{(k)}$

Give $\bar{\mathbf{B}}^{(k)}=0$

   for ic = 1 $\rightarrow$ number of neighbor cell that k edge influence

   + Find $\mathbf{x}_{ic}$ of element $\Omega^{(k)}_{ic}$ contain the edge domain *k

   +* Define reference element $\Omega^{\xi}$ and find ${\bm{\xi}}$ $\rightarrow$ Find Gauss point $\bm{\xi}_{\triangle}$ and its weight $w_{\triangle m}^{\xi}$ of sub-triangle $\triangle$ corresponding to element $\Omega^{(k)}_{ic}$

   + Define function $\phi^{e}_{I}$ for all vertex of element $\rightarrow$ Find $\phi^{\triangle}_{I}$ corresponding to element $\Omega^{(k)}_{ic}$

     for m = 1 {\color[rgb]{0,0,1}\rightarrow} $n_{g}$ % where $n_{g}$ is number of Gauss point in sub-triangle element $\triangle$

      + Compute $N_{b}(\bm{\xi}_{\triangle m}),\bm{\nabla}N_{b}(\bm{\xi}_{\triangle m})$

      + Define $N_{\triangle}(\bm{\xi}_{\triangle m}),\bm{\nabla}N_{\triangle}(\bm{\xi}_{\triangle m})$

      + Compute linear shape function and their derivatives at all vertices:

      ** ** N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}=\displaystyle\sum\limits_{l=1}^{3}N^{l}_{\triangle}(\bm{\xi}_{\triangle m})\phi^{\triangle}_{I} and \bm{\nabla}N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}=\displaystyle\sum\limits_{l=1}^{3}\bm{\nabla}N^{l}_{\triangle}(\bm{\xi}_{\triangle m})\phi^{\triangle}_{I}

      + Compute linear shape function and their derivatives including bubble node:

       N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}=\big{[}N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}\ N_{b}(\bm{\xi}_{\triangle m})\big{]} and \bm{\nabla}N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}=\big{[}\bm{\nabla}N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}\ \bm{\nabla}N_{b}(\bm{\xi}_{\triangle m})\big{]}

      + Compute Jacobian: $\mathbf{J}_{\xi}=\dfrac{\partial\mathbf{x}_{ic}}{\partial\xi}$

      + Compute \bm{\nabla}N_{I}^{e}\big{(}\mathbf{x}_{\triangle Gm}\big{)}=\mathbf{J}^{-1}_{\xi}\bm{\nabla}N_{I}^{e}\big{(}\bm{\xi}_{\triangle m}\big{)}

        B^{*}_{m}=B^{*}_{m}+\bm{\nabla}N_{I}^{e}\big{(}\mathbf{x}_{\triangle Gm}\big{)}\big{|}\mathbf{J}_{\xi}\big{|}w_{\triangle m}^{\xi}

     end for % end loop Gauss point of sub-triangle $\triangle$

     $\bar{\mathbf{B}}^{(k)}=\bar{\mathbf{B}}^{(k)}+\mathbf{B}^{*}_{i}$

   **end for

  ** $\bar{\mathbf{B}}^{(k)}=\dfrac{1}{{A}^{(k)}}\bar{\mathbf{B}}^{(k)}$

   Output $\bar{\mathbf{B}}^{(k)}$

6.3 Algorithm 3: Adaptation on quadtree meshes

Input:$\theta,\Theta,\bm{\Lambda},\bm{\Lambda}_{p}$

   Arrange strain-norm indicator $\Theta$ from largest to smallest

  + Mark elements will be adaption

   $\sum\limits_{e=1}^{n_{e}}\Theta^{e}\geq\theta\Theta$

   while   $\sum_{e=1}^{n_{e}}\Theta^{e}\geq\theta\Theta$

    $\tau^{e}$ = 1  % marking the elements need to be divided

   **end while

  +** Divide parent element into child elements and add new element data

   for e = 1 $\rightarrow n_{e}$

      if $\tau^{e}$ = 1 **then

      +** Divide the parent element into 4 child elements

      + Add 4 new child elements data to $\bm{\Lambda}_{p}$: \bm{\Lambda}_{p}=\bm{\Lambda}_{p}+\big{\{}\Lambda_{p}\big{\}}^{e}_{new} and \bm{\Lambda}=\bm{\Lambda}+\big{\{}\Lambda_{p}\big{\}}^{e}_{new}

      + Delete parent element

      **end if

  ** **end for

  +** Find hanging nodes

   for e = 1 $\rightarrow n_{e}$

      for iedg = 1 $\rightarrow n_{edg}$ % where $n_{edg}$ is the number of edges of element

     + Find nodes $\mathbf{x}_{new}$ on edge

      if $\mathbf{x}_{new}\neq\varnothing$ **then

      ** \big{\{}\Lambda_{p}\big{\}}^{e}=\big{\{}\Lambda_{p}\big{\}}^{e}+\mathbf{x}_{new}\%\emph{ update element data information}\\ \text{\qquad\qquad\qquad}\rightarrow\mathbf{\Lambda}

      **end if

     ** end for % end loop on edges of elements

   **end for

  ** Output: $\bm{\Lambda},\bm{\Lambda}_{p}$

7 Numerical validations

In this section, we examine the performance of the proposed approach through three benchmark problems. The program is complied on Macbook Air (Intel Core I7, 2.0GHz CPU, 8G RAM). The conic interior-point optimizer of the academic MOSEK package is employed.

7.1 Strip footing bearing capacity

The first example considered in our analysis is a smooth strip footing of width $B$ on a weightless cohesive–frictional soil ($c\geq 0$, $\varphi\geq 0$) as given in Fig. 10a. We investigate the collapse plastic load for the structure subjected to loading $Q$. The analytical solution of the limit load factor for the first case ($\varphi=0$) is $\alpha^{*}=2+\pi$ and for the second case is given as $\alpha^{*}=[e^{\pi\tan\varphi}\tan^{2}(\pi/4+\varphi/2)-1]\cot\varphi$ [70]. Applying the symmetry boundary condition, a half of the foundation with length $L$ and height $H$ is considered.

For purely cohesive soil case ($\varphi=0$), the problem is modelled with $L=2.5B$ and $H=B$. A uniform mesh of $40$ quadrilateral (Q4) elements is indicated in Fig.10b. Some adaptive meshes are depicted in Fig.11. Fig. 12 exhibits the convergence of the limit load factor using uniform and adaptive meshes. Table 1 shows the convergence of solution using adaptive meshes. As expected, the accuracy of solution is improved a lot after using adaptive meshes. Regarding computational efficiency, to gain a limit load factor value $\alpha^{+}=5.149$ for bL–Q4, optimization Mosek time takes around 2.7s for 16542 optimization variables. For comparison, Table 2 lists the present result and several other ones. Again, we show that the present method is strongly competitive to other researches [10, 17]. Taking 16542 variables, the present method using adaptive bL–Q4 meshes yields ultra-accurate solution with a small error of 0.14% in comparison with the error of 1.33% reported in [10]. Note that the approach in [10] based on the discontinuous velocity field can leads to the poor performance for unstructured meshes [17]. Makrodimopoulos and Martin [17] introduced a simplex strain element formulation for unstructured meshes. Here, it is numerically proved that our approach works well for unstructured meshes. It is seen that the present limit load value matched well with a strict upper bound solution proposed by Makrodimopoulos and Martin [17] using 18719 6-node triangular meshes (approximately 149752 variables). Also, the present solutions are more accurate than that of mixed FE formulation[9] and are acceptable with that of lower bound formulation [16]. Fig. 16 plots plastic dissipation distribution of a purely cohesive soil problem using an adaptive mesh. Last but not least, the convergence of relative error is shown in Fig. 13. As we can seen from this figure, the present method has a good convergence in comparison with uniform mesh (bL–Q4).

For cohesive–frictional material case ($\varphi=35^{\circ}$), the problem is modeled with $L=13B$ and $H=10B$. Several adaptive meshes are displayed in Fig. 14. The approximate upper bound value against adaptive meshes is listed in Table 1. As expected, the collapse load factor $\alpha^{+}=46.292$ using bL–Q4 for the last adaptive mesh has a quite small error of $0.36\%$ while optimization Mosek time takes approximately $8$s with $43967$ variables. The efficiency of adaptive refinement is that more than $90$ percent of the elements are located in the plastic localization zone. For comparison, Table 2 shows that the present method is very good competitor to the linear strain element [17] using $18719$ 6-node triangular meshes (approximately $149752$ variables).

Additionally, a series of limit load factors are evaluated by both the present and analytical approaches for various internal frictional angles. For the present approach, domain size and initial meshes are the same for all angles and the limit load factors are taken at the eighth step of the successive adaptive mesh refinement loops. The final adaptive meshes are different for each internal frictional angle which lead to the discrepancies in the number of optimization variables. These variables tend to increase correspondingly with the expansion of the internal frictional angles and the final values are $150442$, $204437$, $256457$, $369947$, $449972$, $649247$, $921087$, $1295417$, $1456077$, $1502162$ for the angles of $0^{\circ},5^{\circ},10^{\circ},15^{\circ},20^{\circ},25^{\circ},30^{\circ},35^{\circ},40^{\circ},45^{\circ}$, respectively. These evaluations are then put side by side in Table 3 as well as illustrated in Fig. 17 for comparison purpose. As it pointed out in this table, the relative error values are smaller than $0.1\%$ in all cases except for the last one which is mainly attributed to the computational domain modelled is not large enough to capture the behaviour of plastic dissipation at large frictional angles (see Fig. 18). These analyses once again demonstrate the high accuracy and flexibility of the present method.

7.2 Block with two symmetric holes

The second benchmark is a rectangular domain with two small holes of same radius which has a geometry and applied load as shown in Fig. 19. This problem also investigated by Zouain et al.[12], Makrodimopoulos et al.[16], Munoz et al.[28], H. Nguyen-Xuan et al.[71]. The upper bound solutions of this problem was reported in [17], namely $\alpha^{+}$ = 1.825 for $\varphi$ = 0*∘* and $\alpha^{+}$ = 1.063 for $\varphi$ = 30*∘* using 177426 stress variables and in [71], namely $\alpha^{+}$ = 1.817 for $\varphi$ = 0*∘* using 34948 stress variables and $\alpha^{+}$ = 1.061 for $\varphi$ = 30*∘* using 29768 stress variables are used as reference upper bound values. Also, the lower bound solutions given in [16], namely $\alpha^{-}$ = 1.8089 for $\varphi$ = 0*∘* and $\alpha^{-}$ = 1.0562 for $\varphi$ = 30*∘* and in [28] , namely $\alpha^{-}$ = 1.8119 for $\varphi$ = 0*∘* and $\alpha^{-}$ = 1.0581 for $\varphi$ = 30*∘. All computations are then performed over adaptive mesh refinements. Fig. 20 and Fig. 21 plot several steps of an adaptive mesh process with respect to the internal friction angle $\varphi$ = 0∘* and $\varphi$ = 30*∘, respectively. It is seen that the present method is well captured in accurately capture the plastic zones with adaptive meshes. The present solution versus adaptive meshes is indicated in Table 4. Fig. 22 and Fig. 23 exhibit the limit load factors using uniform and adaptive meshes. Present results agree well a strict upper bound value using the linear strain elements [17]. It is observed that the adaptive bL-Q4 element produces the most accurate solutions for this problem. The global errors are less than 1% for a slightly coarse mesh of $N_{var}$ = 6028 for $\varphi$ = 0∘* and $N_{var}$ = 5713 for $\varphi$ = 30*∘. These do not exceed 0.2% of the global errors for slightly fine meshes with variables around $N_{var}$ = 56818 for $\varphi$ = 0∘* and $N_{var}$ = 67458 for $\varphi$ = 30*∘*. Additionally, a comparison for bL-Q4 listed in Table 5 demonstrates good performance of the present method. Also, Fig. 24 present accurately the representation of plastic collapse regions which is confirmed by other studies [17, 28].

7.3 Slope stability

The next example is a homogeneous slope of cohesive-frictional soil with inclination 70*∘* and height H. The problem model and boundary conditions are given in Fig. 25. The load in this problem is only the soil weight, $\gamma$. The limit load factor or the stability factor is defined as $\alpha^{s}=\gamma H/c$. The analytical solution (upper bound) for the homogeneous slope stability was investigated by Chen[72]. Some authors have also studied numerical solutions by the upper bound models (see the works of Krabbenhoft et al.[13], Makrodimopoulos et al.[17], Lyamin et al.[15]) or the lower bound models (see the works of Makrodimopoulos et al.[16]) for the slope stability in past decades.

In this example, we consider $c=1,H=1$ and two internal friction angles ($\varphi=20^{\circ}$ and $\varphi=35^{\circ}$) in our analyses. Here we demonstrate the high performance of the present approach for unstructured meshes. Several adaptive mesh steps for internal friction angles $\varphi=20^{\circ}$ and $\varphi=35^{\circ}$ are shown in Fig 26 and Fig. 27, respectively. The results of analyses are presented in Table 6 and Fig. 28 for $\varphi=20^{\circ}$ and Fig. 29 for $\varphi=35^{\circ}$. The present solutions agree well with available ones in the literature as listed in Table 7. In detail, obtained solutions are between the upper bound value [17] and lower bound value [16]. Fig.30 displays plastic dissipation distribution. Such plastic collapse zones in structure is detected well by our approach.

7.4 Failure in porous media

Final example is a square specimen measures $L=10$ in length and $H=10$ in width which contains $16$ equal holes measure $r_{c}=0.5$ in radius. The distances between holes’ centers in both direction equal to $l=2$. The applied load and boundary conditions are shown in Fig 31. This example was investigated by Nguyen-Xuan and Liu [19] with aims to analyze the failure mechanism in materials with cavities and showed the applicability of the limit analysis approach to homogenization of strength properties in multiscale modeling of porous materials. Fig 32 and Fig 34 depict several meshes for detecting failure mechanisms using adaptive procedure for two cases of internal friction angles $\varphi=0^{\circ}$ and $\varphi=20^{\circ}$, respectively. The plastic dissipation distribution of the present approach is provided in Figs 33 and 35. The present solution is in good agreement with the result reported in Ref [19]. Furthermore, our approach is capable of producing properly failure mechanisms under several frictional angles for this problem.

8 Conclusion

We presented an efficient and fast numerical scheme over adaptive quadtree mesh generation for limit analysis of structures. The inherent difficulty of the presence of hanging nodes during quadtree mesh generation is simply solved by its definition over arbitrary polygons. The approximate velocity field used piecewise-linear shape functions defined on the primal mesh, while the plastic dissipation was computed over the dual mesh based on the edge-based integration background. The area-averaged edge-based projection was adopted to redistribute the plastic strain rates on the edge-based integration domains. As expected, the flow condition is satisfied at everywhere in the whole structure. The limit analysis model converted into the form of the second-order cone programming (SOCP) was significantly solved by interior-point solvers for the large-scale optimization problems. A quadtree-based adaptive refinement approach controlled by the $\mathcal{L}^{2}$ -norm-based indicator of plastic strain rates showed the high effectiveness in plastic collapse analysis of structures. The numerical results demonstrated excellent agreement with the reference solutions and compared very well with several other numerical ones. The present method is very simple yet efficient to implement into the existing limit analysis packages and enables us to further extend to other plastic models. Finally, another interesting feature of this study will be extensible to a more general form of yield criteria.

Acknowledgements

The author would like to thank Mr. Son Nguyen-Hoang for his comments and suggestions to this work. The financial support provided by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2016.19, and RISE-project BESTOFRAC (734370) is gratefully acknowledged.