Highly accurate smoothed finite element methods based on simplified eight-noded hexahedron elements

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Contents lists available at ScienceDirect

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

Highly accurate smoothed finite element methods based on simplified eight-noded hexahedron elements Y.H. Li a, R.P. Niu a,∗, G.R. Liu b a b

College of Mathematics, Taiyuan University of Technology, Taiyuan, Shanxi 030024, China School of Aerospace Systems, University of Cincinnati 2851 Woodside Dr, Cincinnati, OH 45221, USA

a r t i c l e

i n f o

Keywords: Eight-noded hexahedron element Coordinate mapping Simplified integration technique S-FEM-H8

a b s t r a c t Compared with the tetrahedron elements, hexahedron elements are preferred for their high accuracy. However, coordinate mapping required in the hexahedron elements of FEM formulation costs huge running time, leading to poor performance. Besides, the high quality of Jacobian matrix and mesh is required, which affects the accuracy of the strain results greatly. In order to solve these problems, we propose a novel simplified integration technique based on the smoothed finite element method (S-FEM) for the eight-noded hexahedron elements, where coordinate mapping is not demanded. The proposed new S-FEM-H8 models include simplified NS-FEM-H8 (using node-based smoothing domains) and simplified FS-FEM-H8 (using face-based smoothing domains). In the work, we divide a quadrilateral surface segment of a smoothing domain into two triangular sub-segments, so that the strain-displacement matrix can be calculated using a simple summation in the S-FEM theory instead of the integration in FEM. Then we conduct the Gauss integration scheme in each triangular surface sub-segment in order to avoid the coordinate mapping required in quadrilateral surface segments. The rest solving algorithm is the same as the standard S-FEM. Intensive numerical examples demonstrate that the simplified S-FEM-H8 possess the following features: (1) The strain energy of simplified NS-FEM-H8 is an upper bound of the exact solutions; (2) The simplified NS-FEM-H8 can overcome the volume locking problems for incompressible materials; (3) The method of dividing boundary surface into two triangular surfaces in smoothing domain keeps nearly the same accuracy as the standard S-FEM-H8.

1. Introduction Since the late 1950s, the finite element method (FEM) [1–3], as a tool for modeling and simulating practical problems, is an important method in computational mechanics. The FEM was originally applied in engineering science and technology to simulate and solve physical problems such as engineering mechanics, thermals, and electromagnetics. In recent years, FEM has become a powerful and versatile numerical analytic technique for practical engineering problems. However, when using FEM we encounter some problems, such as overly stiff issues and the significant loss of precision caused by distortions of meshes. In the process of searching effective numerical methods, it is found that mesh-free methods [4–7] can overcome some shortcomings of FEM when solving problems of large deformations, advanced materials, complex geometry, nonlinear material behavior, discontinuities and singularities. However, the methods bring other problems, such as determining the proper shape parameter and demanding more computing resources. A recent review on this can be found in [8]. In recent years, a new numerical method-smoothed finite element method

∗

(S-FEM) [9,10] was presented by G.R. Liu to overcome the disadvantages of the standard FEM, which combines FEM and mesh-free techniques. S-FEM always produces models that are softer than FEM and even softer than the exact model, which is used for bounding the range of the exact solutions for a practical problem. S-FEM can obtain relatively accurate strain solutions for the distorted mesh while FEM cannot work. Chang-Kye Lee [11] proposed a robust and effective smoothed finite element method (S-FEM), which can effectively overcome the volume locking phenomenon for nearly-incompressible materials and is significantly more robust than the standard FEM on highly distorted meshes. In addition, compared with the traditional weak form method (FEM), SFEM reduces the consistency requirement of the shape function because it is based on the W2 form [12–14]. Different smoothing finite element models will be generated by dividing the various forms of the smoothing domains. As presented in [10], different S-FEM models have different properties. For example, ES-FEM (Edge-based Smoothed Finite Method) [15–17] possesses a close-to-exact stiffness, which gives super-accurate results for static problems and stable, accurate results for dynamic problems. NS-FEM (Node-based Smoothed Finite Method) can obtain the

Corresponding author. E-mail address: [email protected] (R.P. Niu).

https://doi.org/10.1016/j.enganabound.2019.03.020 Received 5 December 2018; Received in revised form 20 March 2019; Accepted 21 March 2019 Available online 30 April 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

upper bound of the strain energy solutions [18,19]. Therefore, S-FEM has been used for solving engineering problems in different fields. In practical problems, engineers prefer to use linear elements [20], such as linear triangular (T3) and tetrahedron (T4) elements, since the elements can be easily obtained using automatically dividing programs. However, the poor accuracy of the strain solutions is the biggest problem for linear elements. Although FEM based on bilinear elements [21], such as quadrilateral element (Q4) and eight-noded hexahedron element (H8), can overcome the shortcoming of the poor precision in linear elements [22–25], the coordinate mapping and the strict requirement of the quality of the mesh lead to the limitation of the bilinear elements. To solve these problems, we present a novel simplified integration technique for H8 elements in S-FEM to avoid the coordinate mapping. In our work, we divide a quadrilateral surface segment of a smoothing domain into two triangular surface sub-segments. Then we use one Gauss point for each triangular surface sub-segment where the boundary integral is conducted in S-FEM, which makes the coordinate mapping removed. At the same time, the calculation of the Jacobian matrix required for the bilinear isoparametric element is avoided. According to the division rule, the proposed S-FEM-H8 includes simplified NS-FEMH8 and simplified FS-FEM-H8. The paper is outlined as follows. In Section 2, the idea of SFEM-H8 is briefly reviewed. In Section 3, a simplified integration technique is presented. Section 4 gives an overview on the patch test. Some numerical examples are analyzed using simplified NS-FEM-H8 and FS-FEM-H8 in Section 5. Some conclusions are provided in the last section.

a

: field nodes

