Polytopal composite finite elements

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 355 (2019) 405–437 www.elsevier.com/locate/cma

Polytopal composite finite elements H. Nguyen-Xuan ∗, Khanh N. Chau, Khai N. Chau CIRTech Institute, Ho Chi Minh City University of Technology (HUTECH), Ho Chi Minh City, Viet Nam Received 11 December 2018; received in revised form 18 June 2019; accepted 20 June 2019 Available online xxxx

Highlights • • • • •

We propose a polytopal composite finite element approach for solid mechanics problems. The method relies on a polynomial projection of strain field performed over macro-elements. The present method passes a linear patch test and verifies the inf–sup condition for nearly incompressible materials. It works well with arbitrary polygonal/polyhedral meshes. An applicable capability of the proposed method to industrial problems is confirmed.

Abstract In recent years, polygonal and polyhedral finite elements have gained much interest from researchers as an alternative approach for discretizing complicated domains. Findings to enhance their performance are still desirable to solve a broad class of mechanics problems. In this study, we propose a new formulation based on arbitrary polytopes, which we name as Polytopal Composite finite Element (PCE). The key idea is to devise a polynomial projection of compatible strain fields through the least-squares approximation. For nearly incompressible problems, we design a pair of projection operators of volumetric and deviatoric strains resulting in stability of pressure solution. The present approach fulfills a patch test and satisfies the inf–sup stability. Through several numerical investigations, we show that the proposed method reaches the theoretical convergence rate and significantly improves the accuracy of polygonal element based solutions. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Finite elements; Polytope; Elasticity; Incompressibility; Patch test; Inf–sup

1. Introduction 1.1. State-of-the-art review Polytopes, which are polygons in a two-dimensional (2D) domain and polyhedrons for a three-dimensional (3D) domain, are often found extensively in nature and real life such as giraffe markings, turtle shell, honeycomb, pineapple, fish scales, etc. For example, we see several illustrations of polytopes as shown in Fig. 1. In recent decades, polytopal (i.e., polygonal/polyhedral) finite elements have experienced rapid development with ∗ Corresponding author.

E-mail addresses: [email protected] (H. Nguyen-Xuan), [email protected] (K.N. Chau), [email protected] (K.N. Chau). https://doi.org/10.1016/j.cma.2019.06.030 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

406

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 1. Illustration of several polytopal shapes in nature.

applications to a vast range of mechanics problems. They provide more flexibility in meshing procedure for arbitrary geometries as well as mitigate mesh-distortion sensitivity encountered in conventional FEM. In the state-of-the-art development of advanced finite elements, polytopal finite element methods (PFEM) have been applied to a variety of problems such as micromechanical analysis [1,2], second-order elliptic boundary-value problems [3,4], topology optimization [5,6], incompressible fluid flow [7,8], fracture modeling [9–11], finite elastic deformations [12–15], plates [16], limit analysis [17] and so on. Here we focus on several relevant contributions related to our current study. More details on polytopal finite elements can be found in the recent review [18]. Regarding the displacementbased polytopal finite element approaches, non-polynomial basis functions [19–21] are mostly used. However, due to non-polynomial nature of basis functions, direct application of traditional numerical integration methods leads to certain integration errors when evaluating the weak form [3]. The popular integration scheme is to subdivide the polygonal/polyhedral element into sub-triangles/sub-tetrahedrons (or sub-cells) followed by applying the standard triangular quadrature rules on these sub-cells. Unfortunately, this procedure is sub-optimal and does not exactly produce the patch test [22]. Of course, the improvement can be achieved if a very large number of integration points is employed but that consequently requires high computational expense, especially for 3D problems and thus is not applicable. Various advanced techniques have been devised in recent years to alleviate quadrature errors in numerical integration of the weak form such as Schwarz–Christoffel conformal mapping [23], generalized Gaussian quadrature rules [24], constant strain smoothing [25], harmonic shape functions [4], homogeneous function theorem [26,27], nonconforming basis functions [28,29], gradient correction [30,31], scaled boundary polygon [10,32], linear smoothing [33], virtual elements [34,35] and the combined PFEM–VEM approach [36]. Also, Astaneh et al. [37] proposed ultraweak variational formulations for high-order polygonal discontinuous Petrov—Galerkin (PolyDPG) methods. Very recently, a fast numerical integration on polytopal meshes has been presented for discontinuous Galerkin finite element methods [38]. According to authors’ knowledge, only some of the aforementioned methods are able to pass the patch test to machine precision. Furthermore, some of them are too complicated in implementation and are difficult or impossible to generalize for 3D problems. In addition, there are a few methods addressing a well-posed and stable formulation for incompressible media. 1.2. Related work Our approach relates to composite finite elements, each of which is commonly subdivided into macro-elements and approximations are then performed over these macro-elements in order to achieve a desired accuracy of solution. The early approach addressed the establishment of the C 1 -conforming composite triangle element (HCT) [39] or composite quadrilateral element (CQ) [40] for plate bending analysis. It was also used to formulate the quadrilateral statically admissible finite element [41]. Hackbusch and Sauter [42] formulated composite finite elements for modeling complicated micro-structures. Guo et al. [43] proposed six-node composite triangular elements assembled from triangular/quadrilateral micro-elements for 2D compressible and nearly incompressible elasticity. This approach was then extended to 3D solids [44,45]. A composite mixed triangular finite element [46],

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

407

in which the displacement field is approximated by a quadratic interpolation and the stress field is a piece-wise constant description through micro-elements of the triangular element. Dai et al. [25] applied a constant strain smoothing to micro-elements or smoothing domains in the polygonal finite element method (nSFEM) for solving solid problems. Francis et al. [33] developed linear strain smoothing technique to overcome shortcomings of the constant strain smoothing used in nSFEM, which does not pass the patch test when rational shape functions are employed. 1.3. Objective In this paper, we construct a new computational scheme on polytopal finite element formulations through a polynomial projection of compatible strains being derived from rational basis functions. This projection technique allows us to retain moderate number of quadrature points over sub-cells, without loss of accuracy of solutions when evaluating the weak form. For nearly incompressible problems, we propose a pair of polynomial projections such that the stability of pressure is ensured. Our numerical tests show that the present method satisfies the patch test and accomplishes the inf–sup condition. Through numerical examples, we show excellent performance of the proposed approach in comparison with several existing methods. 1.4. Outline The paper is outlined as follows: The statement of the polytopal finite element method is briefed in Section 2. Next, a polytopal composite finite element formulation is presented in Section 3. Then, Section 4 coins numerical implementation. Subsequently, numerical examples are described in Section 5. Finally, Section 6 closes some concluding remarks. 2. Brief on the finite elements with arbitrary polytopes 2.1. Problem definition and weak form Let us consider a d-dimensional (d = 2 or 3) elastic domain Ω in Rd bounded by a Lipschitz continuous boundary Γ = Γu ∪ Γt , Γu ∩ Γt = ∅. The linear elastic body is acted by body forces b ∈ L2 (Ω )d and external tractions t ∈ H1/2 (∂Ω )d , and the boundary Γu is prescribed by the displacement field u. The relations between the displacement, strain and stress fields can be found in the literature [47]. We now find u = [u 1 , . . . , u d ]T ∈ V which satisfies the weak form of equilibrium equations as follows a(u, v) = ⟨ℓ, v⟩, ∀v ∈ V0

(1)

in which a(., .) is a bilinear operator with the internal work described by ∫ ∫ a(u, v) = λ εv (u)εv (v)dΩ + 2µ εT (u)Dµ ε(v)dΩ Ω

(2)

Ω

where ε v (u) = ∇ · u, λ, µ denote Lam´e constants, Dµ is a 6-by-6 identity matrix with 3 last entries on the main diagonal line replaced by 0.5, and ℓ is a linear functional with the external work defined by ∫ ∫ ⟨ℓ, v⟩ = b · vdΩ + t · vdΓ (3) Ω

Γt

in which V and V0 are spaces of admissible displacement fields defined by V = {u ∈ H1 (Ω )d , u = u¯ on V0 = {v ∈ H (Ω ) , v = 0 1

d

Γu }

(4)

on Γu }

(5)

In addition, we also denote a space S of arbitrary strains or stresses, and a subspace Sc of compatible strains which are defined by S = {ϵ, σ ∈ L2 (Ω )3(d−1) }, Sc = {ε ∈ S|ϵ = ∇ s u, u ∈ V} where ∇ s denotes the symmetric differential operator.

(6)

408

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

2.2. Finite element approximation We consider a bounded domain Ω ∈⋃Rd discretized into a set T of n e non-overlapping polytopal elements with e a set of n n nodes such that Ω ≈ Ω h = ne=1 Ωe . Let V h ⊂ V be a finite element approximation space of compatible displacement fields based on barycentric coordinates. The discrete form of an elastic problem is to find uh ∈ V h such that a(uh , v) = ⟨ℓ, v⟩, ∀v ∈ V0h

