Adaptive selective ES-FEM limit analysis of cracked plane-strain structures

Abstract

This paper presents a simple and efficient approach for predicting the plastic limit loads in cracked plane-strain structures. We use two levels of mesh repartitioning for the finite element limit analysis. The master level handles an adaptive primal-mesh process through a dissipation-based indicator. The slave level performs the subdivision of each triangle into three sub-triangles and constitutes a dual mesh from a pair of two adjacent sub-triangles shared by common edges of the primal mesh. Applying a strain smoothing projection to the strain rates on the dual mesh, the incompressibility constraint and the flow rule constraint are imposed over the edge-based smoothing domains and everywhere in the problem domain. The limit analysis problem is recast into the compact form of a second-order cone programming (SOCP) for the purpose of exploiting interior-point solvers. The present method retains a low number of optimization variables. It offers a convenient way for designing and solving the large-scale optimization problems effectively. Several benchmark examples are given to show the simplicity and effectiveness of the present method.

Graphical abstract

Keywords

cracked structure / limit analysis / von Mises criterion / edge-based strain smoothing / second-order cone programming / adaptive

1 Introduction

Limit analysis furnishes a direct tool for evaluating elastic−plastic fracture toughness and safety of fracture failure structures. The slip-line field (SLF) theory has commonly been used for evaluating analytically the load bearing capacity of cracked structures [1−3]. However, it heavily depends on geometry and loading conditions of the underlying problem. The earlier approaches for estimating the limit collapse load rely on the upper and lower bound theorems. The Koiter’s kinematic (upper bound) theorem [4] uses a kinematically admissible displacement field yielding the minimum load factor, while the lower bound theorem (or Melan’s static theorem) [5] employs a statically admissible stress field to determine the maximum load factor. It is well known that the analytical methods following these two approaches are unsuitable for general problems in engineering practice [6,7]. Finite element approaches for the limit load of cracked structures have therefore been studied, e.g., [8−11]. Alternative formulations include meshfree methods [12], smoothed finite elements [13,14] and isogeometric analysis (IGA) [15,16]. It was shown in Ref. [16] that higher-order methods produce sufficiently accurate solutions, but an adaptive mesh refinement of such higher-order elements is more challenge.

In this study, we focus on lower-order finite elements based on triangular meshes owing to its simplicity. Due to the enforcement of the incompressibility phenomenon in cracked plane-strain problems, these conventional FE models often work poorly. Hence various advanced techniques have been developed [17−26]. Among kinematic formulations based on lower-order finite elements, the approach based on discontinuity velocity fields [18,22,24] is popular in the finite element limit analysis. The discontinuous velocity [4] in constant finite elements is a simple and efficient approach, but it performs poorly for unstructured meshes [4]. The linear strain finite elements (or quadratic velocity fields) [26] can be acted as an alternative way to deal with certain shortcomings of constant strain elements.

On the other hand, the localized plastic deformations in the limit analysis lead to the slow convergence of the numerical solutions [27]. Therefore, the adaptive mesh refinement should be enhanced. Due to the absence of a priori error estimate, posteriori approaches are therefore the most suitable way to derive the adaptive mesh refinement. The local and global errors based on a posteriori error estimate have been studied intensively. A directional error estimate with the recovery gradients and/or Hessian of a finite element solution was studied in Refs. [27,28]. The other effective error estimate based on the global error of solving the combined solutions of both the lower and upper bound problems can be found in Refs. [29,30]. A dissipation-based mesh adaptivity to accomplish closely bracketed lower and upper bounds of the exact collapse load was recently investigated in Ref. [31]. As a special case of a dissipation-based mesh adaptivity [31], the indicator based on plastic dissipation was also used in the adaptive mesh refinement scheme [32−34]. Plastic dissipation indicator permits to detect plastic regions or slip-line pattern where the adaptive re-meshing is needed. In addition, an adaptivity technique in meshfree methods [35] is also promising for limit analysis problems. Very recently, an adaptive selective ES-FEM, which used a two-level mesh repartitioning strategy [36] in combination with a projection operator of strain rates, for plastic collapse analysis has been proposed in Ref. [37].

In this paper, we exploit an adaptive selective ES-FEM for the finite element limit analysis of cracked plane-strain structures. The two-level mesh repartitioning scheme in this study consists of two mesh levels. The kinematic formulation relies on the usual piecewise linear displacements and is constructed based on a primal mesh or master level. A projection operator [38] of strain rates performed over a dual mesh (slave level) of common edges of adjacent elements is used to fulfill the flow rule constraint and to offer an optimal value of constraint ratio for the volumetric locking-free computation. The resulting non-smooth optimization problem is then recast in the form of minimizing a sum of Euclidean norms with an aim of using an efficient interior-point method. More importantly, an adaptive procedure is added at the master-level mesh repartitioning to enhance the computational efficiency. The method is very simple to implement into a finite element code. It will be shown in numerical examples that our method produces a high-precision solution with a low number of degrees of freedom.