⎧ ⎫ ⎪𝐝1 ⎪ } ∑ { ⎪𝐝 ⎪ 𝐍𝐼 (𝐱)𝐝𝐼 = 𝐍1 (𝐱)𝐍2 (𝐱) ⋯ 𝐍𝑁𝑛 (𝐱) ⎨ 2 ⎬ = 𝐍(𝐱)𝐝, 𝐮̄ (𝐱) = ⋮ ⎪ 𝐼=1 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⎪ ⎪𝐝𝑁 ⎪ 𝐍 ⎩ 𝑛⎭

(4)

where x = (x1 ⋅⋅⋅xd }T and 𝐝𝐼 is the displacement vector at the node I in the S-FEM models. Substituting Eq. (4) into Eq. (3), we have the following form: ( ) 1 1 ̄ ⋅ 𝐝̄. 𝐋 𝐮̄ (𝐱)𝑑 Γ = 𝑠 𝐋𝑛 𝐍(𝐱)𝑑 Γ ⋅ 𝐝̄ = 𝐁 (5) 𝛆̄ 𝐱𝑘 = 𝑠 𝐴𝑘 ∫ 𝑛 𝐴𝑘 Γ𝑠𝑘

The smoothed strain at any point xk can be calculated using the strain ̄𝐼 smoothing technique [26]. The smoothed strain-displacement matrix 𝐁 can be obtained for the 3D problems by Eq. (6).

(1)

Ω𝑠𝑘

⎡𝑏̄ 𝐼𝑥 ⎢ 0 ⎢ 0 ̄𝐁𝐼 = 1 𝐋 (𝐱)N(𝐱)𝑑Ω = ⎢ ⎢ 0 𝑉𝑘𝑠 ∫ 𝑛 ⎢𝑏̄ Ω𝑠𝑘 ⎢ 𝐼𝑧 ⎣𝑏̄ 𝐼𝑦

where Ω𝑠𝑘 is a smoothing domain defined in the local vicinity of xk , Ld 𝐇10,ℎ (Ω; 𝐑𝑑 )and

: center point po of face

𝑁𝑛

In the S-FEM models, a smoothed strain field can be obtained via a proper modification of the compatible strain field as follows: ⌢ ( ) 𝐋𝑑 𝐮̄ (𝐱)𝑊 𝐱𝑘 − 𝐱 𝑑Ω,

: mid-edge point

𝐮̄ ∈ 𝐇10,ℎ (Ω; 𝐑𝑑 ) is assumed as:

2.1. Smoothed strain field

∫

: center point of element

Fig. 1. (a) the face-based smoothing domains; (b) the node-based smoothing domains; (c) the edge-based smoothing domains.

2. Brief on the theory of smoothed finite element methods

( ) 𝛆 𝐱𝑘 =

c

b

⌢

is the matrix of differential operators, 𝐮̄ ∈ 𝑊 (𝐱𝑘 − 𝐱) is a weight function associated with xk and satisfies the unity property. We use the following the simplest form generally: { ⌢( ) 1∕𝐴𝑠𝑘 , 𝐱 ∈ Ω𝑠𝑘 ∪ Γ𝑠𝑘 𝑊 𝐱𝑘 − 𝐱 = , (2) 0, 𝐱 ∉ Ω𝑠𝑘 ∪ Γ𝑠𝑘

0 𝑏̄ 𝐼𝑦 0 𝑏̄ 𝐼𝑧 0 𝑏̄ 𝐼𝑥

0 ⎤ 0 ⎥ ⎥ 𝑏̄ 𝐼𝑧 ⎥ , 𝑏̄ 𝐼𝑦 ⎥ ̄𝑏𝐼𝑥 ⎥ ⎥ 0 ⎦

(6)

where 𝑛𝑠

𝐴𝑠𝑘

in which is the area of the kth smoothing domain. Therefore, using the Green’s divergence theorem and substituting Eq. (2) to Eq. (1), it can become:

Ω ( ) 1 1 ∑ 𝑏̄ 𝐼ℎ = 𝑠 𝑛 (𝐱)𝑁𝐼 (𝐱)𝑑Ω = 𝑠 𝑛 𝑁 𝐱𝐺 𝐴𝑆𝑢𝑟𝑓 , (ℎ = 𝑥, 𝑦, 𝑧) 𝑝 𝑉𝑘 ∫ ℎ 𝑉𝑘 𝑝=1 ℎ,𝑝 𝐼 𝑝

( ) 𝛆 𝐱𝑘 =

and 𝑛𝑠Ω is the total number of boundary surfaces Ω𝑠𝑘 and 𝐱𝑝𝐺 is the Gauss point of the boundary surfaces of Ω𝑠𝑘 , whose area and outward unit nor-

∫

⌢ ( ) 𝑊 𝐱𝑘 − 𝐱 𝐋𝑑 𝐮̄ (𝐱)𝑑Ω

Ω𝑠𝑘

mal are denoted as 𝐴𝑆𝑢𝑟𝑓 and nh,p , respectively. The above surface in𝑝 tegration along Ω𝑠𝑘 can be carried out using the Gauss quadrature technique [27]. 𝑏̄ 𝐼ℎ is computed along the surface segments of the smoothing domain according to Eq. (7), where the surface segments include the quadrilateral and triangular surface segments for S-FEM-H8, as shown in Fig. 1. For quadrilateral surface segments, four Gauss points are needed for Gauss integral technique. When the standard procedure is used to create the shape functions [28] for Q4 elements, we will encounter the problems with the compatibility issue. Thus, the formulation of quadrilateral surface segments requires a coordinate mapping procedure which is time-consuming. In addition, as we know, the Jacobian matrix must be calculated when bilinear isoparametric elements are used for calculating the stiffness matrix. Hence, the quality of meshes will greatly influence the accuracy of solutions of S-FEM-H8.

) 1 = − 𝐋𝑑 𝑊 𝐱𝑘 − 𝐱 𝐮̄ (𝐱)𝑑Ω + 𝑠 𝐋 (𝐱)𝐮̄ (𝐱)𝑑Γ ∫ 𝐴𝑘 ∫ 𝑛 ⌢

(

Ω𝑠𝑘

=

Γ𝑠𝑘

1 𝐋 𝐮̄ (𝐱)𝑑Γ, 𝐴𝑠𝑘 ∫ 𝑛

(7)

Ω𝑠𝑘

(3)

Γ𝑠𝑘

where Ln (x) is the matrix of components of the outward normal vector on boundary Γ𝑠𝑘 . 2.2. Smoothed strain-displacement matrix Following exactly the same way as in the standard FEM, the solutions of S-FEM model in terms of a continuous displacement field function 166

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

However, for a triangular surface segment, one Gauss point is sufficient for linear elements. Using the physical Cartesian coordinate system Oxy directly, the shape functions have the following forms: ) 1 ( 𝑁𝑗 = (8) 𝑎𝑗 + 𝑏𝑗 𝑥 + 𝑐𝑗 𝑦 , 6𝐴𝑆𝑢𝑟𝑓 𝑝

5

1

a

6

2

14

8 7

𝛿 𝐮̄ 𝐛𝑑Ω −

∫

̄ 𝑑Ω ⋅𝐝̄ = 𝐁 𝐜𝐁

Ω

3(x3,y3)

4(-1,1)

⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟ ̄ 𝐊

3(1,1)

𝛿 𝐮̄ 𝐭 𝑑Ω = 0,

ξ 1(x1 ,y1 )

2(x2 ,y2 ) x

1(-1,-1)

Physical coordinates

2(1,-1)

Natural coordinates

2x2 Gauss points

(9) Fig. 3. The coordinate mapping for a quadrilateral surface segment in FS domains.

Γ𝑡

∫

2

1

η

where Ns is the number of smoothed domains used for the entire problems domain, 𝑉𝑘𝑠 is the volume of the kth smoothing domain, 𝐮̄ ∈ 𝐇10,ℎ (Ω; 𝐑𝑑 ), c is a matrix of material constants, b ∈ L2 (Ω; Rd ) is the external body force applied over the problem domain, and t is the external traction force applied on the boundary Γt of the problem domain. Substituting Eqs. (4) and (5) into (9), we have the standard discretized algebraic system equations: ̄𝑇

3

y

⎞ ⎛ ⎞ ⎟ ⎜ 1 ⎟ 𝐋𝑛 (𝐱)𝛿 𝐮̄ (𝐱)𝑑Ω⎟ 𝐜⎜ 𝑠 𝐋𝑛 (𝐱)𝛿 𝐮̄ (𝐱)𝑑Ω⎟ ⎟ ⎜ 𝑉𝑘 Ω∫𝑠 ⎟ ⎠ ⎝ ⎠ 𝑘

Ω

∫

12

𝑇

⎛

∫

4

11

The smoothed Galerkin weak form can be written as [10] ⎜ 1 𝑉𝑘𝑠 ⎜ 𝑠 ⎜ 𝑉𝑘 Ω∫𝑠 𝑘=1 ⎝ 𝑘

6

Fig. 2. The face-based smoothing domains: (a) A FS domain for an interior face (b) A FS domain for a boundary face.

4(x4 ,y4 )

−

5 9

3

2.3. Smoothed Galerkin weak form

𝑁𝑠 ∑

7

8

b

10

13

4

where the subscript j varies from 1 to 3, 𝐴𝑆𝑢𝑟𝑓 is the area of the trian𝑝 gular surface segment, and aj , bj , cj are constants. It is clear that any coordinate mapping is not demanded in Eq. (8). In consideration of the difference of the two types of segment surfaces, we try to employ the advantages of the triangular surface segment to simplify the computation of quadrilateral surface segments.

9

𝑇

𝐍 𝐛𝑑Ω +

∫

𝑇

𝐍 𝐭𝑑Γ .

5

y

6

4(x4 ,y4 )

3

1(x1,y1)

3(x3 ,y3 )

2(x2,y2) x

(10)

2

1

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐟̄

Physical coordinates

Fig. 4. A simplified face-based smoothing domain (FS-FEM-H8). (a) a FS domain; (b) two triangular sub-surfaces obtained from a quadrilateral surface.

The above equation can be written as ̄ 𝐝̄ = 𝐟̃, 𝐊

b

9 4

Γ𝑡

Ω

7

8

a

(11)

a

̄ is the smoothed stiffness matrix, 𝐝̄ is the vector of nodal diswhere 𝐊 placements for the all nodes in the S-FEM model and 𝐟̃ is the force vector.

7

6 5

8

4

6

7

8

3

3. Simplified S-FEM-H8 formulations

b

5

2 1

15

In the S-FEM models, a mesh of elements is required, which can be created in exactly the same manner as in the standard FEM. Based on the background mesh, the problem domain Ω is divided into Ns “nonover𝑁𝑠 lap” and “no-gap” smoothing domains Ω𝑠𝑘 such that Ω= ∪𝑘=1 Ω𝑠𝑘 and 𝑠 𝑠 Ω𝑖 ∩ Ω𝑗 = ∅, 𝑖 ≠ 𝑗. In an S-FEM model, the smoothing domain usually consists of a num𝑛𝑠 ber of nonoverlapping and no-gap sub-smoothing domains: Ω𝑠𝑘 = ∪𝑞=1 𝑠 𝑠 𝑠 𝑠 Ω𝑘,𝑞 and Ω𝑘,𝑖 ∩ Ω𝑘,𝑗 = ∅, 𝑖 ≠ 𝑗, where Ω𝑘,𝑞 is the qth sub-smoothing domain and ns is the number of sub-smoothing domains.

13

9

3

11

2 14 10

4

12

1

Fig. 5. A node-based smoothing domain. (a) a NS domain; (b) a quarter of the NS domain.

3.1. 3D simplified FS-FEM-H8

in order to solve compatibility problems when building shape functions, we need to perform coordinate transformation for 2 × 2 Gauss points to calculate the element stiffness matrix in Eq. (7). At the same time, the Jacobian matrix must be calculated because of the coordinate mapping. Compared with the quadrilateral surface segments, the triangular surface segments just need one Gauss point, which is the center point of the triangle, to compute the element stiffness matrix. Therefore, coordinate mapping and Jacobian matrix can be omitted. Taking Fig. 2(b) as an example, when 𝑏̄ 𝐼ℎ is computed along the surface of the smoothing domain according to Eq. (7), the quadrilateral segment 1-2-3-4 should be transferred from the physical coordinate to the natural coordinate, as shown in Fig. 3. Then the positions of the

In 3D problems, the face-based smoothing domain is the simplest smoothing domain. A face-based smoothing domain is created generally by joining the center points of the two adjacent elements to the four nodes of the face. As shown in Fig. 2(a), we separately connect the center points 13 and 14 of the two adjacent background elements to the nodes of the face 5-6-7-8. For boundary faces, we just need connect the center point 9 of the background element to the four nodes of the boundary face 1-2-3-4 to form a quadrangular pyramid, as shown in Fig. 2(b). From Fig. 2, it is easily found that the surfaces of smoothing domains consist of quadrangles and triangles. For quadrilateral surface segments, 167

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

3.2. 3D simplified NS-FEM-H8

6

7 8

Compared with face-based smoothing domains, a node of H8 element is associated with more elements, faces and edges. Thus node-based smoothing domain is more complex. A standard NS domain based on H8 is shown in Fig. 5(a). Taking one eighth of the NS domain as an example, we should successively connect the center point 9 of the element (1-2-3-4-5-6-7-8), two centroids (14 and 15) of the faces (1-2-3-4 and 1-5-8-4) and the middle point 12 of the edge (1-4) to create a quadrilateral surface segment (9-15-12-14), as shown in Fig. 5(b). We construct the NS domain for node 1 by repeating the procedure through all the elements, the faces and the edges including the node. It is clearly seen that the node will be wrapped up by a number of quadrilateral surface segments shown in Fig. 5(a). Similar to FS-FEM-H8, we divide quadrilateral surface segments to triangular sub-segments. Because all the surfaces of NS domains are quadrilateral, the simplified NS domains are a little more complex, as shown in Fig. 6. However, comparing with the coordinate mapping and solving the Jacobian matrix in quadrilateral surface segments (9-15-1214, 9-15-11-13, 9-14-10-13), we just need to calculate the center points of the triangular sub-segments (9-15-12, 9-12-14, 9-15-13, 15-11-13, 914-10, 9-10-13), which is very simplify and fast. Taking the Fig. 5(a) as an example, we study the internal node 1 of the element. It is obvious that the smoothing domain based on the node 1 has 24 smoothing sub-segments, and four Gauss points need to be obtained for each subsmoothing domain. In the process of finding the smoothed strain matrices, we need to calculate 96 Gauss points, coordinate mapping 24 times and Jacobian matrices 24 times for the smoothing domain based on node 1. However, we just need 48 Gauss points, no coordinate mapping and no Jacobian matrices in our work. Therefore, the simplified NS-FEM-H8 can greatly improve the computational efficiency of the program and reduce the computational complexity.

5 13

9 15

11 2

3

14 10

4

1

12

Fig. 6. Simplified node-based smoothing domains (NS-FEM-H8).

7

6 5

8 11 3 4

10 9

2 1

Fig. 7. A edge-based smoothing domain.

four Gauss points need to be calculated in the square. In our present method, we divide the quadrilateral segment 1-2-3-4 into two triangular sub-segments 1-2-4 and 2-3-4, as shown in Fig. 4. It is clearly that the coordinate mapping from physical coordinates to natural coordinates is not needed. Next, in each triangular sub-segment, we choose one Gauss point to calculate the boundary integral for 𝑏̄ 𝐼ℎ .

3.3. 3D ES-FEM-H8 For 3D problems, edge-based smoothing domain is formed by sequentially connecting the center point 9 of the element (1-2-3-4-5-6-

a

b

Fig. 8. The meshes with 64 nodes and 27 H8 elements used in the patch test: (a) eight-noded hexahedron regular elements (b) eight-noded hexahedron irregular elements. 168

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Fig. 9. The problem domains, problem settings and a H8 mesh. (a) A 3D cubic cantilever subjected to a uniform pressure on the top surface; (b) A mesh with H8 elements.

Table 1 The relative error of the displacement obtained by standard patch tests of simplified S-FEM-H8 models. Relative error of displacement 𝛼 ir Regular element Irregular element

0 0.1 0.2 0.3 0.4

Simplified FS-FEM-H8

Simplified NS-FEM-H8

5.0114e−17 8.5006e−17 1.0298e−16 4.3850e−16 2.3744e−16

1.3454e−16 1.5764e−16 1.5352e−16 1.2493e−16 1.1242e−16

Pass Pass Pass Pass Pass

Pass Pass Pass Pass Pass

7-8) to the two nodes of this edge (1-5), the centroids (10,11) of the faces (1-2-6-5, 1-5-8-4) containing the edge to the two nodes (1,5) of this edge and the centroids (10,11) of the faces containing the edge to the center point 9 of the element, as shown in Fig. 7. The mesh plotted by gray dotted lines is an eight-noded hexahedron background mesh, over which ES domains are created, as plotted in black lines. We can easily find that the ES smoothing surface segments (1-9-11, 1-10-9, 5-9-11, 5-10-9) are all triangles for standard ES-FEM-H8. Therefore, it is not necessary to simplify the standard ES-FEM-H8, which is directly employed in our numerical examples.

Fig. 10. Convergence of the strain energy solutions obtained using different methods for the 3D cubic cantilever problem subjected to a uniform pressure.

4. Standard patch test for 3D problems Firstly, we test our present algorithm by the patch test proposed by MacNeal and Harder [29]. The patch test examines the stabilization and accuracy of the present method under linear displacements imposed on the Dirichlet boundary conditions. Fig. 8(a) shows the cubic patch with dimensions of 3 by 3 by 3. Now we study the effect of mesh distortions on the solutions. The irregular spatial discretization of a cubic patch using 27 eight-noded hexahedron elements is shown in Fig. 8(b). The parameters are taken as E = 100, m = 0.3 and linear displacement field is given by

are computer-generated random numbers between −1 and 1. 𝛼 ir ∈ [0.0, 0.5) is the degree of distortion of the control elements. When 𝛼 ir =0.0, the interior node locates at the center point of the cubic patch, and when 𝛼 ir > 0.0, the interior node moves randomly inside the cubic patch. The larger the value of 𝛼 ir , the more irregularly the shape of elements are generated. When 𝛼 ir =0.49, the interior node is almost reaching to the surface of the cube. The following error norm in displacements is used in our paper to examine the computed results:

𝑢 = 𝑥, 𝑣 = 𝑦, 𝑤 = 𝑧.

√ ) √ ∑𝑁𝑑𝑜𝑓 ( 2 2 √ 𝑢𝑖 − 𝑢ℎ𝑖 √ 𝑖=1 𝑒𝑑 = √ × 100%, ∑𝑁𝑑𝑜𝑓 2 𝑢𝑖 𝑖=1

(12)

We can obtain the Dirichlet boundary conditions from the analytical solutions. The coordinates of interior nodes are perturbed as follows: 𝑥 = 𝑥 + 𝑟𝑥 𝛼𝑖𝑟 Δ𝑥

where ui is the exact solution and 𝑢ℎ𝑖 is the numerical solution. From Table 1, it is found that the simplified FS-FEM-H8 and NSFEM-H8 can pass the standard first-order patch test within the machine precision regardless of the value of 𝛼 ir ∈ [0.0, 0.5). This shows that the simplified FS-FEM-H8 and NS-FEM-H8 can work well with severely distorted meshes.

𝑦 = 𝑦 + 𝑟𝑦 𝛼𝑖𝑟 Δ𝑦, 𝑧 = 𝑧 + 𝑟𝑧 𝛼𝑖𝑟 Δ𝑧

(14)

(13)

where x, y, z are the centroid coordinate of the cubic patch, Δx, Δy, Δz are the length along the x-, y- and z-directions, respectively, and rx , ry rz 169

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Table 2 Coordinate mappings and the calculations of Jacobian matrix using the standard NS-FEM-H8 model for 10 × 10 × 10 cube.

Number of smoothing sub-surfaces Number of Gauss points Number of coordinate mapping Calculate the number of Jacobian matrices The number of nodes of this type

Point A

Point B

Point C

Point D

10 × 10 × 10 cube

24 96 24 24 729

16 64 16 16 486

10 40 10 10 108

6 24 6 6 8

29,784 119,136 29,784 29,784 1331

Note: Point A is the internal node; point B is the node on the boundary surface; point C is the node on the boundary line; point D is the node that is all exposed.

a

which make the algorithm of S-FEMs simpler, more efficient and more practical. Fig. 10 also shows the remarkable characteristic [31] of the ES-FEM-H8: it is not only softer than the hard FEM model, but also harder than the soft NS-FEM-H8 model with a close-to-exact stiffness. Taking a 10 × 10 × 10 cube background mesh as an example, we require four Gauss points and perform a coordinate mapping and calculate a Jacobian matrix for each quadrilateral smoothing surface in the standard NS-FEM-H8 model. As shown in Fig. 11 and Table 2, for 729 internal nodes (Fig. 11(a)), there are 24 quadrilateral smoothing sub-surfaces for each smoothing domain; for 486 nodes (Fig. 11(b)) on the boundary surfaces, there are 16 smoothing sub-surfaces for each smoothing domain; for 108 nodes (Fig. 11(c)) on the boundary lines, each smoothing domain has 10 smoothing sub-surfaces; for 8 nodes (Fig. 11(d)) that are all exposed to outside, each smoothing domain has 6 smoothing subsurfaces. Therefore, for such a simple cube with only 1331 nodes, we need to perform coordinate mapping and the calculations of Jacobian matrix 29,784 times and calculate 119,136 Gauss points, which seriously affect the computational efficiency of the program, greatly consume huge computation resources and spend a lot of running time. In fact, we just need a half of the Gauss points that is easier to obtain than standard methods to complete the same task. The conclusion verifies the efficiency of the simplified NS-FEM-H8 for the stiffness matrix calculation. Then, the displacement and its convergence are analyzed. The vertical displacement solutions v of FEM-H8, S-FEM-H8 and the reference solutions are plotted in Fig. 12. The reference solutions are calculated using FEM with the 15,625 elements through ABAQUS®. It is easily observed that all the computed displacements using different methods are in good agreement with the reference solutions. With the refinement of the mesh, the more accuracy results are performed and the approximate solutions greatly approach the reference solutions. Particularly, when the large size of the mesh h reaches 1/25, the solutions are almost completely coincided with the reference solution. Besides, it is again found that the FEM-H8 model is very stiff, while the simplified NS-FEM-H8 model is very soft compared with the reference model, which can bound the range of the solutions for practical problems. Besides, we plot the distribution of displacement solutions v obtained using FEM-H8, simplified FS-FEM-H8, simplified NS-FEM-H8 and the reference solutions in Fig. 13. From these figures, we find that the solutions of simplified FS-FEM-H8 and simplified NS-FEM-H8 agree well with the reference solutions gained with ABAQUS®.

b

A

A

c

d

A

A Fig. 11. Smoothing domain at different locations of node A in NS-FEM-H8 model. (a) Point A is the internal node; (b) Point A is the node on the boundary surface; (c) Point A is the node on the boundary line; (d) Point A is the node that is all exposed.

5. Numerical example 5.1. 3D cubic cantilever As shown in Fig. 9(a), consider a three-dimensional cantilever block in which the upper surface is uniform pressure and the left surface is fixed. The parameters used in the example are E = 1.0 and v = 0.2. Fig. 9(b) illustrates the discretization based on a regular mesh of H8 elements. Fig. 10 shows that the convergence of the strain energy solutions obtained using different methods for the 3D cubic cantilever problem subjected to a uniform pressure. The exact solution for this problem is unknown. However, Almeida Pereira [30] gave a reference solution of strain energy of 0.950930 by combining hexagonal elements with Richardson’s extrapolation method. The figure confirms that the upper bound property on the strain energy of simplified NS-FEM-H8 and the lower bound property of simplified FS-FEM-H8. For easy comparisons, we list the solutions of the standard ES-FEM-H8, the standard FS-FEMH8, the standard NS-FEM-H8, FEM-H8 and FEM-T4 for this 3D problem. It is seen that the FEM-H8 model is more accurate than the FEM-T4 model in general. We obviously find that the strain energy solutions of NS-FEM-H8 and FS-FEM-H8 with a new simplified integration technique are almost the same as the solutions of the standard S-FEM-H8,

5.2. 3D lame problem In this example, we study the three-dimensional Lame problem consisting of a hollow sphere with inner radius a = 1 m, outer radius b = 2 m, which is subjected to an internal pressure P = 1 N/m2 . Since the sphere is symmetric about the three axes in the problem, only one-eighth of the sphere will be studied instead of the whole problem domain, as shown in Fig. 14. Symmetry conditions are imposed on the three symmetric planes. 170

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

h=1/4

h=1/5

h=1/10

h=1/25

Fig. 12. Comparison of displacement solutions v obtained using different methods and the reference solution for the 3D cubic cantilever problem (h is the maximum length of a cube element in the background grid).

The exact stresses [23] of the problem are given as: ] [ 𝑃 𝑎3 𝑟 𝑏3 𝑢𝑟 = ( ) (1 − 2𝑣) + (1 + 𝑣) 3 , 2𝑟 𝐸 𝑏3 − 𝑎3 ( 3 ) ( 3 ) 3 3 3 3 𝑃𝑎 𝑏 − 𝑟 𝑃 𝑎 𝑏 + 2𝑟 𝜎𝜏 = ( ) , 𝜎𝜃 = ( ) , 𝑟3 𝑎3 − 𝑏3 2𝑟3 𝑏3 − 𝑎3

Then, the distributions of the radial displacement, radial and tangential stresses using FEM-H8, simplified FS-FEM-H8, simplified NS-FEMH8 and standard ES-FEM-H8 compared with the analytical solutions are presented in Fig. 16, where h is the characteristic length It shows that when h reaches to 1/10 in which the mesh is relatively dense and fine, the approximate radial displacement, radial and tangential stress are very close to the exact solutions. Hence, the numerical solutions of the above methods all converge to the analytical solutions with the increase of DOFs. Next, we analyze a volumetric locking problem using introduced S-FEM-H8 models. We set Poisson’s ratio of the hollow sphere changing from 0.4 to 0.4999999. Fig. 17 plots the errors in the displacement norm obtained different numerical methods using eight-noded hexahedron elements. For the almost incompressible problems, the results show that simplified NS-FEM-H8 is naturally immune from volumetric locking, while standard FEM-H8, ES-FEM-H8 and simplified FS-FEM-H8 are

(15)

(16)

where r is the radial distance from the center of the sphere to any point on the surface of the sphere. The material parameters of the problem are E = 1.0 and v = 0.3. The exact strain energy of the problem is 6.3e-04. Fig. 15 shows that the convergence of the strain energy solutions obtained using different methods for the hollow sphere subjected to inner pressure. We can easily find that the strain energy solutions of simplified FS-FEM-H8 and NS-FEM-H8 are both more accurate than that of FEMH8. Besides, it confirms the upper bound property on the strain energy of simplified NS-FEM-H8 for this 3D problem, which illustrates that simplified NS-FEM-H8 keeps the properties of NS-FEM-H8. The same phenomenon can be found in simplified FS-FEM-H8. 171

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

a

b

c

d

y

z x

x

Fig. 13. Distributions of displacement solutions v and reference solutions obtained using different methods for the 3D cubic cantilever problem. (a) the displacement solutions of FEM-H8; (b) the displacement solutions of simplified FS-FEM-H8; (c) the displacement solutions of simplified NS-FEM-H8; (d) the displacement solutions of the reference solutions. Table 3 The relative error of the displacement obtained using FEM-H8 and simplified S-FEMH8 models by the hollow sphere subjected to inner pressure. Relative error of displacement 𝛼 ir regular element irregular element

0 0.1 0.2 0.3 0.4 0.49

FEM-H8

Simplified FS-FEM-H8

Simplified NS-FEM-H8

0.0064 0.0092 0.0240 0.0557 0.0993 0.1598

0.0039 0.0076 0.0222 0.0554 0.0893 0.1162

0.0140 0.0198 0.0449 0.0780 0.1124 0.1459

subjected to the volumetric locking [32], resulting in a drastic accuracy loss in the numerical solutions. Finally, we will study the stability [33] of the simplified method using different irregular coefficients. Fig. 18 shows that the regular mesh and irregular mesh of the hollow sphere using 99 H8 elements. When the mesh is extremely deformed, the shape of the H8 element varies greatly, but it is still a hexahedral element, and each face is a quadrilateral.

The relative displacement errors obtained using different methods are listed in Table 3 for meshes of irregular coefficients which are changed from 0 to 0.49 for the hollow sphere model. It is observed that the relative errors of the simplified FS-FEM-H8 are smaller than that of the FEM-H8, while the errors of the simplified NS-FEM-H8 are larger. When the irregular coefficient is gradually increased, the error will increase, but the amplitude of the increase is not very large, especially for

172

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Fig. 14. The problem domain and a H8 mesh. (a) a hollow sphere problem; (b) one-eighth model discretized using eight-noded hexahedron elements. Table 4 Strain energy (Nm) obtained using different methods for a 3D spanner problem subjected to a uniform pressure.

DOFs FEM-H8 simplified FS-FEM-H8 ES-FEM-H8 simplified NS-FEM-H8

Mesh1

Mesh2

Mesh3

Mesh4

Mesh5

Mesh6

856 × 3 0.4803 0.5065 0.5075 0.5524

1325 × 3 0.4979 0.5181 0.5207 0.5508

2298 × 3 0.5093 0.5227 0.5253 0.5480

2646 × 3 0.5124 0.5238 0.5255 0.5429

3577 × 3 0.5145 0.5244 0.5262 0.5422

5608 × 3 0.5197 0.5259 0.5271 0.5371

Table 5 Error in displacement norm using different methods for a 3D spanner problem subjected to a uniform pressure.

Mesh1 Mesh2 Mesh3 Mesh4 Mesh5

DOFs

FEM-H8

FS-FEM-H8

ES-FEM-H8

NS-FEM-H8

856 × 3 1325 × 3 2298 × 3 2646 × 3 3577 × 3

0.0406 0.0119 0.0077 0.0049 0.0032

0.0205 0.0040 0.0032 0.0019 0.0013

0.0214 0.0034 0.0028 0.0016 0.0011

0.0024 5.1732e−5 4.6592e−5 5.6238e−5 5.6250e−5

problem is unknown. It can be found that the simplified NS-FEM-H8 gives the upper bound solutions in the strain energy, while the strain energies of FEM-H8, ES-FEM-H8, and simplified FS-FEM-H8 are always less than the exact one. The errors of displacement norm is shown in Table 5 and plotted in Fig. 21. It is seen that the simplified NS-FEM-H8 stands out clearly. The error of displacement norm of the simplified NS-FEM-H8 is minimal in these methods. The convergence of the errors in energy norm is shown in Fig. 22. It is seen again that the simplified NS-FEM-H8 model stands out obviously compared with the other methods. The error of energy norm of the simplified NS-FEM-H8 is minimal in the four models.

Fig. 15. Convergence of the strain energy solutions obtained using different methods for the hollow sphere subjected to inner pressure.

NS-FEM-H8, which can still illustrate that these presented methods are relatively more stable when the meshes are deformed.

5.4. 3D pulley As shown in Fig. 23(a), we consider a three-dimensional pulley model whose material is aluminum, E = 76,000 MPa and v = 0.33. The back of the pulley model is fixed in this example and a uniform pressure P = 10 N/m2 is applied to the inner wall of the structure. Since the displacement of the pulley model is periodically symmetrical about the central axis, only a part of the structure (1/12 of the model) is considered in our numerical simulation, as shown in Fig. 23(b). The convergence of strain energy solutions and the computational efficiency of the simplified methods are studied in this example, compared with other methods.

5.3. 3D spanner Consider a 3D spanner with complex geometry problem, which is labeled in Fig. 19. The numerical parameters in the experiment are used: E = 210,000 Pa, v = 0.3 and the pressure is set as P = 1 N/m2 . The numerical results of strain energy are presented in Table 4 and plotted in Fig. 20 against the DOFs, revealing the convergence of the solutions of all the models used. The reference solution is obtained under very refined mesh by the ABAQUS® since the exact solution for this 173

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

(a) the radial displacement when h=1/4

(b) the radial displacement when h=1/10

(c) the radial and tangential stresses when h=1/4

(d) the radial and tangential stresses when h=1/10

Fig. 16. Radial displacement when (a) h = 1/4, (b) h = 1/10 and radial and tangential stresses when (c) h = 1/4, (d) h = 1/10 for the hollow sphere subjected to inner pressure.

Fig. 18. The meshes with 176 nodes and 99 H8 elements: (a) eight-noded hexahedron regular elements (b) eight-noded hexahedron irregular elements.

Fig. 17. Displacement norm versus different Poisson’s ratios for the hollow sphere subjected to inner pressure. 174

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Fig. 19. The problem geometry of the 3D spanner: (a) Geometric model of a 3D spanner; (b) a mesh with eight-noded hexahedron elements.

Fig. 22. Comparison of energy error norm using different methods for a 3D spanner problem. Fig. 20. Convergence of the strain energy solutions obtained using different methods for a 3D spanner problem subjected to a uniform pressure.

Fig. 24 depicts the convergence of strain energy solutions obtained using different methods for the 3D pulley problem. The reference solution is obtained by the ABAQUS® using a very fine mesh since the exact solution for this problem is unknown. From this figure, it is observed that the NS-FEM-H8 and simplified NS-FEM-H8 can obtain the upper bound solutions, and other methods can obtain the lower bound solutions. It is noted that the simplified methods preform a little well in strain energy solution, compared with their previous methods. Because the aim of the simplified methods is to improve the efficiency of S-FEM-H8, we next discuss the computational efficiency of the simplified methods for this example. Fig. 25 shows the CPU time versus degree of freedom and Fig. 26 shows the displacement error norms versus the CPU time, both using the different methods for a range of meshes. It is found that with the same sets of nodes, the computation time of the simplified methods is shorter than that of the previous method, as shown in Fig. 25. The CPU time of the simplified methods for the same accuracy is still relatively less, compared with the previous methods, as shown in Fig. 26. The efficiency is estimated using the following formulation, 𝐸𝑓 𝑓 𝑖𝑐 𝑖𝑒𝑛𝑐 𝑦 =

1 𝑟𝑒 × 𝑟𝑡

where re and rt represent, respectively, the displacement error norm and CPU time of the other methods compared with those of the FEM-H8 method. Table 6 lists the displacement error norms and the computational time obtained using different methods under the same mesh. Obviously,

Fig. 21. Comparison of displacement error norm using different methods for a 3D spanner problem.

175

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

Fig. 23. Pulley structure: (a) complete model (b) 1/12 of the model.

Fig. 26. Computation time for the solutions of the same accuracy using different methods for a 3D pulley problem.

Fig. 24. Convergence of the strain energy solutions obtained using different methods for a 3D pulley problem.

Table 6 Estimated computational efficiency of FEM-H8, Simplified FS-FEM-H8, FS-FEMH8, ES-FEM-H8, Simplified NS-FEM-H8 and NS-FEM-H8, measured in displacement norm error for the pulley problem using the same mesh (904 nodes).

1 2 3 4 5 6

Numerical method

Solution error

re

Estimated CPU time (s)

rt

Efficiency

FEM-H8 Simplified FS-FEM-H8 FS-FEM-H8 ES-FEM-H8 Simplified NS-FEM-H8 NS-FEM-H8

0.0315 0.0171

1.00 0.54

7.043556 2.999374

1.00 0.43

1.00 4.31

0.0172 0.0238 0.0190

0.55 0.76 0.60

4.749683 8.709973 17.834356

0.67 1.24 2.53

2.71 1.06 0.66

0.0191

0.61

24.834356

3.53

0.46

compared with the FS-FEM-H8, the error of the simplified FS-FEM-H8 is not significantly decreased, but the computational time is greatly reduced, and the computational efficiency is about 1.5 times that of the previous method for this example. Similarly, the computational efficiency of the simplified NS-FEM-H8 has been improved in comparasion with previous methods. Therefore, when the computational time is taken into account, the computational efficiency of the simplified FS-FEM-H8 and simplified NS-FEM-H8 can be significantly improved.

Fig. 25. The CPU time versus the degrees of freedom using different methods for a 3D pulley problem. 176

Y.H. Li, R.P. Niu and G.R. Liu

Engineering Analysis with Boundary Elements 105 (2019) 165–177

6. Conclusion

[10] Liu GR, ThoiTrung Nguyen. Smoothed finite element methods. CRC Press; 2010. p. 569–80. [11] Lee CK, Angela Mihai L, Hale JS, Kerfriden P, Bordas SPA. Strain smoothing for compressible and nearly-incompressible finite elasticity. Comput Struct 2017;182:540–55. [12] Liu GR. A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: part I theory. Int J Numer Methods Eng 2010;81(9):1093–126. [13] Liu GR. A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: part II applications. Int J Numer Methods Eng 2010;81(9):1127–56. [14] Liu GR, Zhang GY. Smoothed point interpolation methods-G space theory and weakened weak forms. World Scientific; 2013. [15] Liu GR, Nguyen-Thoi T, Lam KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses in solids. J Sound Vib 2009;320:1100–30. [16] He ZC, Li GY, Zhong ZH, et al. An edge-based smoothed tetrahedron finite element method (ES-T-FEM) for 3D static and dynamic problems. Comput Mech 2013;52(1):221–36. [17] Nguyen-Xuan H, Liu GR, Nguyen-Thoi T, Nguyen-Tran C. An edge-based smoothed finite element method (ES-FEM) for analysis of two dimensional piezoelectric structures. Smart Mater Struct 2009;18(12pp.):065015. [18] Nguyen-Thoi T, Liu GR, Nguyen-Xuan H. Additional properties of the node-based smoothed finite element method (NS-FEM) for solid mechanics problems. Int J Comput Methods 2009;6(4):633–66. [19] Zhang ZQ, Liu GR. Upper and lower bounds for natural frequencies: a property of the smoothed finite element methods. Int J Numer Methods Eng 2010. [20] Francis A, Ortiz-Bernardin A, Bordas SP, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. Int J Numer Methods Eng 2016;109(9):1263–88. [21] Natarajan S, Bordas SP, Ooi ET. Virtual and smoothed finite elements: a connection and its application to polygonal/polyhedral finite element methods. Int J Numer Methods Eng 2015;104(13):1173–99. [22] Puso MA. A highly efficient enhanced assumed strain physically stabilized hexahedral element. Int J Numer Methods Eng 2000;49(8):1029–64. [23] Fredriksson M, Ottosen NS. “Accurate eight-node hexahedral element. Int J Numer Methods Eng 2007;72(6):631–57. [24] Krysl P. Mean-strain eight-node hexahedron with stabilization by energy sampling. Int J Numer Methods Eng. 2014. http://dx.doi.org/10.1002/nme.4721 n/a–n/a. [25] Belytschko T, Bendeman LP. Assumed strain stabilization of the eight node hexahedral element. Comput Meth Appl Mech Eng 1993;105(2):225–60. [26] Zerroukat M, Power H, Chen CS. A numerical method for heat transfer problems using collocation and radial basis functions. Numer Methods Eng 1998;42(7):1263–78. [27] Chakraborty S, Natarajan S, Singh S, Roy Mahapatra D, Bordas SPA. Optimal numerical integration schemes for a family of polygonal finite elements with Schwarz–Christoffel conformal mapping. Int J Comput Methods Eng Sci Mech 2018:1–22. [28] Yue JH, Li M, Liu GR. Proofs of the stability and convergence of a weakened weak method using PIM shape functions. Comput Math Appl 2016;72:933–51. [29] MacNeal RH, Harder RL. A proposed standard set of problems to test finite element accuracy. Finite Elem Anal Des 1985;1(1):1–20. [30] Almeida Pereira OJB. Hybrid equilibrium hexahedral elements and super-elements. Commun Numer Methods Eng 2006;24(2):157–65. [31] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee CK, Zi G, Zhuang X, Rabczuk T. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. J Appl Math 2013;2013:1–12. [32] Nguyen-Xuan H, Bordas S, Nguyen-Dang H. Addressing volumetric locking and instabilities by selective integration in smoothed finite elements. Commun Numer Methods Eng 2008;25:19–34. [33] Ong TH, Heaney CE, Lee C-K, Liu GR, Nguyen-Xuan H. On stability, convergence and accuracy of bES-FEM and bFS-FEM for nearly incompressible elasticity. Comput Meth Appl Mech Eng 2015;285:315–45.

In the paper, we put forward simplified NS-FEM-H8 and FS-FEM-H8 models to avoid coordinate mapping and the calculations of Jacobian matrix for 3D problems, which can greatly improve the efficiency of the S-FEM-H8 algorithm. At the same time, the novel methods can maintain the high accuracy of the bilinear elements. Through numerical examples, some main conclusions are presented as follows: (1) The strain energy of simplified NS-FEM-H8 is an upper bound of the exact strain energy, while the solutions of simplified FS-FEMH8 and standard FEM-H8 are the lower bound of the exact strain energy. (2) For the almost incompressible problems, simplified NS-FEM-H8 has the characteristic of freeing from volume locking. (3) Compared with FEM-T4, it is easy to verify the high accuracy of simplified FS-FEM-H8 and simplified NS-FEM-H8. (4) The simplified S-FEM-H8 makes the standard S-FEM-H8 more efficient and more practical. Meanwhile, they have the same properties and nearly the same accuracy as that of standard S-FEM-H8. Acknowledgments The project is supported by the National Natural Science Foundation of China (Grant no.11771321) and the project of senior foreign experts (Grant no. GDW20181400420). Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.enganabound.2019.03.020. References [1] Liu GR, Quek SS. The finite element method - a practical course. 2nd ed. Elsevier (BH); 2013. [2] Hughes TJR. The finite element method. Linear static and dynamic finite element analysis. NJ: Prentice Hall; 1987. [3] Oliveira Eduardo R, De Arantes E. Theoretical foundations of the finite element method. Int J Solids Struct 1968;4:929–52. [4] Liu GR, Liu MB. Smoothed particle hydrodynamics: aa meshfree particle method, a bestseller of world scientific. [5] Liu GR. Mesh-free methods: moving beyond the finite element method. Boca Raton: CRC Press; 2009. [6] Chen JS, Wu CT, Yoon S, You Y. Astabilized conforming nodal integration for Galerkin meshfree method. Int J Numer Methods Eng 2001;50:435–66. [7] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79(3):763–813. [8] Liu GR. An overview on meshfree methods: for computational solid mechanics. Int J Comput Methods 2016;13(Issue 05):1630001. [9] Liu GR, Nguyen-Thoi T, Dai KY, Lam KY. Theoretical aspects of the smoothed finite element method (SFEM). Int J Numer Methods Eng 2007;71:902–30.

177