(7)

It was proved that Eq. (7) has a unique solution with respect to a mesh size h. Numerical integration of Eq. (7) has been reported in many researches [3,48]. Of course, we have no means to provide a complete list of references here. In addition, for incompressible and nearly incompressible problems, it is also known that polytopal finite elements need to be constructed carefully to imitate volumetric locking free condition, i.e., ∇ · uh = 0.

(8)

This condition requires a special treatment in the polytopal finite element formulation. There are a few approaches found in the literature such as u/ p mixed formulation [7], virtual elements [49,50], and bubble functions [17]. Our approach forms a bubble-like scheme which extra degrees of freedom located at the centroid of elements are taken into account in the stiffness formulation and they can be then eliminated at the element level. To do this, we define a finite element space U h on each polytopal element defined as ⨁ Uh = Vh Vbh , (9) where the finite dimension spaces V h and Vbh are described as ⏐ ⏐ ⏐ { } ⏐ ⏐ ⏐ V h = uh ∈ H1 (Ω )d , uh ⏐ ∈ [Q(Ωe )]d , uh ⏐ ∈ C ∞ , u˜ h ⏐ ∈ C 0 , ∀Ωe ∈ T Ωe ⏐ Ωe ⏐ { } ∂Ωe h h 1 d h⏐ ∞ h⏐ Vb = ub ∈ H (Ω ) , ub ⏐ ∈ C , ub ⏐ = 0, ∀Ωe ∈ T , Ωe

∂Ωe

(10) (11)

in which Q(Ωe ) is basis functions, for example, Lagrange basis functions [47] or Wachspress basis functions [19], defined over Ωe . An explicit form of the displacements uh ∈ U h over Ωe ∈ T can be expressed as follows n nod ⏐ ∑ h⏐ u ⏐ = (Nei (x)Id )die + (Neb (x)Id )dbe , (12) Ωe

i=1

where n nod is the number of element’s vertices, Id is the unit matrix of dnd rank, die is the vector of unknown nodal values of the ith vertex of the element, Nie is the ith standard nodal basis shape function of element, dbe is the corresponding unknown of the centroid node, and Neb is the bubble shape function of element’s centroid node. Our approach is valid for various types of shape functions. To demonstrate this applicability, we have implemented and investigated two prevalent shape functions, namely Lagrange and Wachspress basis functions with the corresponding piecewise-linear and rational bubble shape functions. Generally, the present approach can be straightforwardly applied to other basis functions used in meshfree methods [51] or isogeometric analysis [52]. Now we illustrate Wachspress shape functions for a 2D domain as follows [53] A(vi−1 , vi , vi+1 ) wi (x) , wi (x) = = det (pi−1 , pi ), Nei (x) = ∑n j (x) w A(x, v i−1 , vi )A(x, vi , vi+1 ) j=1

(13)

where A(x1 , x2 , x3 ) is the signed area of the triangle constructed from three vertices [x1 , x2 , x3 ], pi is defined by pi = ni / h i (x) in which ni is the outward unit normal vector to the edge ei = [vi , vi+1 ] and h i (x) is the perpendicular distance from x to the edge ei as shown in Fig. 2. Additionally, its associated bubble shape function is defined by [54] Neb (x) =

1 ,  2 ∑n b j   1+ x−x j=1 w (x)

(14)

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

409

Fig. 2. Wachspress coordinates defined over subtriangles.

Fig. 3. An illustration of bubble shape functions at centroid node of some n-gon elements.

which is a rational function on Ωe , where xb is coordinates of internal (bubble) node. A 3D view of bubble functions over some typical n-gon elements are illustrated in Fig. 3. Derivatives of these functions are given by ⎛ ⎞ n nod ∑ ∇wi (x) ∇ Nei (x) = Nei (x) ⎝Ri − Nej (x)R j ⎠ , Ri = = pi−1 + pi , wi (x) j=1  2 ∑n nod j (15) ∑n nod j j w (x) + x − xb  2(x − xb ) j=1 j=1 w (x)R b ∇ Ne (x) = − . ( )2  2 ∑n nod j 1 + x − xb  j=1 w (x) Owing to the rational nature of Wachspress shape functions, numerical integration of the weak form integrals using standard Gaussian quadrature entails errors that prevent the convergence of solution when the mesh refinement reaches a certain level as pointed out in [30]. To tackle this burden, a technique called Gradient Correction is adopted. The core idea of this technique is to correct the gradients of these rational functions as follows (∮ ) ∫ 1 ∇ Nei (x) = ∇ Nei (x) + Nei (x)ndΓ − ∇ Nei (x)dΩ , |Ωe | ∂Ωe Ω (∮ ) ∫e (16) 1 b b b ∇ Ne (x) = ∇ Ne (x) + Ne (x)ndΓ − ∇ Neb (x)dΩ , |Ωe | ∂Ωe Ωe where |Ωe | is the area of the polygonal element and n is the unit normal vector along the element boundary. The consistency of this technique is illustrated through several numerical examples in Section 5.

410

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 4. The shape function definition of piece-wise linear basis.

For non-convex element cases, piecewise linear shape functions are built through children triangular sub-domains of the polygonal element as shown in Fig. 4. The boundaries of each T3 sub-domain are defined by two adjacent vertices of an edge and a virtual point at the centroid. The shape functions inherit some properties from polynomial functions such as Kronecker-delta, the partition of unity, non-negative. From this definition, the piecewise linear basis functions comply with the following rule to adopt the Kronecker delta property: ⎧ xi ≡ x j ⎨ 1 0 xi ̸= x j φei (x) = δi j ≡ (17) ⎩ 1/n nod xi ≡ x¯ Afterwards, the poly-shape functions and the first corresponding derivatives are evaluated over sub-domains as follows 3 ∑ j Nei (x) = NT3 (x)φei (x j ) j=1

∇ Nei (x) =

3 ∑

, x ∈ Ω T3

(18)

j ∇ NT3 (x)φei (x j )

j=1 j i where NT3 (x) and ∇ NT3 (x) are the basis functions and T3 i Ω with φe (x j ) obtained from Eq. (17), respectively.

first derivatives of linear triangular elements on sub-domains

In addition, to enhance the performance of the present approach, especially incompressible elasticity, we will later introduce a new assumed strain field and a weak form formulation derived from the Hu–Washizu variational principle. 3. Polytopal composite finite elements 3.1. An assumed strain field In this section we construct an assumed strain field which is obtained by a projector from the space of the ˜ compatible strains Sc into the[ space of the polynomial strains ] in S. Let us consider a 3D general form of the T polynomial strain field ϵ˜ = ϵ˜11 ϵ˜22 ϵ˜33 γ˜12 γ˜23 γ˜31 being described by the complete polynomial of an arbitrary degree k as follows ϵ˜ T (x) = β 0 + β 1 x1 + β 2 x2 + β 3 x3 + · · · ≡ S(x)β, where x = (x1 , x2 , x3 ) ⎡ β 00 β 01 ⎢β 10 β 11 β=⎢ ⎣β 20 β 21 β 30 β 31

(19)

is an arbitrary point in Ω , and β is a matrix of strain unknowns, which is defined as ⎤ β 02 β 03 ... β 12 β 13 ...⎥ ⎥, (20) β 22 β 23 ...⎦ β 32 β 33 ... h

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

411

Fig. 5. Illustration of 2D integration mappings associated with the proposed approach based on Wachspress and/or piecewise linear basis functions for two selected polygonal elements. The one highlighted in teal is a convex six-gon and another one highlighted in orange is a concave seven-gon. Depending on which type of polygon (i.e. convex or concave), the integration procedure could be designed to switch between Wachspress and piecewise linear bases. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and S(x) denotes a vector of the Pascal triangle polynomials, which is described by [ ] S(x) = 1 x1 x2 x3 ... .

(21)

Note that, in general case, higher-order polynomial strain modes can be utilized when using higher-order polytopal finite elements. In the current work, we restrict to using lower-order displacement approximations and hence a firstorder assumed strain field can be employed. In polytopal finite element formulation, the most popular approach is to subdivide the polytopal element into the number of sub-cells. The elemental stiffness matrix is contributed by whose performed over each sub-cell using the standard quadrature. This scheme is simple yet requires a large number of quadrature points that helps to achieve acceptable accuracy of solution. However, this approach does not pass patch test [48]. Several approaches in the literature have been devised to conduct the existing shortcomings of PFEMs. In ˜ verifying the current study, we propose a simple approach based on findings of an assumed strain field belonging S, ˜ Sc ⊂ S ⊂ S, which fits well every compatible strain field, ϵ ∈ Sc over sub-cells. To enhance active degrees of freedom for incompressible problems, we implement internal nodes (or bubble nodes) located at the centroid of the polytopal elements as depicted in Fig. 5. We then compute compatible strains on each sub-cell derived from the displacement field defined over barycentric coordinates [21]. Our solution is to make a bridge between Sc and S˜ via a polynomial projection Πk : L2 (Ω )3(d−1) ↦−→ S˜ = Pk (Ω )3(d−1) , where k = 0, 1, 2, . . . is the⋃degree of polynomial. nT Now, we find ϵ˜ = Πk ϵ ∈ S˜ with the strain unknowns (20) over sub-cell, Ωei , and Ωe = i=1 Ωei leading to the following error minimization Πk ϵ =

1 argmin E, 2    ∀˜ϵ ∈S˜

(22)

412

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

where the functional E can be expressed through contribution of error on sub-cells as follows ∫ nT ∫ ∑ ( )T ( ) ϵ˜ − ϵ i ϵ˜ − ϵ i dΩ , E = ∥˜ϵ − ϵ∥2L2 (Ωe ) = (˜ϵ − ϵ)T (˜ϵ − ϵ) dΩ = Ωe

(23)

Ωei

i=1

in which εi ∈ Sc is computed over sub-cell Ωei . Minimizing Eq. (22), i.e, ∂E = 0, ∀k, l = 0, 1, 2, 3, ∂βkl

(24)

we obtain the equations system with respect to unknowns β Mβ = Q

(25)

yielding β = M−1 Q

(26) k=1

where ⎡

M=

 1

⎢ nT n G ⎢x1 j ∑ ∑⎢ ⎢x2 j S(x)T S(x)dΩ = ⎢ i Ωe i=1 j=1 ⎢ x ⎣ 3j

nT ∫ ∑ i=1



x2 j

 x3 j

x12 j

x1 j x2 j

x1 j x3 j

x1 j x2 j

x22 j

x2 j x3 j

x1 j x3 j

x2 j x3 j

x32 j

...

...

...

x1 j

...

⎤ ... ...⎥ ⎥ ⎥ i ...⎥ ⎥ w j ∥J j ∥, ⎥ ...⎦

(27)

...

and Q=

nT ∫ ∑ i=1

Ωei

(

)T ϵ i (x)Si (x) dΩ

⎡

i ϵ11

⎢ i n T n G ⎢ x 1 j ϵ11 ∑ ∑⎢ ⎢x2 j ϵ i = 11 ⎢ i=1 j=1 ⎢ x ϵ i ⎣ 3 j 11

i ϵ22

i ϵ33

i γ12

i γ23

i γ31

i x1 j ϵ22

i x1 j ϵ33

i x1 j γ12

i x1 j γ23

i x1 j γ31

i x2 j ϵ22 i x3 j ϵ22

i x2 j ϵ33 i x3 j ϵ33

i x2 j γ12 i x3 j γ12

i x2 j γ23 i x3 j γ23

i x2 j γ31 i x3 j γ31

...

...

...

...

...

...

where n G , w j , Jij are a number sub-cell, Ωei , respectively.

⎤ ... ⎥ ...⎥ ⎥ i ...⎥ ⎥ w j ∥J j ∥, ⎥ ...⎦ ...

(28)

of quadrature points, quadrature weight and Jacobian determinant used in each

Now, by substituting Eq. (26) to Eq. (19), one writes a first order (k = 1) polynomial strain field as follows ϵ˜ T = Π1 ϵ = S(x)M−1 Q.

(29)

By substituting Eq. (12) to Eqs. (29) and (38), we express the displacement strain relation on each polytopal element as nT ∫ ∑ ( i )T ϵ˜ T = Π1 ϵ = S(x)M−1 Q = S(x)M−1 ϵ (x)Si (x) dΩ Ωei

i=1 ⎛ ⎞ n T n G1 ∑ ∑( ) ( ) T T = S(x)M−1 ⎝ Si (x j ) Bie (x j )de w j ∥Jij ∥⎠ i=1 j=1

⎛ ⎞ n T n G1 ( )( ∑ ∑ ( ) ) T T = dTe ⎝ S(x)M−1 Si (x j ) Bie (x j ) w j ∥Jij ∥⎠ i=1 j=1

=

dTe B˜ ke (x),

(30)

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

413

where Bie is the compatible strain–displacement matrix, de is the vector of nodal displacements of the polytopal ˜ e (x) evaluated at the lth Gaussian point (l = 1, 2, . . . , n G2 as specified in Eq. (43)) of the kth element, and B sub-cells is defined as n T n G1 ( ∑ ∑ ( )T ) i B˜ ke (xl ) = Sk (xl )M−1 Si (x j ) Be (x j )w j ∥Jij ∥, (31) i=1 j=1

where n G1 is the number of Gaussian points that depends on which type of basis is used. For instance, if the basis is piecewise linear, one Gaussian point per sub-cell is enough. However, if the basis is Wachspress, then at least three Gaussian points per sub-cell are required (as depicted on the left of Fig. 5). Note that ∫such an assumed strain field (29) satisfies the orthogonality condition from the Hu–Washizu principle [47], i.e, Ω δσ T (˜ϵ − ϵ)dΩ = 0, ∀δσ ∈ S ∩ P1 (Ω )3(d−1) , where P1 (Ω ) denotes the set of polynomials of degree 1. We will numerically prove that an approximate solution derived from Eq. (12) is convergent and passes a linear patch test. However, the aforementioned polytopal finite element formulations still do not work well for nearly incompressible elasticity in 3D.1 To overcome this shortcoming, we devise a pair (Π0 , Π1 ) of projection operators in accordance with a pair of (Π0 ϵ v , Π1 ϵ d ) ∈ S˜ × S˜ such that (Π0 ϵ v , Π1 ϵ d ) =

1 2

E,

argmin   

(32)

˜ S˜ ∀(˜ϵ v ,˜ϵ d )∈S×

where E=

nT ∫ ∑ i=1

( Ωei

ϵ˜ v −

)2 ϵ iv

dΩ +

nT ∫ ∑

(

Ωei

i=1

ϵ˜ d − ϵ id

)T (

) ϵ˜ d − ϵ id dΩ .

To end this, we explicitly write a constant projector Π0 , k = 0 such as nT ∫ nT ∫ 1 ∑ mT ∑ Π0 ϵ v = ∇ · ui h dΩ = εi dΩ , |Ωe | i=1 Ωei |Ωe | i=1 Ωei where |Ωe | denotes the “volume”2 of the element. Hence, the functional (33) reduces to nT ∫ ∑ ( )T ( ) E= ϵ˜ d − ϵ id ϵ˜ d − ϵ id dΩ , i=1

(33)

(34)

(35)

Ωei

where ϵ id is considered as a deviatoric strain field defined on Ωei . Now ϵ d can be derived from Eq. (19). Following same procedures as what given in Eqs. (24)–(26), one obtains ) ( 1 T T −1 T (36) Π1 ϵ d = Qd M S = I − mm Π1 ϵ, 3 where Qd =

n

nT ∫ ∑ i=1

(

Ωei

T ∑ )T ϵ id S dΩ =

i=1

∫ Ωei

(

( ) ( ) )T 1 1 ϵ i S dΩ I − mmT = Q I − mmT , 3 3

( ) 1 T with =m ϵ , = I − mm ϵ i , m = [1, 1, 1, 0, 0, 0]T and the identify matrix I. 3 Therefore, we obtain an assumed strain field for nearly incompressible elasticity as follows ( ) 1 1 1 T ϵ˜ = mΠ0 ϵ v + Π1 ϵ d = mΠ0 ϵ v + I − mm Π1 ϵ. 3 3 3 ϵ iv

1 2

T i

(37)

ϵ id

The present formulation still performs well for nearly incompressible 2D elasticity. We use this term as the area for 2D and the volume for 3D.

(38)

414

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

3.2. A stiffness matrix formulation From the Hu–Washizu principle, we obtain a weak form of the aforementioned strain field and a discrete solution uh ∈ U h that satisfies: a(uh , v) = ⟨ℓ, v⟩, ∀v ∈ U0h ,

(39)

where a(uh , v) = λ

∫

ϵ˜ v (uh )˜ϵ v (v)dΩ + 2µ Ω

∫

ϵ˜ T (uh )Dµ ϵ˜ (v)dΩ .

(40)

Ω

Substituting Eq. (12) into Eq. (39), we obtain systems of algebraic equations as follows ˜ = F, Kd

(41)

where d is the global nodal displacement vector, and K is the global stiffness matrix: ˜ = K

ne ⋀

˜ e, K

(42)

e=1

with the elemental stiffness, which is computed using three quadrature points on each sub-cell: n T n G2 ( ) ∑ ∑ ˜e = ˜ Te (x j )mmT B˜ e (x j ) + 2µB ˜ Te (x j )Dµ B˜ e (x j ) w j ∥Jij ∥, K λB

(43)

i=1 j=1

where n G2 is the number of Gaussian points used for computing stiffness matrix which is determined by the degree of the polynomial utilized to approximate the strain field, e.g. for first-order approximation (k = 1), three integration points are in demand (as depicted on the right of Fig. 5) which can be attributed to the product of two linear bases that produces a second-order basis presented in Eq. (31), and F is the global load vector denoted by: ∫ ∫ NT tdΓ (44) NT bdΩ + F= Ω

Γt

□ Remark: The discrete solution uh ∈ U h derived from Eq. (39) is convergent and stable. It is important to note that the approximation based on the assumed strain field helps to enhance significantly the accuracy of the numerical solution but it does not improve the convergence rate. This means that the present method produces convergence of order (h 2 ) in the displacement error norm and (h) in the energy error norm. 4. Numerical implementation In this section, we provide key algorithms for straightforwardly implementing the present method. For the sake of clarity of numerical implementation, we illustrate some main algorithms for 2D problems using polygonal elements, but these algorithms can be readily adapted to polyhedral elements. Algorithm 1: Compute matrix M 1 2 3 4 5 6 7 8 9 10 11 12 13

Input: P of size (n T + 1) × 2 containing nodal coordinates of polygonal element including center node Output: M Create M as an 3-by-3 zero-initialized matrix Define quadrature rule of the sub-cell in natural coordinates: ξ , w for i = 1 to n T do for j = 1 to n G do Evaluate 3-node shape functions and derivatives at ξ j : NT3 (ξ j ), ∇NT3 (ξ j ) Compute Jacobian and physical coordinates: J = ∇NT3 (ξ j )Pi , xij = [x1 j , x2 j ] = NT3 (ξ j )Pi , S = [1, x1 j , x2 j ] Evaluate matrix M: M = M + ST S det(J)w j end end Return M

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

415

Algorithm 2: Compute gradient correction matrix C 1 2 3 4 5 6 7 8 9 10 11 12 13

14 15 16 17 18 19

Input: P of size (n T + 1) × 2 containing nodal coordinates of polygonal element including center node Output: C Create Cs and Cv as (n T + 1) × 2 zero-initialized matrices Compute area of the polygonal element: Ae Define n G E quadrature rule for edges of sub-cell on natural coordinates: ξ , wξ Define n G quadrature rule of the sub-cell on natural coordinates: η, wη for i = 1 to n T do Extract coordinates of 2 nodes associated to ith edge: Ei Use Ei to compute normal vector of ith edge: n for js = 1 to n G E do Evaluate 2-node shape functions and derivatives at ξ js : NL2 (ξ js ), ∇NL2 (ξ js ) Compute Jacobian and physical coordinates: Js = ∇NL2 (ξ js )Ei , xijs = NL2 (ξ js )Ei Evaluate Wachspress and bubble basis functions at xijs : N(xijs )    1×(n T +1) √ Compute Cs = Cs + N(xijs )T nT Js JTs (wξ ) js end for jv = 1 to n G do Evaluate 3-node shape functions and derivatives at η jv : NT3 (η jv ), ∇NT3 (η jv ) Compute Jacobian and physical coordinates: Jv = ∇NT3 (η jv )Pi , xijv = NT3 (η jv )Pi Evaluate derivatives of Wachspress and bubble basis functions at xijv : ∇N(xijv )    2×(n T +1)

20 21 22 23

Compute Cv = Cv + ∇N(xijv )T det (Jv )(wη ) jv end end Return C = 1/Ae (Cs − Cv )

5. Numerical examples In this section, some benchmark problems are studied. For reference, some elements under investigation are labeled as follows MINI: Mixed displacement/pressure finite element enriched with cubic bubble functions [55]. PEn-PL: Polygonal elements with n-sides based on piecewise linear shape functions. PEn-WP: Polygonal elements with n-sides based on Wachspress shape functions [19]. PEn-GCWP: Polygonal elements with n-sides based on Wachspress shape functions with gradient correction [30]. • PEn-LS: Polygonal elements with n-sides based on Wachspress shape functions with linear smoothing technique [33]. • PCEn-PL: Polygonal composite elements (or the present element) with n-sides based on assumed strain fields and piecewise linear shape functions. • PCEn-GCWP: Polygonal composite element with n-sides based on assumed strain fields and Wachspress shape functions with gradient correction.

• • • •

5.1. Patch test:linear reproducibility/convergence This example numerically examines the convergence of the present method in the context of the standard patch test. We consider a patch as shown in Fig. 6. The exact displacements and stresses are given by u 1 = 10−3 (x1 + x2 /2), u 2 = 10−3 (x1 /2 + x2 ), σ11 = σ22 = 1333, τ12 = 400

(45) (46)

416

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Algorithm 3: Compute strain–displacement array B and H = M−1 ST for all sub-triangles of polygonal element. 1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

Input: P, n G1 , f l P of size (n T + 1) × 2 contains nodal coordinates of polygonal element including center node, n G1 is the number of quadrature points, f l is a flag to specify which type of basis functions to be used Output: B, H Call Algorithm 1 to compute M using P Create G as an 3 × 2(n T + 1) × n G1 × n T array Create H as an 3 × n G1 × n T array Define n G1 quadrature rule of the sub-cell on natural coordinates: ξ , w if fl is GCWP then Call Algorithm 2 to compute gradient correction matrix C end for i = 1 to n T do for j = 1 to n G1 do Evaluate 3-node shape functions and derivatives at ξ j : NT3 (ξ j ), ∇NT3 (ξ j ) Compute Jacobian, physical gradients and physical coordinates: J = ∇NT3 (ξ j )Pi , xij = [x1 j , x2 j ] = NT3 (ξ j )Pi → S = [1, x1 j , x2 j ] Compute M−1 ST and distribute its entries to H if fl is T3 then ∇N(xij ) = (J)−1 ∇NT3 (ξ j )    2×3

18 19

else if fl is WP or GCWP then Evaluate derivatives of Wachspress and bubble basis functions at xij : ∇N(xij )   

2×(n T +1)

20 21 22 23 24 25 26 27 28

if fl is GCWP then ∇N(xij ) = ∇N(xij ) + C end end Create Be as an 3 × 2(n T + 1) matrix and distribute ∇N(xij ) to corresponding positions in Be Compute and distribute Be det(J)w j to B end end Return B, H

Here we define the patch with one element entirely surrounded by neighboring elements. As resulted in Table 1, PEn-GCWP, PEn-LS, PCEn-GCWP and PCEn-PL pass the patch test with machine precision, while PEn-WP fails. Next we investigate a 3D patch test for the present method. A linear patch test is prescribed by the displacement field on the boundary as follows ⎛ ⎞ ⎛ ⎞ u1 0.1 + 0.1x + 0.2y + 0.2z ⎝u 2 ⎠ = ⎝0.05 + 0.15x + 0.1y + 0.2z ⎠ (47) u3 0.05 + 0.1x + 0.2y + 0.2z Several polyhedral elements are displayed in Fig. 7. Table 2 describes displacements at an arbitrary node (let us say 25th node, for instance) and error norms in displacement and energy. We observe that PEn-GCWP, PEn-LS and PCEn-GCWP nearly pass the patch test. PCEn-PL passes the patch test with machine precision, while PEn-WP fails. Based on the above examinations, in the subsequent Sections 5.4 and 5.5, we only take the PCEn-PL element into the discussion about the present approach.

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

417

Algorithm 4: Compute element stiffness matrix Ke 1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

Input: P, s f l, λ, µ, f l P of size (n T + 1) × 2 contains nodal coordinates of polygonal element including center node, s f l is a flag to indicate which stress state is assumed, f l is a flag to specify which type of basis functions to be used Output: Ke Create Kl as an 2(n T + 1) × 2(n T + 1) zero-initialized matrix Specify number of quadrature points to be used to evaluate strain–displacement matrix: n G2 Call Algorithm 3 to compute matrices B and H using P, n G2 and f l if sfl is PLANESTRESS then Define λ = 2λµ/(λ + 2µ) else if sfl is PLANESTRAIN then Define m = [1, 1, 0]T and compute area of the polygonal element Ae Create B˜ v as 1 × 2(n T + 1) zero-initialized matrix for i = 1 to n T do for j = 1 to n G1 do ˜ v + Bij Compute B˜ v = B end end ˜v B˜ v = (mT /Ae )B end Define 3 points quadrature rule of the sub-cell on natural coordinates: ξ , w for k = 1 to n T do for l = 1 to n G2 do Evaluate 3-node shape functions and derivatives at ξ l : NT3 (ξ l ), ∇NT3 (ξ l ) Compute J = ∇NT3 (ξ l )Pk , xlk = [x1l , x2l ] = NT3 (ξ l )Pk → S = [1, x1l , x2l ] Create B˜ lk as an 3 × 2(n T + 1) zero-initialized matrix for i = 1 to n T do for j = 1 to n G1 do ( ) Compute B˜ lk = B˜ lk + SHij Bij end end if sfl is PLANESTRAIN then ( ) Compute B˜ lk = 1/3mB˜ v + I − 1/3(mmT ) B˜ lk end ( ) Compute Kloc = Kloc + λ(B˜ lk )T mmT B˜ lk + 2µ(B˜ lk )T Dµ B˜ det(J)wl end end [

35

Define Kloc

Kll = Kql

] Klq −1 and return Ke = Kll − Klq Kqq Kql Kqq

5.2. Study of convergence and accuracy

5.2.1. Cantilever beam Consider a cantilever beam subjected to a parabolic load P = 1000 at the free end as shown in Fig. 8. The geometry dimensions are taken as: length L = 48, height D = 12 and thickness t = 1. The isotropic material is utilized with Young’s modulus E = 3 × 107 .

418

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 6. Several arbitrary polygonal elements for the patch test (E = 1 × 106 , ν = 0.25, thickness t = 0.001). Table 1 A 2D patch test. PEn-WP PEn-GCWP PEn-LS PCEn-PL PCEn-GCWP a Exact

a u 10 1

a u 10 2

Displacement error norm

Energy error norm

Patch test

1.6246 × 10−4

1.4295 × 10−4

2.8095 × 10−8

1.65 × 10−4 1.65 × 10−4 1.65 × 10−4 1.65 × 10−4

1.5 × 10−4 1.5 × 10−4 1.5 × 10−4 1.5 × 10−4

8.1043 × 10−20 8.8442 × 10−19 4.3582 × 10−21 9.4094 × 10−20

0.0032 2.1909 × 10−6 2.1909 × 10−6 2.1909 × 10−6 2.1909 × 10−6

Fail Pass Pass Pass Pass

displacements of 10th node: u 1 = 1.65 × 10−4 , u 2 = 1.5 × 10−4 .

Table 2 A 3D patch test. PEn-WP PEn-GCWP PEn-LS PCEn-PL PCEn-GCWP a Exact

a u 25 1

a u 25 2

a u 25 3

Displacement error norm

Energy error norm

Patch test

0.3888 0.3885 0.3885 0.3885 0.3885

0.3156 0.3125 0.3125 0.3125 0.3125

0.3379 0.3385 0.3385 0.3385 0.3385

0.0010 5.9926 × 10−12 6.7006 × 10−12 6.0799 × 10−16 1.4951 × 10−11

0.0126 6.005 × 10−11 1.0159 × 10−10 4.4761 × 10−15 9.5187 × 10−11

Fail Nearly pass Nearly pass Pass Nearly pass

displacements of 25th node: u 1 = 0.3885, u 2 = 0.3125, u 3 = 0.3385.

The analytical solution of this benchmark is available in [56]. The displacements in the x1 and x2 directions are given as [ ( )] D2 P x2 2 u1 = (6L − 3x1 )x1 + (2 + ν˜ ) x2 − 4 6 E˜ I (48) [ ] 2 P D x 1 2 2 u2 = − 3˜ν x2 (L − x1 ) + (4 + 5˜ν ) + (3L − x1 )x1 , 4 6 E˜ I { { E, ν for plane stress ˜ where E = ν˜ = . 2 E/(1 − ν ), ν/(1 − ν) for plane strain and the stress components are ( ) P(L − x1 )x2 P D2 2 σ11 (x1 , x2 ) = (49) , σ22 (x1 , x2 ) = 0, τ12 (x1 , x2 ) = − − x2 I 2I 4 with I = D 3 /12.

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

419

Fig. 7. Several polyhedral elements for the 3D patch test (E = 1, ν = 0.25). The data to construct this mesh is taken from the Linear Patch test example of the paper [33]. All elements are convex and polygon faces are planar.

Fig. 8. Cantilever beam’s geometry and boundary conditions.

Fig. 9. An illustration of polygonal meshes for the cantilever beam.

Fig. 9 illustrates coarse and slightly fine meshes using polygonal elements. In Fig. 10, we plot the numerical solutions of the vertical displacements along the neutral axis for both compressible and the nearly incompressible materials using the first mesh. For compressible case (under plane stress conditions), all methods produce reliable solutions in comparison with the exact value. For incompressible case (under plane strain conditions), the present

420

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 10. Vertical displacement of nodes associated to the neutral axis of the cantilever beam under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

Fig. 11. A series of polygonal meshes used for error norm estimations of the cantilever beam.

elements, PCEn-PL and PEn-GCWP, overcome volumetric locking while PEn-LS and PEn-GCWP perform poorly. It is interesting that the PCEn-PL and PEn-GCWP elements yield very accurate results. To reveal the good performance of the present method, we investigate convergences of strain energy as well as displacement and energy error norms versus a series of refined meshes displayed in Fig. 11. The “characteristic length” parameter of the meshes, h, in the following figures is defined as a dimensionless length given by √ h = 1/ n do f , where n do f is the total number of degree of freedoms of the entire domain. Convergences of strain energy versus the number of degree of freedoms n do f for compressible and the nearly incompressible materials are plotted in Fig. 12 with the corresponding exact values computed analytically as 4.47467 (Nm) and 3.44 (Nm), respectively. For compressible material, various elements are tested as depicted in Fig. 12a. It is found that all elements produce a lower bound of the exact strain energy when the mesh is refined. From the results we see that two group of elements gives nearly the same convergence rates and accuracy. The first

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

421

Fig. 12. Strain energy curves of the cantilever beam problem evaluated by various approaches under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

group which consists of PEn-LS and PEn-GCWP elements gives the worst convergence and accuracy. In contrast, the second group which comprises PCEn-PL and PCEn-GCWP elements provides very highly accurate strain energy as expected. More notably, the standard four-nodal quadrilateral element (FEM-Q4) using full integration produces an intermediate level of accuracy and convergence where its convergence curve lies between those of the first and second group. For nearly incompressible material, numerical tests are conducted based on PCEn-PL, PCEn-GCWP, and MINI elements. From the convergence curves shown in Fig. 12b, we see that the PCEn-PL and PCEn-GCWP elements again gain almost identical strain energy values and achieve a higher accuracy compared to that using MINI element. Figs. 13a and 14a plot the convergence rates for the compressible case. We observe that results of PEn-LS and PEn-GCWP are almost identical. The same observation is also found for PEn-LS and PEn-GCWP. As expected, solutions of the present elements are more accurate than those of PCEn-PL and PEn-GCWP. The accuracy of the present elements is about 28.5 times in the displacement error norm and about 3 times in the energy error norm higher than that of PEn-LS and PEn-GCWP. More importantly, we compare the present elements with the MINI element [55]. As shown in Figs. 13b and 14b, the present elements demonstrate excellent performance in incompressibility limit. Their accuracy is approximately 84 times higher in the displacement error norm and 10 times higher in the energy error norm than that of MINI. In terms of the convergence rate, the convergence rate of the present method is asymptotically the theoretical value, i.e., 2.0 in the displacement norm and 1.0 in the energy norm. 5.2.2. A cylindrical pipe subjected to an inner pressure In this subsection, we consider the thick-wall cylindrical pipe problem having internal and external radii of a = 1 and b = 2, respectively. The internal face of the pipe is subjected to a constant pressure of p = 8 as depicted in Fig. 15. By taking advantage of the axisymmetric property of the problem, only the upper right quadrant of the pipe is analyzed. Fig. 16 illustrates the discretizations of the domain using n-node polygonal elements. Plane strain condition is assumed and nearly incompressible elastic material is adopted with Poisson ratio ν = 0.4999999 and Young’s modulus E = 21000. The left and bottom edges are constrained by symmetric conditions meanwhile the outer boundary is assigned traction free. The analytical solution of radial and tangential displacement [56] is ( ) (1 + ν)a 2 p b2 u r (r ) = (1 − 2ν)r + ; u ϕ = 0, (50) E(b2 − a 2 ) r

422

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 13. Displacement error norms of the cantilever beam problem under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

Fig. 14. Energy error norms of the cantilever beam problem under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

while the stress components are given by ( ) ( ) a2 p b2 a2 p b2 σr (r ) = 2 1 − 2 ; σϕ (r ) = 2 1 + 2 ; σr ϕ = 0, b − a2 r b − a2 r

(51)

where (r, ϕ) are the polar coordinates from the origin O, and ϕ is determined counter-clockwise from the positive x-axis (Fig. 15). For the convergence rate analysis, a series of sequentially refined meshes with 128, 288, 648, 1225 and 2048 elements are utilized as depicted in Fig. 16. The “characteristic length” parameter of the meshes, h, is taken as described in the first test on cantilever beam. MINI element and the present approach are examined in this problem. Convergence in displacement norm, energy norm and pressure norm are shown in Figs. 17a, 17b, 17c, respectively. It is seen that MINI element, the proposed PCEn-PL, and PCEn-GCWP elements exhibit optimal convergence (with rates of 2 and 1 in the displacement norm

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

423

Fig. 15. A cylindrical pipe subjected to an inner pressure and its symmetric quarter model with boundary conditions.

Fig. 16. A series of polygonal meshes used for error norm estimations of the cylindrical pipe.

and energy norm, respectively). Nearly the same levels of error are accomplished for the PCEn-PL and PCEn-GCWP elements while MINI element again gives a much larger error. The averaged accuracy in displacement norm, energy

424

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 17. Error norms of the cylindrical pipe problem under nearly incompressible condition (ν = 0.4999999).

norm and pressure norm of the present elements are, in turn, 8.8, 5.8 and 10 times higher than that of MINI element. These findings again confirm that the proposed approach gives reliable and highly accurate solutions. As shown in Fig. 17c, the pressure norm error given by the PCEn-PL element is slightly lower than that of the PCEn-GCWP element, but overall the convergence rates of both elements are the same and equal to 1.5. It is noteworthy that these two elements outperform MINI element both in accuracy and convergence rate. This test problem again demonstrates that the proposed approach is stable and highly accurate for nearly incompressible materials. 5.2.3. Cook membrane in incompressibility limit: stability test In this example, we numerically examine the inf–sup (or Ladyzhenskaya–Babuˇska–Brezzi) condition [57] for the present method. This aims to prove that the present method is stable in incompressibility limit. Then, we also demonstrate excellence performance of the present elements for a structure in deformation dominated by both a bending deformation and nearly incompressible materials. To end this, we consider the well-known Cook’s membrane problem [58] under plane strain conditions as illustrated in Fig. 18a. The membrane is a clamped tapered panel subjected to an in-plane shearing load. This leads to deformation dominated by a bending deformation. For computation, we use Young’s modulus E = 1 and Poisson ratio ν = 0.4999999. The exact solution of the problem is unknown. The best reference value of the vertical displacement at center tip section (indicated by the point A in Fig. 18a) and strain energy are found to be 18.5002 and 9.27769, respectively [59].

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

425

Fig. 18. Cook’s membrane problem and the inf–sup test.

First of all, the numerical inf–sup values are evaluating using 5 gradually finer meshes depicted in Fig. 19. √ The ultimate goal is to inspect the tendency of the inf–sup values, βh , when the characteristic length h = 1/ n do f decreases. If β asymptotically approaches a value greater than zero when h decreases, the satisfaction of the inf–sup condition is guaranteed. Otherwise, the approach does not pass the LBB condition. As shown in Fig. 18b, the inf– sup values for the PCEn-PL and PCEn-GCWP elements remain bounded by a positive constant with refinement of meshes. The same observation is found for MINI element as expected. On the other hand, inf–sup values calculated for the PEn-GCWP element approach to zero when the mesh is finer which indicates that this element is failed to satisfy the LBB condition. This numerical study confirms that the proposed approach is LBB compliant irrespective of which kind of basis functions are utilized. Second, the convergence behavior of the center tip section displacement and strain energy versus the number of degree of freedoms are investigated. The numerical study is performed by making use of PCEn-PL, PCEn-GCWP elements from the present approach as well as MINI element. The obtained convergence curves of displacement and strain energy solution are depicted by Figs. 20a and 20b, respectively. From Fig. 20b, it is shown clearly that all three elements produce lower bound solutions in terms of strain energy. In general, with the increase of DOFs, these two quantities obtained from three elements reveal good convergent trends toward the reference values. However, it is observed that PCEn-PL and PCEn-GCWP elements give much more accurate displacement and strain energy solutions that far better than those of MINI element at the same number of degree of freedoms. More interestingly, it can seen that PCEn-PL and PCEn-GCWP elements can produce very close solutions to the references even when the coarsest mesh is used. Finally, the pressure distribution inside the membrane is judged by the present approach utilizing PCEn-PL element with the fifth mesh as illustrated in Fig. 19. The produced result is demonstrated in Fig. 21b which renders smooth pressure over entire computational domain without any spurious oscillation and compatible with that of MINI element in Fig. 21a. This observation indicates that the proposed approach is stable under incompressibility limit. 5.3. 2D bracket: Capacity of method to complicated geometry The following problem is investigated with the aim of demonstrating the capacity of the proposed approach in dealing with a geometry having intricate cutouts. With this aim, we consider a 2D bracket subjected to a distributed load on the top edge with clamped constraint conditions on the right edge as described in Fig. 22a. The geometric dimensions of the problem domain are L = 0.7 m, T = 0.1 m and the distribute force f = 103 N/m. Linear elastic

426

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 19. A series of polygonal meshes used for convergent evaluations of the Cook’s membrane problem.

material with plane stress condition is assumed with Poisson ratio ν = 0.3, Young’s modulus E = 2.1 × 1011 . The n-gon polygonal mesh of the domain employed for computation is shown in Fig. 22b. To illustrate the accuracy and to estimate the errors of the results obtained from the present method, the results are compared with solutions that are computed with commercial finite element software ANSYS [60]. To quantify the error of the present method, von-Mises’ stress value is chosen as comparison quantity. To this end, a computational mesh with 2033 quadratic elements of type Q8 is employed for getting reference results from ANSYS. For the present method, the problem domain is decomposed into 1000 polygonal elements and analyzed using PCEn-PL and PCEn-GCWP elements. The comparison is accomplished by extracting von-Mises’ stress values along a curved-edge of the bracket as depicted in Fig. 23a. As shown in Fig. 23b, the von-Mises’ stress values obtained from PCEn-PL and PCEn-GCWP elements are nearly identical and capture really well with those from ANSYS. Fig. 24 illustratively depicts the von-Mises’ stress distribution in the bracket induced by the distributed loading by both ANSYS and the present method with the PCEn-GCWP element. It can be seen that the present method produces von-Mises’ stress distribution in compliant with that of ANSYS although a fewer number of degree of freedoms are used (2000 nodes of PCEn-GCWP element vs 6455 nodes of Q8 element). These studies imply that the present approach allows yielding highly accurate results with a significantly fewer number of DOFs in comparison with traditional elements.

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

427

Fig. 20. Displacement and strain energy convergent of Cook’s membrane problem (ν = 0.4999999) obtained from various methods.

Fig. 21. Nodal pressure of Cook’s membrane problem in nearly incompressible state (ν = 0.4999999).

5.4. Convergence study for 3D-dimensional elasticity We consider a 3D cantilever beam subjected to shear load at the free end. The beam’s geometry is shown in Fig. 25 where the dimensions are chosen to be L = 5, 2a = 1, and 2b = 1. The beam is constrained weakly by the analytical solution on the face x3 = L and subjected to a shear force F on the face x3 = 0. Isotropic and

428

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 22. A bracket subjected to a distributed load on the top edge with clamped constraint conditions on the right edge and its corresponding polygonal mesh used for computation. Source: http://hpfem.org.

Fig. 23. Comparison of von Mises stresses evaluated along the curved-edge of the bracket. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

homogeneous material is assumed with Young’s modulus E = 1.0, Poisson’s ratio ν = 0.3 for compressible and ν = 0.4999999 for nearly-incompressible conditions, respectively. The analytical solution of the Cauchy stress field

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

429

Fig. 24. Contour plot of von Mises stresses obtained from ANSYS and present method.

Fig. 25. The prismatic beam of length L with rectangular cross-section measured 2a in width and 2b in height.

and the associated displacement field is given in [4] as follows σ11 = σ21 = σ22 = 0 F σ33 = x2 x3 I ∞ F 2a 2 ν ∑ (−1)n sinh(nπ x2 /a) σ31 = sin(nπ x1 /a) 2 2 I π 1 + ν n=1 n cosh(nπ b/a) [ ] ∞ 3x12 − a 2 F b2 − x22 F ν 2a 2 ∑ (−1)n cosh(nπ x2 /a) σ32 = + − 2 cos(nπ x1 /a) , I 2 I 1+ν 6 π n=1 n 2 cosh(nπ b/a)

(52)

430

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 26. Illustration of four structured hexahedral meshes used in the convergent study. The number of elements in each mesh is annotated under each sub-figure.

Fν x1 x2 x3 E[ I ] ) 1 3 F ν( 2 2 u2 = x − x2 x3 − x3 EI 2 1 6 [ ( ) ( ) F 1 1 1 1 u3 = x2 νx12 + x32 + νx23 + (1 + ν) b2 x2 − x23 − a 2 νx2 EI 2 6 3 3 ] ∞ 4a 3 ν ∑ (−1)n sinh(nπ x2 /a) − 3 cos(nπ x1 /a) , π n=1 n 3 cosh(nπ b/a) u1 = −

(53)

where E is Young’s modulus, ν is Possion’s ratio, and I = 4ab3 /3 is the second moment of area around x1 -axis. For convergence study, two types of meshes are generated for numerical simulation. The first one is the structured hexahedral mesh with four levels of progressively refined meshes (2 × 2 × 10, 4 × 4 × 20, 8 × 8 × 40,

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

431

Fig. 27. Illustration of four unstructured polyhedral meshes used in the convergent study. The number of elements in each mesh is annotated under each sub-figure. These meshes are generated by using the open-source software namely cfMesh [61] which generally produces meshes with non-planar polygonal faces. However, by using Piecewise Linear basis functions for these non-planar face elements, our formulation still provides accurate results. Note that, although the four surface meshes of the beam look “structured”, the corresponding volume meshes are actually unstructured as they can clearly be seen through the transparent skin of the first mesh illustrated in the sub-figure (a).

16 × 16 × 80) and the second one is unstructured polyhedral mesh with four different levels including 100, 289, 664, and 1720 elements. The open-source meshing software cfMesh [61] is used to generate these polygonal meshes. Since this software mostly produces non-planar polygonal faces, only piecewise linear basis functions can be used for investigation. The number of quadrature points follows the rule established in Section 3, where for this specific case, one quadrature point is employed for computing Be matrix and four quadrature point are required for forming B˜ e matrix. Figs. 26 and 27 in turn illustrate these two types of meshes used in computation. For comparison purpose, PEn-PL and FEM-H8 elements are used as references while for the present approach we only consider PCEn-PL element, since only this element is quite simple in implementation and can pass the patch test to machine precision as pointed out in Section 5.1. More importantly, PCEn-PL can capture well with concave polyhedral elements, which becomes too difficult or even impossible for other elements using Wachspress basis functions. Figs. 28 and 29 present the errors in displacement and energy norms for both material conditions when meshes are refined. It is observed that the present PCEn-PL element can attain the same accuracy and convergence rates as canonical FEM-H8 element in case of compressible material. Meanwhile, for nearly-incompressible material, PCEn-PL shows optimal convergence rates for both error norms. This indicates that the present element is already

432

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 28. Displacement error norms of the prismatic beam problem under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

Fig. 29. Energy error norms of the prismatic beam problem under compressible (ν = 0.3) and nearly incompressible (ν = 0.4999999) conditions.

immune to volumetric locking while PEn-PL is subjected to volumetric locking, resulting in severe accuracy loss in energy norm and non-convergence in displacement norm. 5.5. Static analysis of hydraulic press machine: from model design to manufacture This example aims to prove the high capability of the proposed method by applying it to static analysis of a complicated 3D solid structure, namely Hydraulic press machine. One of the most popular applications of this machine is found in the automobile industry where many parts of a car are produced using forming and stamping operations. The hydraulic press machine with given geometry and boundary conditions is displayed in Fig. 30. Its basic dimensions are 3.5 m, 2.2 m and 3.8 m corresponding to length, width and height, respectively. The machine is made from steel with Poisson’s ratio ν = 0.3 and Young’s modulus E = 2.1 × 1011 . The fixed supports are placed at the feet (A). An assumed vertical load F = 106 N is applied at the red region (B) on the working bench. As a result, a corresponding reaction force is generated at the piston rod (C).

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

433

Fig. 30. CAD and mesh models of the hydraulic press machine: (a) its geometry with boundary conditions and (b) the corresponding polyhedral mesh.

Again, to evaluate the accuracy and efficiency of the present approach, we use the results obtained from ANSYS software as references. For the computation using ANSYS, a tetrahedral (T4) element mesh with 2,444,637 DOFs is employed while for the present approach, a mesh of unstructured polyhedral elements with 587,868 DOFs is used as given in Fig. 30b. The analysis based on the present approach is accomplished via PCEn-PL element. Numerical solutions of displacement and von-Mises fields produced from ANSYS and the PCEn-PL element are correspondingly presented in Figs. 31 and 32. From these two figures, it is clear that the PCEn-PL element is capable of attaining almost equivalent displacement and stress distribution over the whole structure while using a significantly fewer number of degree of freedoms as compared to the canonical T4-element (587,868 vs. 2,444,637 DOFs). This observation reconfirms the versatility and high accuracy of the proposed approach. The diagram in Fig. 33 illustrates the process of manufacturing this machine. Overall, there are five stages in total. Firstly, the frame of the machine is designed in a CAD software following input requirements about dimension. Next, the geometry is imported to the mesher to generate the mesh model. In the following step, this mesh is employed to analyze the behavior of the machine structure under assumed working conditions as well as to check the durability of the structure. After that, the design can be modified to satisfy the durable standards. Then, a small-scale 3D prototype of the design is fabricated thanks to 3D printer to plan the order for construction. Finally, the machine is welded from different parts which are cut from steel plates. 6. Conclusions We presented a new computational scheme for polytopal finite elements. By applying a polynomial projection, an assumed strain field was obtained from the sense of the least-squares approximation. The underlying stiffness formulation can be accurately evaluated without requirement of using large number of quadrature points. For nearly incompressible problems, we applied the constant projection to volumetric strain components while retaining the first-order projection to deviatoric strains. Through numerical tests, substantial remarks can be drawn as follows • The proposed approach passed the linear patch test in both 2D and 3D cases. • The present method verified the inf–sup condition, which proves the stability of the pressure field in incompatibility limit. • The method was less sensitive to concave and irregular elements, which demonstrates its good performance for arbitrary polytopal meshes.

434

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

Fig. 31. Profile of total displacements obtained from ANSYS and the present method.

Fig. 32. Distribution of von Mises stresses obtained from ANSYS and present method.

• The overhead time of computation included the projections during process of establishment of elemental stiffness formulation. • The present approach yielded highly accurate solutions and achieved optimal rate of convergence in both displacement and energy error norms. • An application of the proposed method to industrial problems has been confirmed with high reliability in comparison with commercial software.

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

435

Fig. 33. Cycle from model design to manufacture for hydraulic press machine.

The extension of the present method to nonlinear problems can be straightforwardly formulated. Furthermore, it can be extended to other numerical methods like meshfree or isogeometric analysis in which rational basis functions are often used. Finally, the extension of the proposed approach to quadratic serendipity polytopal elements will most likely be a topic of future study. Acknowledgments This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.01-2018.32. References [1] S. Ghosh, S. Moorthy, Elastic-plastic analysis of arbitrary heterogeneous materials with the voronoi cell finite element method, Comput. Methods Appl. Mech. Eng. 121 (1995) 373–409. [2] J. Zhang, N. Katsube, A polygonal element approach to random heterogeneous media with rigid ellipses or elliptical voids, Comput. Methods Appl. Mech. Eng. 148 (1997) 225–234. [3] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, Internat. J. Numer. Methods Engrg. 61(12) (2004) 2045–2066.

436

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

[4] J.E. Bishop, A displacement-based finite element formulation for general polyhedra using harmonic shape functions, Internat. J. Numer. Methods Engrg. 97 (2014) 1–31. [5] C. Talischi, G.H. Paulino, A. Pereira, I.M. Menezes, Polygonal finite element for topology optimization: a unifying paradigm, Internat. J. Numer. Methods Engrg. 82(6) (2007) 671–698. [6] C. Talischi, G.H. Paulino, A. Pereira, I.M. Menezes, Polytop: a matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes, J. Struct. Multidiscipl. Optim. 45(3) (2012) 329–357. [7] C. Talischi, A. Pereira, G.H. Paulino, I.M. Menezes, M.S. Carvalho, Polygonal finite elements for incompressible fluid flow, Internat. J. Numer. Methods Fluids 74 (2014) 134–151. [8] A. Pereira, C. Talischi, G.H. Paulino, I.F. M. Menezes, M.S. Carvalho, Fluid flow topology optimization in polytop: stability and computational implementation, Struct. Multidiscip. Optim. 54 (5) (2016) 1345–1364. [9] A. Tabarraei, N. Sukumar, Extended finite element method on polygonal and quadtree meshes, Comput. Methods Appl. Mech. Engrg. 197(5) (2008) 425–438. [10] E.T. Ooi, C.M. Song, F. Tin-Loi, Z.J. Yang, Polygon scaled boundary finite elements for crack propagation modeling, Internat. J. Numer. Methods Engrg. 91(3) (2012) 319–342. [11] A.R. Khoei, R. Yasbolaghi, S.O.R. Biabanaki, A polygonal finite element method for modeling crack propagation with minimum remeshing, Int. J. Fract. (2015) 123—148. [12] H. Chi, C. Talischi, O. LopezPamies, G.H. Paulino, Polygonal finite elements for finite elasticity, Internat. J. Numer. Methods Engrg. 101 (2015) 305–328. [13] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, Comput. Methods Appl. Mech. Engrg. 306 (2016) 216–251. [14] S. Nguyen-Hoang, D. Sohn, H.G. Kim, A new polyhedral element for the analysis of hexahedral-dominant finite element models and its application to nonlinear solid mechanics problems, Comput. Methods Appl. Mech. Eng. 324 (2017) 248—277. [15] A. Rajagopal, M. Kraus, P. Steinmann, Hyperelastic analysis based on a polygonal finite element method, Mech. Adv. Mater. Struct. 25 (2018) 930–942. [16] H. Nguyen-Xuan, A polygonal finite element method for plate analysis, Comput. Struct. 188 (2017) 45–62. [17] H. Nguyen-Xuan, S. Nguyen-Hoang, T. Rabczuk, K. Hackl, A polytree-based adaptive approach to limit analysis of cracked structures, Comput. Methods Appl. Mech. Engrg. 313 (2017) 1006–1039. [18] L. Perumal, A brief review on polygonal/polyhedral finite element methods, Math. Probl. Eng. (2018) 5792372, 22 pages. [19] E.L. Wachspress, A rational finite element basis, Academic Press, London, 1975. [20] M. Floater, Mean value coordinates, Comput. Aided Geom. Design 20 (2003) 19–27. [21] M.S. Floater, K. Hormann, G. Kos, A general construction of barycentric coordinates over convex polygons, Adv. Comput. Math. 24 (2006) 311–331. [22] N. Sukumar, E.A. Malsch, Recent advances in the construction of polygonal finite element interpolants, Arch. Comput. Methods Eng. 13 (1) (2006) 129. [23] S. Natarajan, S. Bordas, D.R. Mahapathra, Numerical integration over arbitrary polygonal domains based on Schwarz–Christoffel conformal mapping, Internat. J. Numer. Methods Engrg. 80(1) (2009) 103–134. [24] S.E. Mousavi, H. Xiao, N. Sukumar, Generalized Gaussian quadrature rules on arbitrary polygons, Internat. J. Numer. Methods Engrg. 82 (1) (2010) 99–113. [25] K.Y. Dai, G.R. Liu, T.T. Nguyen, An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics, Finite Elem. Anal. Des. 43(11-12) (2007) 847–860. [26] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons, Comput. Mech. 47 (2011) 535–554. [27] E.B. Chin, J.B. Lasserre, N. Sukumar, Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra, Comput. Mech. 56 (2015) 967–981. [28] D. Wirasaet, E.J. Kubatko, S.T. C. E. Michoski, J.J. Westerink, C. Dawson, Discontinuous Galerkin methods with nodal and hybrid modal/nodal triangular, quadrilateral, and polygonal elements for nonlinear shallow water flow, Comput. Methods Appl. Mech. Eng. 270 (2014) 113—149. [29] A. Cangiani, E.H. Georgoulis, P. Houston, Hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes, Comput. Methods Appl. Mech. Eng. 24 (2014) 2009—2041. [30] C. Talischi, A. Pereira, I.F.M. Menezes, G.H. Paulino, Gradient correction for polygonal and polyhedral finite elements, Internat. J. Numer. Methods Engrg. 102 (2015) 728–747. [31] H. Chi, C. Talischi, O. Lopez-Pamies, G.H. Paulino, A paradigm for higher-order polygonal elements in finite elasticity using a gradient correction scheme, Comput. Methods Appl. Mech. Engrg. 306 (2016) 216–251. [32] S. Natarajana, E.T. Ooi, I. Chiong, C. Song, Convergence and accuracy of displacement based finite element formulations over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation, Finite Elem. Anal. Des. 85 (2014) 101–122. [33] A. Francis, A. Ortiz-Bernardin, S. Bordas, S. Natarajan, Linear smoothed polygonal and polyhedral finite elements, Internat. J. Numer. Methods Engrg. 109 (2017) 1263–1288. [34] L.B.D. Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 199–214. [35] L.B.D. Veiga, F. Brezzi, L.D.M. and, Virtual elements for linear elasticity problems, SIAM J. Numer. Anal. 5 (2013) 794–812. [36] G. Manzini, A. Russo, N. Sukumar, New Perspectives on Polygonal and Polyhedral Finite Element Methods, 2014. [37] A. Astaneh, F. Fuentes, J. Mora, L. Demkowicz, High-order polygonal discontinuous petrovgalerkin (polydpg) methods using ultraweak formulations, Comput. Methods Appl. Mech. Eng. 332 (2018) 686–711.

H. Nguyen-Xuan, K.N. Chau and K.N. Chau / Computer Methods in Applied Mechanics and Engineering 355 (2019) 405–437

437

[38] P. Antonietti, P. Houston, G. Pennesi, Fast numerical integration on polytopic meshes with applications to discontinuous Galerkin finite element methods, J. Sci. Comput. 77 (2018) 1339—1370. [39] R.W. Clough, J.L. Tocher, Finite element stiffness matrices for analysis of plate bending, in: Proc. Conf. on Matrix Methods in Structural Mechanics, WPAFB, Ohio, 1965, pp. 515–545. [40] B.M.F. de Veubeke, A conforming finite element for plate bending, Int. J. Solids Struct. 4 (1968) 95–108. [41] B.M.F. de Veubeke, Displacement and equilibrium models in the finite element method, in: O. Zienkiewicz, G. Holister (Eds.), Stress Analysis, Wiley, New York, 2001, pp. 145–197, Chapter 9, Reprinted in International Journal for Numerical Methods in Engineering 52 (2001) 287–342. [42] W. Hackbusch, S.A. Sauter, Composite finite elements for the approximation of PDEs on domains with complicated micro-structures, Numer. Math. 75 (1997) 447–472. [43] Y. Guo, M. Ortiz, T. Belytschko, E.A. Repetto, Triangular composite finite elements, Internat. J. Numer. Methods Engrg. 47 (2000) 287–316. [44] P. Thoutireddy, J.F. Molinari, E.A. Repetto, M. Ortiz, Tetrahedral composite finite elements, Internat. J. Numer. Methods Engrg. 53 (2002) 1337—1351. [45] J.T. Ostien, J.W. Foulk, A. Mota, M.G. Veilleux, A 10-node composite tetrahedral finite element for solid mechanics, Internat. J. Numer. Methods Engrg. 107 (2016) 1145–1170. [46] L. Leonetti, M. Aristodemo, A composite mixed finite element model for plane structural problems, Finite Elem. Anal. Des. 94 (2015) 33–46. [47] O.C. Zienkiewicz, R.L. Taylor, The finte element method. Vol.1, Fifth ed., McGraw-Hill, 2000. [48] C. Talischi, G.H. Paulino, Addressing integration error for polygonal finite elements through polynomial projections: A patch test connection, Math. Models Methods Appl. Sci. 24 (2014) 1701–1727. [49] P. Wriggers, B.D. Reddy, W. Rust, B. Hudobivnik, Efficient virtual element formulations for compressible and incompressible finite deformations, Comput. Mech. 60 (2017) 253–268. [50] E. Artioli, S. de Miranda, C. Lovadina, L. Patruno, A stress/displacement virtual element method for plane elasticity problems, Comput. Methods Appl. Mech. Engrg. 325 (2017) 155–174. [51] G.R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method, second ed., CRC Press, 2009. [52] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [53] M.S. Floater, A. Gillette, N. Sukumar, Gradient bounds for wachspress coordinates on polytopes, SIAM J. Numer. Anal. 52 (1) (2014) 515–532. [54] E.A. Malsch, G. Dasgupta, Shape functions for polygonal domains with interior nodes, Internat. J. Numer. Methods Engrg. 61(8) (2004) 1153–1172. [55] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344. [56] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, 1970. [57] D. Chapelle, K.J. Bathe, The inf-sup test, Comput. Struct. 47 (1993) 537–545. [58] R. Cook, Improved two-dimensional finite element, J. Struct. Div. 100 (1974) 1851–1863. [59] G.R. Liu, H. Nguyen-Xuan, T. Nguyen-Thoi, A variationally consistent αfem (vcαfem) for solution bounds and nearly exact solution to solid mechanics problems using quadrilateral elements, Internat. J. Numer. Methods Engrg. 85 (4) (2011) 461–497. [60] Ansys, http://www.ansys.com. [61] cfmesh, https://sourceforge.net/projects/cfmesh/, Accessed: 18-12-02.