The rest of the paper is outlined as follows: the kinematic theorem is reviewed in the next section. Section 3 presents a two-level mesh finite element formulation and an edge-based strain smoothing. Section 4 summarizes a solution procedure of the discrete problem. An adaptive mesh refinement is given in Section 5. The numerical procedure is presented in Section 6. Section 7 illustrates three numerical examples. Finally, we close with some concluding remarks.

2 A background on the kinematic theorem

Let us consider a rigid-perfectly cracked plastic body of a two-dimensional problem domain bounded by a boundary of continuous and discontinuous parts \mathrm{\Gamma }={\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\cup {\mathrm{\Gamma }}_{\mathbf{t}}\cup {\mathrm{\Gamma }}_c,{\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\cap {\mathrm{\Gamma }}_{\mathbf{t}}\cap {\mathrm{\Gamma }}_c=\varnothing, where {\mathrm{\Gamma }}_c stands for the crack boundary. The body force f acts within the domain and the surface traction g is enforced on the free portion {\mathrm{\Gamma }}_{\mathbf{t}} while the boundary {\mathrm{\Gamma }}_{\dot{\mathbf{u}}} is prescribed by the velocity vector \dot{\mathbf{u}}. The weak form can be stated as follows

(1)a\left(\mathbf{\sigma },\dot{\mathbf{u}}\right)=f\left(\dot{\mathbf{u}}\right),\forall \dot{\mathbf{u}}\in V,

where the external work rate reads

(2)f(\dot{\mathbf{u}})={\displaystyle {\int }_{\mathrm{\Omega }}f\cdot \dot{\mathbf{u}}}\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathbf{t}}}g\cdot \dot{\mathbf{u}}}\mathrm{d}\mathrm{\Gamma },

and the internal work rate is

(3)a\left(\mathbf{\sigma },\dot{\mathbf{\epsilon }}(\dot{\mathbf{u}})\right)={\displaystyle {\int }_{\mathrm{\Omega }}\mathbf{\sigma }}\dot{\mathbf{\epsilon }}(\dot{\mathbf{u}})\mathrm{d}\mathrm{\Omega }.

We define a space of kinematically admissible velocity fields indicated by

(4)\mathrm{V}=\left\{\dot{\mathbf{u}}\in {({\mathrm{H}}^{}(\mathrm{\Omega }))}^{\mathrm{2}},\dot{\mathbf{u}}\mathbf{=}\dot{\bar{\mathbf{u}}}\hspace{0.17em}\mathrm{on}\hspace{0.17em}{\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\right\},

and a convex set \mathrm{S} which contains admissible stress fields such that

(5)\mathrm{S}=\{\mathbf{\sigma }\in \mathrm{\Sigma }|\psi (\mathbf{\sigma })\le \mathrm{0}\},

where the yield function \psi (\mathbf{\sigma }) is convex, \mathbf{\sigma } satisfies the yield condition and \mathrm{\Sigma } denotes a space of symmetric stress tensors.

The limit analysis problem is to seek the actual collapse multiplier {\alpha }^{*} such that [20]

(6)\begin{array}{l}{\alpha }^{*}=\max \left\{\exists \mathbf{\sigma }\in S|a\left(\mathbf{\sigma },\mathbf{ϵ}\left(\dot{\mathbf{u}}\right)\right)=\alpha f\left(\dot{\mathbf{u}}\right),\forall \dot{\mathbf{u}}\in V\right\}\hfill \\ =\underset{\sigma \in S}{\max }\underset{\dot{\mathbf{u}}\in W}{\min }a\left(\mathbf{\sigma },\dot{\mathbf{ϵ}}\left(\dot{\mathbf{u}}\right)\right)\hfill \\ =\underset{\dot{\mathbf{u}}\in W}{\min }\underset{\mathbf{\sigma }\in S}{\max }f\left(\mathbf{\sigma },\dot{\mathbf{ϵ}}\left(\dot{\mathbf{u}}\right)\right)\hfill \\ =\underset{\dot{\mathbf{u}}\in W}{\min }D\left(\dot{\mathbf{u}}\right),\hfill \end{array}

where the set \mathrm{W} is defined by \mathrm{W}=\left\{\dot{\mathbf{u}}\in \mathrm{V}|f(\dot{\mathbf{u}})=\mathrm{1}\right\} and the dissipation is D(\dot{\mathbf{u}})=\underset{\mathbf{\sigma }\in \mathrm{S}}{\max }a\left(\mathbf{\sigma },\dot{\mathbf{\epsilon }}(\dot{\mathbf{u}})\right).

Regarding the von Mises criteria presented for plane strain cracked structures, one writes

(7)\psi (\mathbf{\sigma })=\sqrt{J_{\mathrm{2}}(\mathbf{s})}-k\le \mathrm{0},

where

(8)J_{\mathrm{2}}(\mathbf{s})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{s}:\mathbf{s},{\sigma }_m={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{tr}(\mathbf{\sigma }),\mathbf{s}=\mathbf{\sigma }-{\sigma }_m\mathbf{I},k={\displaystyle \frac{{\sigma }_y}{\sqrt{\mathrm{3}}}},

and {\sigma }_y is the yield stress.

Now the explicit form of the plastic dissipation can be formulated as a function of the plastically admissible strains (those obeying the associated flow rule) [18]:

(9)D(\dot{\mathbf{u}})={\displaystyle {\int }_{\mathrm{\Omega }}\mathrm{2}k\sqrt{J_{\mathrm{2}}(\dot{\mathbf{e}})}}d\mathrm{\Omega },

where {\dot{e}}_{ij}={\dot{\epsilon }}_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\nabla \cdot \dot{\mathbf{u}}{\delta }_{ij}, and J_{\mathrm{2}}(\dot{\mathbf{e}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\dot{\mathbf{e}}:\dot{\mathbf{e}}.

Theoretically, the incompressibility condition \nabla \cdot \dot{\mathbf{u}}=\mathrm{0} is required to ensure the finite plastic dissipation [20]. Numerically, this causes difficulties to lower-order finite elements [17,18]. For this reason, we introduce an alternative formulation based on the edge-based strain smoothing technique [38] and the two-level mesh repartitioning scheme [36].

3 On a selective edge-based strain smoothing

3.1 A two-level mesh repartitioning scheme for the lower-order finite element formulation

Assume that a bounded domain \mathrm{\Omega } is discretized into a set \mathfrak{T}(at master level) of N_e triangular elements having a set \partial \mathfrak{T} of edges N_s and N_n nodes such that \mathrm{\Omega }\approx {\mathrm{\Omega }}^h={\displaystyle {\sum }_{e=\mathrm{1}}^{N_e}{\mathrm{\Omega }}_e}. Each element {\mathrm{\Omega }}_e is then implemented with one node at its centroid and is further divided into three sub-triangles (or macro-elements) as shown in Fig. 1. Let us denote a set \mathfrak{M}(at slave level) of \mathrm{3}\times N_e macro-elements involving a set \partial \mathfrak{M} of N_s+\mathrm{3}\times N_e edges and N_n+N_e nodes such that \mathrm{\Omega }\approx {\mathrm{\Omega }}^h={\displaystyle {\sum }_{e=\mathrm{1}}^{\mathrm{3}\times N_e}{\mathrm{\Omega }}_e^m} as illustrated in Fig. 2. Let {\mathrm{V}}^h\subset \mathrm{V} be a finite element approximation space of kinematically admissible velocity fields. The approximate displacement space {\mathrm{V}}^h involves piecewise linear functions denoted as

(10){\mathrm{V}}^h=\left\{{\dot{\mathbf{u}}}^h\in {({\mathrm{H}}^{\mathrm{1}}(\mathrm{\Omega }))}^{\mathrm{2}},{{\dot{\mathbf{u}}}^h|}_{{\mathrm{\Omega }}_e}\in {\mathrm{P}}_{\mathrm{1}}{({\mathrm{\Omega }}_e^m)}^{\mathrm{2}}\right\},

where {\mathrm{P}}_{\mathrm{1}}({\mathrm{\Omega }}_e^m) stands for the set of polynomials of degree 1 for each displacement component.

The discrete form is to seek a collapse multiplier {\alpha }^{+} that satisfies the following optimization problem:

(11){\alpha }^{+}=\min D^h\left({\dot{\mathbf{u}}}^h\right),s\mathrm{.t}.\{\begin{array}{l}\nabla \cdot {\dot{\mathbf{u}}}^h=\mathrm{0},\\ {\mathbf{u}}^h\in {\mathrm{W}}^h\mathrm{.}\end{array}

It is evident from finite element limit analysis that constant strain elements lead to poor numerical performance in the incompressibility limit due to volumetric locking under plastic yield conditions [1−18]. Among various advanced numerical techniques, constant strain elements enhanced by discontinuous velocity fields [18,22,24] and the linear strain elements [26] are very efficient. It was proven in Ref. [36] that a {\mathrm{V}}^h space defined in Eq. (10). can be combined with the two-level mesh repartitioning scheme so that such constant strain elements pass the inf-sup condition [39] which is also required in the plastic limit analysis. In the following, we generalize the two-level mesh repartitioning scheme in the selective ES-FEM for accurately computing the plastic limit loads of plane-strain fracture problems.

3.2 Selective edge-based smoothed strain rates

The basic idea is to construct a dual mesh for imposing the incompressibility condition [36]. Therefore, we first define a slave mesh \mathfrak{M} of \mathrm{3}\times N_e macro-elements with N_n+N_e nodes and N_s+\mathrm{3}\times N_e edges. A dual mesh \bar{\mathfrak{T}} is then obtained based on smoothing domains {\mathrm{\Omega }}_k^s with the area A_k^s. Each smoothing domain is created by embracing a pair of macro-elements having a common edge k\in \{\mathrm{1},\mathrm{2},\cdots ,N_s\} such that \mathrm{\Omega }={\displaystyle \underset{k=1}{\overset{N_s}{\cup }}{\mathrm{\Omega }}_k^s} and {\mathrm{\Omega }}_i^s\cap {\mathrm{\Omega }}_j^s=\varnothing, i\ne j as shown in Fig. 2. More details for the mesh repartitioning scheme can be founded in Ref. [36].

Next we reconstruct constant strain rates via a smoothing or constant (average) projector P_k^h of the compatible strain rates on a dual mesh of edge-based smoothing domains [38]:

(12){\dot{\bar{\epsilon }}}^{(k)}=\stackrel{‾}{{\mathbf{\nabla }}_s{\dot{\mathbf{u}}}^h}=P_k^h{\mathbf{\nabla }}_s{\dot{\mathbf{u}}}^h={\displaystyle \frac{\mathrm{1}}{A_k^s}}{\displaystyle {\int }_{{\mathrm{\Omega }}_k^s}{\mathbf{\nabla }}_s{\dot{\mathbf{u}}}^{h}d\mathrm{\Omega }},

where {\mathbf{\nabla }}_s is the matrix of differential operators.

For the smoothed volume expansion rate, one has

(13)\stackrel{‾}{\nabla \cdot {\dot{\mathbf{u}}}^h}=P_k^h\nabla \cdot {\dot{\mathbf{u}}}^h={\displaystyle \frac{\mathrm{1}}{A_k^s}}{\displaystyle {\int }_{{\mathrm{\Omega }}_k^s}\nabla \cdot {\dot{\mathbf{u}}}^{h}d\mathrm{\Omega }}.

Fulfilling the incompressibility condition requires

(14)\forall {\mathrm{\Omega }}_k^s\in \bar{\mathfrak{T}},\stackrel{‾}{\nabla \cdot {\dot{\mathbf{u}}}^h}=\mathrm{0.}

It was shown in Ref. [36] that volumetric locking can be avoided by enforcing the incompressibility condition on the dual mesh through the two-level mesh repartitioning scheme.

3.3 Limit load statement for sES-FEM

The discrete problem is to seek an approximated collapse multiplier {\alpha }_{}^{+}

{\alpha }^{+}=\min {\bar{D}}^h\left(\Vert {\overline{\dot{\mathbf{e}}}}_s^h\left({\dot{\mathbf{u}}}^h\right)\Vert \right),

such that

(15)\{\begin{array}{l}\sqrt{J_{\mathrm{2}}\left({\overline{\dot{\mathit{e}}}}^h\right)}=\Vert {\overline{\dot{\mathit{e}}}}_s^h\Vert \hfill \\ \nabla \cdot {\dot{\mathbf{u}}}^h=\mathrm{0}\hfill \\ {\dot{\mathbf{u}}}^h\in W^h,\hfill \end{array}

where

(16){\overline{\dot{\mathit{e}}}}_s^h={\left[\begin{array}{ll}\mathrm{2}{\dot{e}}_{\mathrm{11}}^h\hfill & \mathrm{2}{\dot{e}}_{\mathrm{12}}^h\hfill \end{array}\right]}^{\mathrm{T}}\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}{\overline{D}}^h\left({\dot{\mathbf{u}}}^h\right)={\displaystyle {\int }_{\mathrm{\Omega }}k\Vert {\overline{\dot{\mathit{e}}}}_s^h\Vert }\mathrm{d}\mathrm{\Omega }={\displaystyle \sum _{k=\mathrm{1}}^{{\mathit{{\rm N}}}_s}k{A}_k^s\Vert {\overline{\dot{\mathit{e}}}}_{sk}^h}.

Due to the smoothed strain rates being constant over the smoothing domain {\mathrm{\Omega }}_k^s, the incompressibility constraints is well enforced over each smoothing domain, and is hence fulfilled over the entire problem domain.

4 The discrete problem with the second-order cone programming

The limit analysis problem given in Eq. (15) is a nonlinear optimization problem [19] with equality constraints. As stated in Ref. [20], some problems can be solved by minimizing a sum of norms. Hence the underlying optimization problem can be rewritten as follows:

{\alpha }^{+}=\min {\displaystyle \sum _{k=\mathrm{1}}^{N_s}k{A}_k^s{t}_k},

such that

(17)\{\begin{array}{l}{t}_k\ge \Vert {\overline{\dot{\mathit{e}}}}_{sk}^h\Vert ,\hfill \\ P_k^h\nabla \cdot {\dot{\mathbf{u}}}^h=\mathrm{0},\hspace{0.17em}k=\mathrm{1},\cdots ,N_s,\hfill \\ {\dot{\mathbf{u}}}^h={\dot{\mathbf{u}}}_{\mathrm{0}}\hspace{0.17em}\mathrm{on}\hspace{0.17em}{\mathrm{\Gamma }}_{\dot{\mathbf{u}}},\hfill \\ f\left({\dot{\mathbf{u}}}^h\right)=\mathrm{1}.\hfill \end{array}

The optimization problem Eq. (17) belongs to a second-order cone programming (SOCP) form [40,41]. Let N_var be the total number of variables of the optimization problem. For plane strain problems, N_var equals NoDofs + 3N_s where NoDofs is the total number of the degrees of freedom (DOFs) of the discrete problem.

5 An adaptive mesh refinement

In plastic collapse analysis, the plastic strain rate varies substantially over the problem domain, especially along the yield lines. Adaptive h-refinement is on ideal choice to reduce the high computational cost. Ideally, a mesh should be fine along the yield lines and gradually coarse when it is away from the plastic zones. More importantly, this mesh size adjustment should be done automatically, as illustrated in Ref. [29]. In elasticity, an effective error estimator or indicator is often used to quantify the error arising from numerical methods. In plastic collapse analysis, such an error indicator is a non-trivial task. Alternatively, an indicator based on the plastic dissipation through the element can be used for the activation of adaptive re-meshing.

In the present method, the performance of the automatic local refinement is quite simple and efficient. Having this advantage is due to utilizing three-node triangular elements.

5.1 Dissipation-based indicator

The localized plastic deformations heavily affect the accuracy of numerical solutions in limit analysis. At limit state, the plastic dissipation will concentrate in the regions characterized by high strain rates. We here use the plastic dissipation of each element regarding as an indicator in the adaptive h-refinement procedure. It is calculated as

(18){\eta }_e={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}({\eta }^{(i)}+{\eta }^{(j)}+{\eta }^{(k)}),

where {\eta }^{(i)},{\eta }^{(j)} and {\eta }^{(k)} are the plastic dissipations on the three edges i, j and k of the element e, respectively, and the plastic dissipation on edge k is given as

(19){\eta }^{(k)}=k{A}_k^s{t}_k.

5.2 Plastic dissipation indicator and refinement strategy

The global indicator, \eta, can be computed as the sum of the local refinement indicator for all the individual elements:

(20)\eta ={\displaystyle \sum _{e=\mathrm{1}}^{nel}{\eta }_e}.

Through the local refinement indicator, {\eta }_e, one marks the elements {\mathrm{\Omega }}_e\in \mathfrak{T} for refinement following the Dorfler criterion [42]. A minimal set \mathfrak{M}\subseteq \mathfrak{T} is determined such that

(21){\displaystyle \sum _{{\mathrm{\Omega }}_e\in \mathfrak{M}}{\eta }_e}\ge \vartheta \eta ,\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{some}\hspace{0.17em}\vartheta \in (\mathrm{0},\mathrm{1}).

A new mesh \mathfrak{T} is then generated from \mathfrak{T} by refining (at least) the marked elements {\mathrm{\Omega }}_e\in \mathfrak{M} though the algorithm of newest vertex bisection reported in Refs. [43,44].

6 Numerical procedure

In this section, we describe the numerical implementation of the present method. The algorithm is summarized as follows:

1) Define the problem domain;

2)Discretize the problem domain with the coarse primal mesh based on 3-node triangular elements and:

- Add one node at center of each element subdividing into 3 sub-triangles;

- Create element connectivity of sub-triangles (subsection 3.1).

3)Loop over each adaptive step // (iter= 1:MaxIt)

{

- Construct a dual mesh of edge-based smoothing domains (subsection 3.2);

- Loop over edge-based smoothing domains

{

+ Evaluate smoothed strain rates defined in Eq. (12);

+ Store smoothed strain rates for each flow rule point k.

}// End the loop over the smoothing domains

Enforce the incompressibility constraints in Eq. (14);

Introduce auxiliary variables and define a second-order cone programming (SOCP) in Eq. (17);

Exploit several SOCP solvers, e.g., MOSEK, for optimization solution;

Evaluate element-based dissipation indicator {\eta }_e and the global indicator (\eta)

While condition Eq. (21)

{

Call adaptive procedures in Section 5

}

- Add one node at center of each element and subdivide each element into three sub-triangles;

}// End the loop over final adaptive step

7 Numerical examples

This section aims to examine the performance of the present method through three benchmark problems. For comparison, the elements used in this paper are denoted as follows:

T3 − the linear triangular elements.

ES-T3 − the edge-based smoothed finite linear triangular elements [38,46].

sES-T3 − the selective edge-based smoothed finite linear triangular elements.

The program is compiled by a desktop computer with Intel® Xeon (2.4GHz CPU, 32G RAM). The conic interior-point optimizer of the academic MOSEK package [45] is employed.

7.1 Notched tensile specimen

The first example is a double notched tensile specimen problem [20] as shown in Fig. 3. Two typical meshes are plotted in Figs. 4 (a) and (b). Only the upper-right quarter of the specimen is modeled. The exact solution is not available, and a reference upper bound value of α* = 1.1316 has been provided by Ref. [20].

To prove the performance of the present method in the incompressibility regime, we also consider T3 and ES-T3 elements. The computation of limit loads for plane strain problems using ES-T3 has been already pointed out in Ref. [46]. However, the previous method required a special mesh arrangement [17,18]. As numerically validated below, it can lead to poor performance with unstructured meshes under the incompressibility constraints. In other words, the original ES-T3 [46] is not capable of producing accurate solutions for arbitrary meshes. Fig. 5 shows that the solution using both T3 and ES-T3 does not converge to the reference value with a primal mesh as given in Fig. 4(a). As expected, the sES-T3 works well.

We next evaluate the computational cost based on the total variables N_var, iterations, and optimization Mosek times. The results are displayed in Fig. 6. As seen, the present method solves well for the large-scale optimization problem with the maximum variables around N_var = 140000 after 19 step iterations. Although the total number of variables is quickly increased for fine meshes, the computing time is about 16 s only.

Now we use the adaptive re-meshing procedure provided in Section 5 to enhance the solution accuracy with a low number of variables. The mesh refinement process begins at a uniformly coarse mesh as plotted in Fig. 4(b). Several adaptive meshes are presented in Fig. 7. The mesh is refined primarily in the regions characterized by high strain rates. Fig. 8 shows the error in limit load factor with respect to the number of variables. The h-adaptive method reduces significantly the number of variables and reaches accurate solutions. Figure 9 demonstrates that the plastic dissipation distribution is produced well by the adaptive sES-T3. We also solved this problem for various crack lengths. As shown in Fig. 10, the present solutions match well with the reference values.

7.2 Centre cracked plate subjected to tension

Next we consider a plate of dimension b× H having a center crack of length a as shown in Fig. 11(a). The plate is subjected to tension. The analytical limit load multiplier is given as Ref. [3]

(22)\alpha ={\displaystyle \frac{\mathrm{2}}{\sqrt{\mathrm{3}}}}(\mathrm{1}-a/b).

Due to the symmetry of the problem, the upper-right quarter of the plate is modeled. Similar to the previous example, we verify the performance of the method for the case a/b = 0.5. A coarse mesh and several h-adaptive meshes are shown in Figs. 11(b)−(f). Limit load factor using uniform and adaptive refined meshes, respectively, are depicted in Fig. 12. As expected, the adaptive sES-T3 gains α⁺ = 0.581 (with the error of 0.7%) with a low number of 25312 variables. On the other hand, the collapse limit load value α⁺ = 0.583 (with the error of 1%) requires a uniformly fine mesh of 166532 variables. It is worth to mention that our approach employs only the constant strain rate and nonetheless produces accurate results for the unstructured meshes. The plastic dissipation distribution in the full model using an adaptive meshing algorithm is displayed in Fig. 13.

The present method is also applied for the problem with various crack lengths. The results are shown in Fig. 14. An excellent agreement to the analytical solutions is observed.

7.3 Single-edge cracked plate subjected to tension

The last example is a plate of length H width b and a single edge cracked of length a. It is subjected to tension as depicted in Fig. 15(a). The analytical solution with respect to x=a/b was reported in Ref. [2] as

(23)\alpha ={\displaystyle \frac{\mathrm{3.404}}{\sqrt{\mathrm{3}}}}\{{[{(\mathrm{2.06}-x)}^{\mathrm{2}}+\mathrm{0.5876}{(\mathrm{1}-x)}^{\mathrm{2}}]}^{\mathrm{1}/\mathrm{2}}+(\mathrm{2.06}-x)\},\forall x\ge \mathrm{0.545},

and in case of , the solution is bounded by Ref. [3]

(24){\displaystyle \frac{\mathrm{2}}{\sqrt{\mathrm{3}}}}(\mathrm{1}-x-\mathrm{1.232}{x}^{\mathrm{2}}+x^{\mathrm{3}})\le \alpha \le {\displaystyle \frac{\mathrm{2}}{\sqrt{\mathrm{3}}}}(\mathrm{1}-x-\mathrm{1.232}{x}^{\mathrm{2}}+x^{\mathrm{3}}+\mathrm{22}{x}^{\mathrm{3}}{(\mathrm{0.545}-x)}^{\mathrm{2}}).

Due to its symmetry, only the upper haft of the plate is modeled. A coarse mesh is given in Fig. 15(b). For x=\mathrm{0.5}, the exact value of the limit load factor is bounded by 0.366≤α≤0.373. For comparison, the lower bound solution, namely α⁻ = 0.366, is used to measure the global error (%)

(25)error(\%)={\displaystyle \frac{{\alpha }^{+}-{\alpha }^{-}}{{\alpha }^{+}+{\alpha }^{-}}}\times \mathrm{100.}

Several adaptive re-meshing steps are given in Fig. 15(c)−(f). The mesh refinement is primarily performed along plastic zones. Fig. 16 plots the load factor versus the number of variables. As seen, the computed value using the adaptive mesh is quite close to the actual value. The error is less than 4% for N_var = 1087, while it requires a large number of variables (N_var = 16502) for the uniform mesh to reach an error below 4%. To reduce the error to 0.8%, the present method requires a low number of variables (N_var = 27467) for the adaptive refined mesh, while for a uniform mesh with N_var = 166532, the error is approximately 1.4%. Figure 17 shows plastic dissipation zones which agrees well with results in Ref. [10]. Finally, Fig. 18 shows the load factor for various crack lengths.

8 Conclusions

We have presented a constant strain element formulation and an adaptive h-refinement for plastic collapse analysis of plane-strain fracture structures. The present method displaced the constant strain rate over a dual mesh by the two-level mesh repartitioning scheme. The average operator over edge-based strain smoothing domains of the dual mesh was used to enforce the incompressibility condition, which is then held everywhere in the entire structure. The method is capable of handling a large size of optimization problem, which is formulated by the second-order cone programming (SOCP) and is solved efficiently by interior-point solvers. The adaptive mesh procedure was based on the dissipation indicator and the latest vertex bisection strategy was sufficiently exploited in this study. Several numerical examples were provided and the obtained results showed high accuracy and effectiveness of the present method. To end the paper, the present method can be straightforwardly extended to the 3D fracture problems for the limit and shakedown analysis.

References

[1] Hill R. On discontinuous plastic states, with special reference to localized necking in thin sheets. Journal of the Mechanics and Physics of Solids, 1952, 1(1): 19–30

[2] Ewing D J F, Richards C E. The yield-point loading of singly-notched pin loaded tensile strips. Journal of the Mechanics and Physics of Solids, 1974, 22(1): 27–36

[3] Miller A G. Review of limit loading of structures containing defects. International Journal of Pressure Vessels and Piping, 1988, 32(1−4): 197–327

[4] Koiter W T. General theorems for elastic plastic solids. Progress in Solid Mechanics, Sneddon I N, Hill R, eds. Nord-Holland, Amsterdam, 1960, 1: 165–221

[5] Melan E. Theorie statisch unbestimmter Systeme aus ideal plastischem Baustoff. Sitzber. Akad. Wiss. Wien IIa, 1936, 145: 195–2182

[6] Prager W, Hodge PGJr. Theory of Perfectly Plastic Solids. New York: Wiley, 1951, 3

[7] Chakrabarty J. Theory of Plasticity. 3rd ed. Elsevier Butterworth-Heinemann, 2006, 4

[8] Yan A M, Nguyen-Dang H. Limit analysis of cracked structures by mathematical programming and finite element technique. Computational Mechanics, 1999, 24(5): 319–333

[9] Vu D K. Dual Limit and Shakedown analysis of structures. Dissertation for the Doctoral Degree. Belgium: Université de Liège, 2001, 5

[10] Khan I A, Ghosh A K. A modified upper bound approach to limit analysis for plane strain deeply cracked specimens. International Journal of Solids and Structures, 2007, 44(10): 3114–3135

[11] Khan I A, Bhasin V, Chattopadhyay J, Singh R K, Vaze K K, Ghosh A K. An insight of the structure of stress fields for stationary crack in strength mismatch weld under plane strain mode−I loading−Part II: Compact tension and middle tension specimens. International Journal of Mechanical Sciences, 2014, 87: 281–296

[12] Le C V, Askes H, Gilbert M. A locking-free stabilized kinematic EFG model for plane strain limit analysis. Computers & Structures, 2012, 106−107: 1–8

[13] Tran T N, Liu G R, Nguyen-Xuan H, Nguyen-Thoi T. An edge-based smoothed finite element method for primal-dual shakedown analysis of structures. International Journal for Numerical Methods in Engineering, 2010, 82(7): 917–9386

[14] Nguyen-Xuan H, Rabczuk T, Nguyen-Thoi T, Tran T N, Nguyen-Thanh N. Computation of limit and shakedown loads using a node-based smoothed finite element method. International Journal for Numerical Methods in Engineering, 2012, 90(3): 287–310

[15] Nguyen-Xuan H, Thai C H, Bleyer J, Nguyen P V. Upper bound limit analysis of plates using a rotation-free isogeometric approach. Asia Pacific Journal on Computational Engineering, 2014, 1(1): 12

[16] Nguyen-Xuan H, Tran L V, Thai C H, Le C V. Plastic collapse analysis of cracked structures using extended isogeometric elements and second-order cone programming. Theoretical and Applied Fracture Mechanics, 2014, 72: 13–27

[17] Nagtegaal J C, Parks D M, Rice J R. On numerically accurate finite element solutions in the fully plastic range. Computer Methods in Applied Mechanics and Engineering, 1974, 4(2): 153–177

[18] Sloan S W, Kleeman P W. Upper bound limit analysis using discontinuous velocity fields. Computer Methods in Applied Mechanics and Engineering, 1995, 127(1−4): 293–314

[19] Capsoni A, Corradi L. A finite element formulation of the rigid-plastic limit analysis problem. International Journal for Numerical Methods in Engineering, 1997, 40(11): 2063–2086

[20] Christiansen E, Andersen K D. Computation of collapse states with von Mises type yield condition. International Journal for Numerical Methods in Engineering, 1999, 46(8): 1185–1202

[21] Vu D K, Yan A M, Nguyen-Dang H. A primal-dual algorithm for shakedown analysis of structure. Computer Methods in Applied Mechanics and Engineering, 2004, 193(42−44): 4663–4674

[22] Krabbenhøft K, Lyamin A V, Hjiaj M, Sloan S W. A new discontinuous upper bound limit analysis formulation. International Journal for Numerical Methods in Engineering, 2005, 63(7): 1069–1088

[23] Vicente da Silva M, Antao A N. A non-linear programming method approach for upper bound limit analysis. International Journal for Numerical Methods in Engineering, 2007, 72(10): 1192–1218

[24] Lyamin A V, Sloan S W. Upper bound limit analysis using linear finite elements and nonlinear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 2002, 26(2): 181–216

[25] Makrodimopoulos A, Martin C M. Lower bound limit analysis of cohesive-frictional materials using second-order cone programming. International Journal for Numerical Methods in Engineering, 2006, 66(4): 604–634

[26] Makrodimopoulos A, Martin C M. Upper bound limit analysis using simplex strain elements and second-order cone programming. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(6): 835–865

[27] Borges L A, Zouain N, Costa C, Feijoo R. An adaptive approach to limit analysis. International Journal of Solids and Structures, 2001, 38(10−13): 1707–1720

[28] Lyamin A V, Sloan S W, Krabbenhoft K, Hjiaj M. Lower bound limit analysis with adaptive remeshing. International Journal for Numerical Methods in Engineering, 2005, 63(14): 1961–1974

[29] Ciria H, Peraire J, Bonet J. Mesh adaptive computation of upper and lower bounds in limit analysis. International Journal for Numerical Methods in Engineering, 2008, 75(8): 899–944

[30] Munoz J, Bonet J, Huerta A, Peraire J. Upper and lower bounds in limit analysis: adaptive meshing strategies and discontinuous loading. International Journal for Numerical Methods in Engineering, 2009, 77(4): 471–501

[31] Martin C M. The use of adaptive finite-element limit analysis to reveal slip-line fields. Géotechnique Letters, 2011, 1(April-June): 23–29

[32] Van-Phuc P, Nguyen-Thoi T, Nguyen C H, Le V C. An effective adaptive limit analysis of soil using FEM and second-order cone programming. The International Conference on Advances in Computational Mechanics (ACOME), 2012, 177–190

[33] Le V C. A stabilized discrete shear gap finite element for adaptive limit analysis of Mindlin−Reissner plates. International Journal for Numerical Methods in Engineering, 2013, 96: 231–246

[34] Nguyen-Xuan H, Liu G R. An edge-based finite element method (ES-FEM) with adaptive scaled-bubble functions for plane strain limit analysis. Computer Methods in Applied Mechanics and Engineering, 2015, 285: 877–905

[35] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[36] Wu C T, Hu W. A two-level mesh repartitioning scheme for the displacement-based lower-order finite element methods in volumetric locking-free analyses. Computational Mechanics, 2012, 50(1): 1–18

[37] Nguyen-Xuan H, Wu C T, Liu G R. An adaptive selective ES-FEM for plastic collapse analysis. European Journal of Mechanics- A/Solid, submitted, 2014,

[38] Liu G R, Nguyen-Thoi T, Lam K Y. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration, 2009, 320(4−5): 1100–1130

[39] Chapelle D, Bathe K J. The inf-sup test. Computers & Structures, 1993, 47(4−5): 537–545

[40] Andersen K D, Christiansen E, Conn A R, Overton M L. An efficient primal-dual interior-point method for minimizing a sum of Euclidean norms. SIAM Journal on Scientific Computing, 2001, 22(1): 243–262

[41] Andersen E D, Roos C, Terlaky T. On implementing a primal-dual interior-point method for conic quadratic programming. Mathematical Programming, 2003, 95(2): 249–277

[42] Dorfler W. A convergent adaptive algorithm for Poisson’s equation. SIAM Journal on Numerical Analysis, 1996, 33(3): 1106–1124

[43] Rivara M C, Venere M. Cost analysis of the longest-side (triangle bisection) refinement algorithms for triangulations. Engineering with Computers, 1996, 12(3-4): 224–234

[44] Funken S, Praetorius D, Wissgott P. Efficient implementation of adaptive p1-FEM in Matlab. Preprint, 2008 (

[45] Mosek. The MOSEK optimization toolbox for MATLAB manual. Mosek ApS, Version 5.0 Edition.9, 2009

[46] Le C V, Nguyen-Xuan H, Askes H, Rabczuk T, Nguyen-Thoi T. Computation of limit load using edge-based smoothed finite element method and second-order cone programming. International Journal of Computational Methods, 2013, 10(1): 1340004

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg