- Email: [email protected]
A three-dimensional adaptive analysis using the meshfree node-based smoothed point interpolation method (NS-PIM)
Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A three-dimensional adaptive analysis using the meshfree node-based smoothed point interpolation method (NS-PIM) Q. Tang a,b, G.Y. Zhang c,n, G.R. Liu d, Z.H. Zhong a, Z.C. He a,b a
State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, China Centre for Advanced Computations in Engineering Science (ACES), Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117576, Singapore c School of Mechanical and Chemical Engineering, The University of Western Australia. 35 Stirling Highway, WA 6009, Australia d School of Aerospace Systems, University of Cincinnati, Cincinnati, OH 45221-0070, USA b
a r t i c l e i n f o
abstract
Article history: Received 20 March 2009 Accepted 9 May 2011 Available online 7 June 2011
In this paper, a three-dimensional (3-D) adaptive analysis procedure is proposed using the meshfree node-based smoothed point interpolation method (NS-PIM). Previous study has shown that the NS-PIM works well with the simplest four-node tetrahedral mesh, which is easy to be implemented for complicated geometry. In contrast to the displacement-based FEM providing lower bound solutions, the NS-PIM possesses the attractive property of providing upper bound solutions in strain energy norm. In the present adaptive procedure, a novel error indicator is devised for NS-PIM settings, which evaluates the maximum difference of strain energy values among four nodes in each of the tetrahedral cells. A simple h-type local refinement scheme is adopted and coupled with 3-D mesh automatic generator based on Delaunay technology. Numerical results indicate that the adaptive refinement procedure can effectively capture the stress concentration and solution singularities, and carry out local refinement automatically. The present adaptive procedure achieves much higher convergence in strain energy solution compared to the uniform refinement, and obtains upper bound solution in strain energy efficiently for force driven problems. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Numerical methods Meshfree methods Error indicator Adaptive analysis Point interpolation Tetrahedral mesh generation
1. Introduction It is well known that adaptive analysis is important in the computational mechanics and engineering structural designs. The ultimate objective of adaptive analysis is to achieve desired high accuracy with minimum computational cost. Adaptive mesh refinement techniques along with proper error analysis have been well studied in the traditional FEM [1–4], but much less welldeveloped for meshfree method [5,6]. The meshfree node-based point interpolation method (NS-PIM or LC-PIM) [7] has been recently developed using the generalized smoothed Galerkin (GS- Galerkin) weak form [8] and PIM shape functions with a set of small number of nodes located in a local support domain. The PIM (polynomial PIM or radial PIM) shape functions possess Delta function property, which allows straightforward imposition of point essential boundary conditions. The use of node-based smoothing domains [7,9–11] provides the so-called ‘‘soft’’ effect to the discrete model and gives a number of good features. Compared with the linear FEM, NS-PIM can
n
Corresponding author. E-mail address: [email protected] (G.Y. Zhang).
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.05.019
obtain better accuracy and higher convergence rate, especially for stress results. Recently, Liu and Zhang [12] made a thorough theoretical study on the NS-PIM and obtained some important conclusions including the upper bound property. NS-PIM is easy to be implemented and suitable for adaptive analysis. An adaptive analysis using the NS-PIM has been conducted to certify solutions with exact bounds of strain energy for 2-D problems [13]. A reliable and efficient error indicator and the associated refinement strategy are crucial issues for the adaptive procedure, especially for 3-D problems. In general, two distinct types of procedures are currently available for deriving error indicators: the recovery based error indicator and the residual based error indicator. The recovery based error indicator was first introduced by Zienkiewicz and Zhu [1] in 1987, by constructing locally an improved solution from the approximation. The recovery processes play a critical role in the computation of this error indicator and there are many papers published on this recovery methodology [2,14,15]. Residual based error indicators make use of the residual of the numerical approximation, either explicitly or implicitly, which offers a very effective alternative [16–20]. Refinement schemes can also be classified into three categories: h-refinement, p-refinement and r-refinement. The h-refinement scheme changes the size of element in a localized fashion based on
1124
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
the error indicator. The p-refinement scheme is to increase the order of the polynomial, and the r-refinement keeps the total number of nodes unchanged but to adjust their position to obtain an optimal approximation [21,22]. Because the NS-PIM works well with lower order interpolation of displacement field, this work chooses to use the h-refinement in our adaptive analysis. In the framework of element-free Galerkin (EFG) method, Chung and Belytschko [23] proposed local and global error indicator for adaptive analysis. Furthermore, Lee and Zhou [24,25] have done intensive works for refinement scheme in EFG. Duarte and Oden [26] derived an error indicator that involves the computation of the interior residuals and the residuals for Neumann boundary conditions for hp-cloud method. Belytschko et al. [27] have used the adaptive analysis based on the residuals for the reproducing kernel particle method (RKPM). Gan et al. [28] presented the adaptive procedure of RKPM for 3D contact problems with elastic–plastic dynamic large deformation. Liu and Tu [29] have introduced an adaptive procedure for meshfree methods using an error indicator of energy error computed via different order of sampling in Gauss integration based on background cells. Rabczuk and Belytschko [30] proposed adaptive analysis for structured meshfree particle methods in 2-D and 3-D problems. Park and Kwon [31] used a least square mesh free method with Voronoi cells to refine interesting regions. Angulo et al. [32] implemented adaptive produce of meshfree finite point method for solving boundary value problems. In this paper, we propose an efficient error indicator and the associated refinement scheme within the framework of the NS-PIM for adaptive analysis of three-dimensional problems. The new error indicator is defined based on maximum differences of strain energy between the four of each tetrahedral cells connected to the node. A simple h-type refinement scheme is then implemented with an effective strategy for adding in nodes into the regions identified by the error indicator. An automatic three-dimensional mesh generator based on Delaunay technology is next coded to generate high quality tetrahedral meshes for each step in the adaptive process. Adaptive analysis is finally performed for a number of 3-D problems, including ones with stress concentration and singularities. The results demonstrate that the present adaptive procedure performed very well for the NS-PIM to obtain solutions of desired accuracy and with upper bounds to the exact solution. The paper is organized as follows. In Section 2, a brief description and the basic equations of NS-PIM are given. Section 3 describes the effective adaptive procedure, including the definition of error indicator, the calculation of local critical value, the strategy of refinement and the three-dimensional mesh automatic generation. In Section 4, analyses of some 3-D numerical problems are carried out to assess the capabilities of the proposed adaptive procedure. Conclusions are stated in Section 5.
2.1. Basic equations Consider a solid mechanics problem in three-dimensional domain
O bounded by G(G ¼ Gt þ Gu). The standard strong form governing equations can be expressed by the following equations [33]: Equilibrium equation in O,
where L is a differential operator in the following form, 2@ @ @ 3 0 0 @y 0 @z @x 6 7 @ @ @ 0 @y 0 @x 0 7, LT ¼ 6 @z 4 5 @ @ @ 0 0 @z 0 @y @x
o
is the vector containing stress component, u o¼ u v w is the displacement vector, n T and b ¼ bx by bz is the external body force vector. Essential boundary conditions: u ¼ up on Gu ,
ð1Þ
ð2Þ
ð3Þ
where up is the prescribed displacement on the essential boundaries. Natural boundary conditions:
r U n ¼ tp on Gt ,
ð4Þ
where tp is the prescribed traction on the natural boundaries, and n is the vector of unit outward normal on Gt. 2.2. Construction of PIM shape functions PIM shape functions are constructed using a set of small number of nodes located in a local support domain. There are two types of PIM shape functions, which have been developed with different basis functions, i.e. polynomial basis functions [7,9,12,34] and radial basis functions [11,35–38]. In this work we use both the simplest linear polynomial basis functions and effective radial basis functions to construct PIM shape functions. For the polynomial PIM, the formulations start with the following assumption: uðxÞ ¼
n X
Pi ðxÞai ¼ PT ðxÞa,
ð5Þ
i¼1
where u(x) is a field variable function defined in the Cartesian coordinate space xT ¼ x y z , Pi(x) is the basis function of monomials, which is usually built utilizing Pascal’s triangles, ai is the corresponding coefficient, and n is the number of nodes in the local support domain. For the radial PIM, using radial basis functions augmented with polynomials, a field variable function u(x) can be approximated as follows: uðxÞ ¼
n X
Ri ðxÞai þ
i¼1
m X
Pj ðxÞbj ¼ RT ðxÞa þ PT ðxÞb,
ð6Þ
j¼1
where Ri(x) and Pj(x) are radial basis functions and polynomial basis functions, respectively, ai and bj are corresponding coefficient, n is the number of nodes in the local support domain used for constructing RPIM shape functions and m is the number of polynomial terms augmented to the radial basis functions [6]. The coefficients in Eqs. (5) and (6) can be determined by enforcing the field function to be satisfied at the n nodes within the local support domain. Finally, the PIM shape functions can be obtained and a displacement component can be interpolated as uðxÞ ¼ UT ðxÞUs ,
2. Briefing on the NS-PIM
LT r þb ¼ 0
n
rT ¼ sxx syy szz sxy syz szx T
ð7Þ
where U(x) is the row-matrix of PIM shape functions and Us is the vector of nodal displacement [6]. When linear polynomial PIM shape functions are used, four vertexes of the background four-node tetrahedron cell are taken to perform the interpolation of the interest points located inside the cell, which is same as the linear FEM does. This can be easily implemented and can always ensure the invertibility of the moment matrix, as long as the four vertexes of the tetrahedron are not on a plane. When radial PIM shape functions are used, we use eight nodes to perform the approximation: four vertexes of the cell plus four nodes located at the remote vertices of the four neighboring cells. For a boundary cell with three neighboring cells, four vertexes of the cell, three vertexes of its neighboring cells and another closet to the point of interest will be selected for constructing the RPIM shape functions.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
2.3. Node-based strain smoothing operation In the formulation of NS-PIM, the compatible strains are replaced by the generalized smoothed strains which are constructed over each node-based smoothing domain. Based on the four-node tetrahedron background cells, the problem domain O is partitioned into small smoothing domains O ¼ O1[O2[y[ON with corresponding nodes, where N is the total number of field nodes. The smoothing domain for each field node is centered by the node and constructed based on the background cells of fournode tetrahedrons. As shown in Fig. 1, a typical 3D node-based smoothing domain is a polyhedron formed by connecting several hexahedral sub-smoothing-domains which are one quarter of the surrounding cells of the node. To make the illustration easier to follow, Fig. 2 shows the construction of one hexahedral subsmoothing-domain, which is constructed by connecting the midedge-points, the centroids of the surface triangles and the centroid of a tetrahedron cell. The sub-domain of the smoothing domain for node k located in the particular cell j can be obtained by connecting the mid-edge-points, the centroids of the surface triangles and the centroid of cell j. Finding out other sub-domains located in cells, which contain node k and the smoothing domain for node k can be constructed uniting all the sub-domains.
k6
Because the PIM shape functions can be discontinuous, we introduce the ‘‘smoothed’’ boundary flux-approximation strain as follow [8,39]: Z 1 _ e ij ðxk Þ ¼ ðu n þ uj ni ÞdG, ð8Þ 2Vk Gk i j where Gk is the boundary surface of the smoothing domain for node k, i.e. Ok. Substituting Eq. (7) into Eq. (8), we can obtain the smoothed strain as follows: X_ _ e ðxk Þ ¼ ð9Þ B i ðxk ÞUi , i A Gk
where Gk contains a number of nodes whose shape functions support cover node k. In three dimensional space, n_ _ _ _ _ o _ _ e ¼ e xx e yy e zz 2e xy 2e yz 2e zx n o ð10Þ UTi ¼ uxi uyi uzi , 2_ 3 b ix ðxk Þ 0 0 6 7 _ 6 0 0 7 b iy ðxk Þ 6 7 6 7 _ 6 0 0 b iz ðxk Þ 7 _ 6 7 B i ðxk Þ ¼ 6_ 7, _ 6 b iy ðxk Þ b ix ðxk Þ 7 0 6 7 6 7 _ _ 6 0 b iz ðxk Þ b iy ðxk Þ 7 4 5 _ _ b iz ðxk Þ 0 b ix ðxk Þ
k4
_ 1 b il ¼ Vk Field node k1
Mid−edge−point k k3
Centroid of the surface triangle Centroid of the tetrahedron
k2 k5 Fig. 1. A typical 3D node-based smoothing domain for node k, which is a polyhedron formed by connecting several hexahedral sub-smoothing-domains belonging to the surrounding cells of node k.
Cell j
Field node Mid-edge-point Centroid of a surface triangle
Node k Centroid of the tetrahedron cell j
Fig. 2. A background four-node tetrahedral cell j and one of the hexahedral subsmoothing-domain for node k in cell j. The sub-smoothing-domain constructed by connecting the mid-edge-points, the centroids of the surface triangles and the centroid of the tetrahedron.
1125
Z Gk
fi ðxÞnl ðxÞdO ðl ¼ x,y,zÞ:
ð11Þ
ð12Þ
where ji is PIM shape function. Applying Gauss integration among each part of smoothing domain surface Gk, the above equation can be written in algebraic form as " N # Ns g X _ 1 X b il ¼ wn ðfi ðxmn Þnl ðxm ÞÞ , ð13Þ Vk m ¼ 1 n ¼ 1 where Ns is the number of surface areas (quadrangles) of smoothing domain Ok, as shown in Fig. 1, the Ns ¼12; Ng is the number of Gauss points distributed in each area, and wn is the corresponding weight number of Gauss integration scheme. In the present method, Ng ¼4 is adopted which means that 2 2 Gauss points are used for integration on each quadrangular surface area of the smoothing domain. 2.4. Discretized system equations In the NS-PIM, it has been proved that the generalized smoothed Galerkin weak form that allows the use of discontinuous displacement functions can be used [8,39]. Z Z Z dð_ e ðuÞÞT Dð_ e ðuÞÞdO duT bdO duT tp dG ¼ 0, ð14Þ O
O
Gt
where d means the virtual change of displacement and D is the matrix of material constants. Substituting the assumed displacements (Eq. (7)) and the smoothed strains (Eq. (9)) into Eq. (14), the discretized system equation can be obtained in the following matrix form: _ Ku ¼ f, ð15Þ where f is the force vector that can be obtained as Z Z fi ¼ fi tp dG þ fi bdO, Gt
ð16Þ
O
_ and the stiffness matrix K is assembled from the sub-stiffness matrix for all the integration cells, which are exactly the node-based
1126
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
smoothing domains for the present method. N N X X _ ðkÞ _T _ _ Vk B i DB j , K ij ¼ K ij ¼ k¼1
ð17Þ
k¼1
_ where the smoothed strain matrix B can be obtained using Eq. (11).
3. Adaptive scheme In this section, a novel error indicator and refinement strategy are introduced and implemented based on the NS-PIM. The proposed adaptive analysis procedure is illustrated by the flow chart shown in Fig. 3. Adaptive analysis is an iterative process that terminate when the error in all cells falls below the specified critical value or when the number of step meets the maximum number. The proposed adaptive procedure involves the following four main points: 1. Definition of a reliable error indicator for adaptation. 2. Achievement of a local critical value to decide whether a cell should be subdivided. 3. A h-type refinement scheme for adding in nodes into the region where is identified by the error indicator. 4. Generation high quality tetrahedral meshes for three-dimensional domain based on Delaunay technology.
Therefore, the stresses in these four parts of any cell generally are different and a big difference will occur in these cells located in the volume of steep gradient stress. Making use of this solution feature of the NS-PIM, the proposed error indicator uses the maximum difference of the strain energy values among the four parts (nodes) in each tetrahedral cell as the basic error measure. For each four-node tetrahedral cell, we can compute the differences of stress components between two different nodes, namely node 1 and node 2 by the following equation: ð1Þ ð2Þ Wsð12Þ xx ¼ sxx sxx , ð12Þ ð1Þ Wsyy ¼ syy sð2Þ yy , ð1Þ ð2Þ Wsð12Þ zz ¼ szz szz , ð1Þ ð2Þ Wsð12Þ xy ¼ sxy sxy , ð1Þ ð2Þ Wsð12Þ yz ¼ syz syz , ð1Þ ð2Þ Wsð12Þ xz ¼ sxz sxz :
ð18Þ
In a similar way to Eq. (18), we can obtain the difference of strain components between node 1 and node 2 of the corresponding cell. ð1Þ ð2Þ Weð12Þ xx ¼ exx exx , ð12Þ ð1Þ Weyy ¼ eyy eð2Þ yy , ð1Þ ð2Þ Weð12Þ zz ¼ ezz ezz , ð1Þ ð2Þ Weð12Þ xy ¼ exy exy , ð1Þ ð2Þ Weð12Þ yz ¼ eyz eyz ,
3.1. Error indicator
ð1Þ ð2Þ Weð12Þ xz ¼ exz exz :
A reliable error indicator is required in numerical simulations to identify the regions for refinement, and it plays a very crucial role for adaptive analysis to effectively improve solution accuracy. In the NS-PIM, the stress (or strain) obtained is constant in each of the node-based smoothing domains. Any cell in a NS-PIM model contributes four parts, respectively, to four nodes of the cell (see Fig. 2). START
Invoke NS-PIM solver
Obtain the maximum difference in strain energy between nodes for each background cell
YES
If the adaptive steps exceed the maximum steps? for each cell
ð19Þ
We then define the energy for these differences in stresses and strains for each cell as follows: E¼
1 ðWrÞT ðWeÞ, 2
ð20Þ
where Wr ¼ ðWsxx , Wsyy , Wszz , Wsxy , Wsyz , Wsxz ÞT and We ¼ ðWexx , Weyy , Wezz , Wexy , Weyz , Wexz ÞT , substituting Eqs. (18) and (19) into Eq. (20), we can obtain for a pair of two nodes in a cell. 1 ð12Þ ð12Þ ð12Þ ð12Þ ð12Þ Wsð12Þ E12 ¼ xx Wexx þ Wsyy Weyy þ Wsxy Wexy 2 ð12Þ ð12Þ ð12Þ ð12Þ ð12Þ : ð21Þ þ Wsð12Þ zz Wezz þ Wsyz Weyz þ Wsxz Wexz In the same way, we can calculate the E13, E14, E23, E24 and E34 for the cell. Finally, we use the maximum energy for the difference in strains among four nodes as the error indicator. Thus the error indicator Ei for the ith cell can be obtained as follow :Ei : ¼ MaxðE12 ,E13 ,E14 ,E23 ,E24 ,E34 Þ:
NO
ð22Þ
Estimate the error for a cell and calculate the local critical value
2 If error indicator for the
NO
6
cell exceeds Eth ?
7 YES
5
Performing the h-type refinement scheme by adding in new nodes
4 9
Regenerate the new mesh by automatic three-dimensional mesh generator based on Delaunay technique
END Fig. 3. The flow chart of the adaptive procedure for NS-PIM.
8 10
1 Initial nodes
Added new nodes
Fig. 4. Illustration of the h-type refinement strategy.
3
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
that the cell with the sequence number of (NenZ) is chosen and its error indicator value is used as the local critical value. It is clear that Z is the ratio of the cells to be refined to all the cells at a particular step. To control the adaptation and obtain results with desired accuracy, users can change the value of Z. For each adaptive step, a bigger value of Z leads to adding in more nodes; while a smaller Z leads to adding in fewer nodes.
3.2. Local refinement criteria and refinement rate Generally the refinement in each adaptive step is always performed for a small part of the background cells. In our adaptive scheme, we simply use a refinement rate to select the cells for refinement. First, we obtain the error indicator values using Eq. (22) for each background cell, and sort out in descending order according to these values. Then the local critical value for a particular adaptive step (Eth) is determined as Eth ¼ EbNe nZc ,
1127
3.3. Refinement strategy
ð23Þ
In each of adaptive steps, we first compute the error indicator value for each background cell using Eq. (22), and then sort all these values in descending order. A cell (i) will be refined by
where Ne is the number of background cells for this step and ZA[0,1] is a refinement rate prescribed by users. Eq. (23) shows
6
2 6
6
4
7 5
9
7 8
5
7
9
5 8
9
8 10
10
3
10
1 6
6 7
5
5
9 8
9
5
7 5
9
9
8 10 10
Fig. 5. Formation of new cells splitting of a tetrahedron into four tetrahedrons and an octahedron, then splitting the octahedron into four tetrahedrons.
z
H a
a
Face z = 3,
I
Ty = 1.0
G K b
A
E a
Face x = 3, D
Tx = −1.0
F y
B
c
x C
Face z = 0, u = v = w = 0.0
Fig. 6. Problem 1: a biaxial bending of a column footing.
1128
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
adding in nodes if the corresponding value of error indicator (Ei) meets the following equation:
40
Ei 4 Eth ,
ð24Þ
36
where Eth is the local critical value. In three-dimensions, several refinement techniques have been developed in recent years for tetrahedral meshes by means of bisection. Our simple and effective h-type refinement strategy consists of two main steps: firstly, add a new vertex at the midpoint of each edge to the targeted tetrahedron, which should be refined, and then connect the midpoints of the edges, as depicted in Fig. 4. Secondly, as shown in Fig. 5, the targeted tetrahedron is divided into eight sub-tetrahedrons by cutting off the four corners by the midpoints of the six edges, the remaining octahedron is subdivided further into four tetrahedrons by adding a diagonal. If a cell and its neighbor cell are both refined, there will be some duplicated nodes, which are added twice in the sharing face of these two cells. After refinement, the duplicated nodes should be deleted and the boundary conditions must be reset for the new field nodes.
34
Step 1: 819 Nodes 3355 Tetrahedrons
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (NS-RPIM) Uniform Mesh (NS-RPIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
Strain Energy
38
32 30 28 26 24 22 20 18 0
4000
8000
12000
Step 2: 1458 Nodes 6748 Tetrahedrons
Step 3: 2787 Nodes 13847 Tetrahedrons
Step 2: 1466 Nodes 6797 Tetrahedrons
Step 3: 2801 Nodes 13987 Tetrahedrons
Fig. 8. Meshes of each adaptive step using NS-RPIM for the problem 1.
Step 1: 819 Nodes 3355 Tetrahedrons
Step 2: 1388 Nodes 6095 Tetrahedrons
20000
24000
Fig. 10. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 1.
Fig. 7. Meshes of each adaptive step using NS-PIM for the problem 1.
Step 1: 819 Nodes 3355 Tetrahedrons
16000
DOF
Step 3: 4332 Nodes 21606 Tetrahedrons
Fig. 9. Uniform refinement of the problem 1.
Step 4: 7405 Nodes 38227 Tetrahedrons
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
3.4. Mesh generation
40 Relative Error of Strain Energy (%)
1129
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (NS-RPIM) Uniform Mesh (NS-RPIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
30
20
10 9 8 7 6 5 103
104 DOF
105
Fig. 11. Comparison of convergence for the problem 1.
Unit uniform tension
4. Numerical examples
a
c
a
c b
Crack
Crack
Unit uniform tension Face z = 0.5, Tz = 1.0
z E K
a
Face y = 0, v=0
Face x = 0, u=0
B
a
J
A D
x
y b
Face, z = 0, x< = 0.25, w = 0.0
Fig. 12. Problem 2: crack under uniform tensile force. (a) whole model and (b) one-eighth of model.
Step 1: 562 Nodes 2384 Tetrahedrons
In this section, four numerical examples involving stress concentration and singularities have been employed to investigate the performance of the present adaptive procedure. In order to demonstrate the effectiveness of the adaptive refinement procedure, both uniform refinement and the proposed adaptive refinement were implemented in the framework of both NS-PIM and FEM. FEM analyses using the commercial software ABAQUS [47] are also carried out for these problems to provide reference solutions. The units used are based on international standard unit system unless specially mentioned.
4.1. Problem 1
F
a H
C
Crack
The present adaptive procedure now needs an efficient automatic three-dimensional mesh generator. While there are a large number of 2-D mesh generators available, 3-D mesh generators are scare. Due to the complexities associated with generating 3-D meshes, much work has been devoted to the automation of the 3-D mesh generation. In this paper an outstanding automatic 3-D mesh generator called Tetgen [40] is used. Tetgen can generate high quality tetrahedral meshes from three-dimensional domain based on Delaunay technology [41–43]. The geometry of three-dimensional computational domain is described using a set of vertices, segments and facets, which is called piecewise linear complex [44]. To generate the best possible shaped mesh for each adaptive step, Tetgen uses an incremental point insertion approach extended from a classical Delaunay refinement scheme proposed by Shewchuk [45–46]. Tetgen generates the new mesh based on the new nodes configuration, which preserves the information of old nodes in the previous step. Therefore, Tetgen is coded in our adaptive code.
Step 2: 2032 Nodes 9933 Tetrahedrons
In this first example, a problem involving biaxial bending of a column footing is considered. Fig. 6 shows the dimensions of the problem domain, the loading and boundary conditions applied. The column footing is fixed on the bottom face BCFE. Traction boundary conditions are a uniform shear Ty at the top face GHIK and a uniform pressure Tx on the face ABCD. The parameters are taken as Young’s modulus E¼1.0, Poisson’s ratio u ¼0.3, a¼1, b¼2, c ¼3, Ty ¼1.0 and Tx ¼ 1.0. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 564,215 nodes and 3,266,207 elements) with refinement around the region of stress concentration.
Step 3: 4163 Nodes 21846 Tetrahedrons
Fig. 13. Uniform refinement of the problem 2.
Step 4: 7384 Nodes 40163 Tetrahedrons
1130
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
Step 2: 1377 Nodes 6848 Tetrahedrons
Step 1: 562 Nodes 2384 Tetrahedrons
Step 3: 3578 Nodes 19193 Tetrahedrons
Fig. 14. Meshes of each adaptive step using NS-PIM for the problem 2.
Fig. 15. Nodes distribution of each adaptive step using NS-PIM for the problem 2.
Relative Error of Strain Energy (%)
0.106 0.104 0.102 Strain Energy
0.1 0.098 0.096 0.094 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
0.092 0.09 0.088
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
10 9 8 7 6 5 4 3 2
1
0.086 0
4,000
8000
12000 DOF
16000
20000
24000
103
104
105
DOF
Fig. 16. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 2.
Fig. 17. Comparison of convergence for the problem 2.
Adaptive analysis is conducted for this problem in the schemes of both FEM and NS-PIM (including both polynomial PIM and radial PIM). The entire adaptive analysis takes three steps to complete, ended up with 2941 nodes (FEM), 2787 nodes (NS-PIM) and 2801 nodes (NS-RPIM) from 819 nodes initially, with refinement rate Z ¼0.1. The adaptive meshes with NS-PIM and NS-RPIM for each step are shown in Figs. 7 and 8. The figures clearly indicate that the present error indicator effectively implements the refinement of nodes, which are added into the regions for stress concentration. For comparison, four uniform refinement
models with 819, 1388, 4332 and 7405 nodes are also studied, as presented in Fig. 9. Fig. 10 shows the comparisons on convergence of strain energy between the results of uniform and adaptive meshes. It clearly illustrates that both the NS-PIM and NS-RPIM provide upper bound solutions in energy norm, and FEM provides lower bounds to the exact strain energy. The proposed adaptive strategy for FEM, NS-PIM and NS-RPIM obtain much higher convergence rate than the uniformly refined models. The figure shows that the NSRPIM achieves better accuracy and tighter upper bound of
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
z Face z = 1,w = 0.0
Face z = 2, Tz = 100
A
D
B
b
J
C
Face y = 2, Ty = 100
E
Face x = 2, Tx = 100
G
1131
solutions for both uniform and adaptive models than the NS-PIM does, which is consistent with previous studies [11]. Plotting the relative errors of strain energy with respect to the degree of freedom (DOF) in Fig. 11, we can find that the proposed adaptive strategy converges much faster than the uniform model for FEM, NS-PIM and NS-RPIM. Based on the previous study [6], more support nodes are required to create RPIM shape functions and will lead to the increase of bandwidth of system matrix. Considering the time consuming problem and parameters involved in the construction of radial PIM shape functions, we use the simplest linear NS-PIM to study following cases.
a
4.2. Problem 2
F
H I M
K
b
x Face y = 1, v = 0.0
a L
The second problem is a 3-D solid with a crack under uniform tensile force. The dimensions of the problem domain, the loading and boundary conditions applied are presented in the Fig. 12a. Due to symmetry, only one-eighth of model is considered, together with the appropriate boundary and loading conditions are shown in Fig. 12b. Essential boundary conditions are imposed on the faces ABCD, ABEF, BEKJ: on the face ABCD the displacement along the z direction is set
y
Face x = 1, u = 0.0
Fig. 18. Problem 3: a rigid cube without a 1/8 part.
Step 1: 434 Nodes 1662 Tetrahedrons
Step 2: 2011 Nodes 9612 Tetrahedrons
Step 3: 5189 Nodes 27055 Tetrahedrons
Step 4: 7551 Nodes 40485 Tetrahedrons
Fig. 19. Uniform refinement of the problem 3.
Step 1: 434 Nodes 1662 Tetrahedrons
Step 2: 655 Nodes 2674 Tetrahedrons
Step 4: 1683 Nodes 7394 Tetrahedrons
Step 5: 2743 Nodes 12293 Tetrahedrons
Step 3: 1033 Nodes 4422 Tetrahedrons
Step 6: 4503 Nodes 20483 Tetrahedrons
Fig. 20. Meshes of each adaptive step using NS-PIM for the problem 3.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
to zero, on the face ABEF the displacement along the x direction is set to zero and on the face BEKJ the displacement along the y direction is set to zero. Traction boundary condition is uniform tensile force Tz imposed on the face EFHK. The parameters are taken as Young’s modulus E¼1.0, Poisson’s ratio u¼0.3, a¼ 0.5, b¼0.25, c¼1 and Tz ¼1.0. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 153,915 nodes and 852,457 elements) with refinement around the region of crack. This problem is studied by using linear NS-PIM and FEM with both uniform and adaptive meshes. Four uniform refinement models with 562, 2032, 4163 and 7384 nodes are presented in Fig. 13. The adaptive analysis started from uniform model of 562 nodes and is performed for three steps with refinement rate Z ¼0.05. The adaptive mesh and nodal distribution obtained using NS-PIM for each step are plotted respectively in Figs. 14 and 15. These figures illustrate that our error indicator is able to accurately identify the steep gradient of stresses and successfully leads to the refinement of nodes concentrated along the singular line CD, where the crack tip is located. Fig. 16 shows the comparisons on convergence of strain energy between the results of uniform meshes and adaptive meshes. The results show that by using NS-PIM and FEM with adaptive analysis procedure, a much tighter upper and lower bound to the reference energy can be obtained efficiently. It can be found that the NS-PIM has better accuracy than FEM in both uniform and adaptive meshes. As shown in Fig. 17 the relative errors of strain energy are plotted against the DOF for both uniform and adaptive refinements. One can also see that the adaptive models obtain much faster convergence than the uniform models, no matter NS-PIM or FEM.
mesh of 257,404 nodes and 1,318,451 elements) with refinement around the region of stress concentration. An adaptive analysis of this problem is conducted using both NS-PIM and FEM. Four uniform refinement models with 434, 2011, 5189 and 7551 nodes are studied, as presented in Fig. 19. The adaptive analysis started from uniform model of 434 nodes and is performed for six steps with refinement rate Z ¼ 0.05. During the adaptation, problem domain is refined with newly added nodes based on our adaptive strategy, and the adaptive mesh obtained using NS-PIM for each step is plotted in Fig. 20. The figures show that our error indicator can accurately capture the steep gradient of stresses and lead to the refinement of nodes concentrated along the three singular lines GD, GH and GF.
102 Relative Error of Strain Energy (%)
1132
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
101
103
4.3. Problem 3 The third problem is a rigid cube with 1/8 part removed and subjected to uniform tensile forces, as shown in Fig. 18. Essential boundary conditions are imposed on the faces DEFG, CDGH and HIFG: on the face DEFG the displacement along the x direction is set to zero, on the face CDGH the displacement along the y direction is set to zero and on the face HIFG the displacement along the z direction is set to zero. Traction boundary conditions are uniform forces Tx,Ty,Tz imposed on the faces BCHILM, EFILKJ and ABCDEJ, respectively. The parameters are taken as Young’s modulus E¼ 3.0 107, Poisson’s ratio u ¼0.3, a¼2, b¼ 1, Tx ¼100, Ty ¼100 and Tz ¼100. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral
105
104 DOF Fig. 22. Comparison of convergence for the problem 3.
z
Unit uniform tension y
K
G
b E
F x
a b a D
x 10-3 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
6
Strain Energy
5.5
C
A
Unit uniform tension
B b
5 4.5 J
z
K
4
Face z = 5, Tz = 1
3.5 Face x = 0, u=0
3
I
W
2.5
F
0
4000
8000
12000
16000
20000
24000
DOF Fig. 21. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 3.
Face y = 0, v = 0
y
G
D P
E A
C Face z = 0, w = 0
H R B
x
Fig. 23. Problem 4: a cube with a cube hole in the middle. (a) Whole model and (b) one-eighth of model.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
Fig. 21 shows the comparisons on convergence of strain energy between the results of uniform meshes and adaptive meshes. It can be seen that using all the models, no matter using uniform meshes or adaptive meshes, NS-PIM gives upper bound solutions (in strain energy norm) to the reference one; while FEM gives lower bound. This founding is consistent with our previous study [11,12]. The computed upper and lower bound strain energy obtained using the present adaptive strategy converges much faster than that computed using uniform mesh. In Fig. 22, it is also found that the relative error of strain energy of the adaptive model converges much faster than that of the uniform refinement for both NS-PIM and FEM. 4.4. Problem 4 The last example is the cube domain with a cube cavity in the middle, as illustrated in Fig. 23a. Traction boundary conditions are two uniform tensions Tz on the top face EFGK and on the bottom face ABCD. Due to symmetry, only one-eighth of model is considered, together with the appropriate boundary and loading conditions are
Step 1: 529 Nodes 2131 Tetrahedrons
Step 2: 2156 Nodes 10318 Tetrahedrons
shown in Fig. 23b. The parameters are taken as Young’s modulus E¼3.0 107, Poisson’s ratio u¼0.3, a¼1, b¼10, Tz ¼1. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 425,931 nodes and 2,430,905 elements) with refinement around the region of cube hole. Adaptive analyses using the NS-PIM and FEM started from the mesh of 529 nodes, and five adaptive steps are performed for both these two methods with refinement rate Z ¼0.05. For comparison, four uniform meshes with 529, 2156, 4165 and 7294 nodes are also computed, as presented in Fig. 24. The adaptive meshes for each step using NS-PIM are plotted in Fig. 25. The figure illustrates that the proposed error indicator is able to accurately capture the steep gradient of stresses and successfully leads to the refinement concentrated along the singular area, i.e. the cube hole. Fig. 26 shows that NS-PIM and FEM provide upper and lower bounds to the exact strain energy, respectively, and the adaptive models provide much tighter bounds than the uniform models. In Fig. 27, one can see that the adaptive models obtain much faster convergence rate than the uniform models for both NS-PIM and FEM.
Step 3: 4165 Nodes 21727 Tetrahedrons
Step 4: 7294 Nodes 39609 Tetrahedrons
Fig. 24. Uniform refinement of the problem 4.
Step 3: 1055 Nodes 4683 Tetrahedrons
1133
Step 1: 529 Nodes
Step 2: 738 Nodes
2131 Tetrahedrons
3110 Tetrahedrons
Step 4: 1586 Nodes 7432 Tetrahedrons
Step 5: 2484 Nodes 11945 Tetrahedrons
Fig. 25. Meshes of each adaptive step using NS-PIM for the problem 4.
1134
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
NS-RPIM achieves better accuracy and tighter upper bound for both uniform and adaptive models than the NS-PIM. (6) Compared with uniform refinements, the proposed adaptive strategy based on NS-PIM improves convergence significantly. It clearly indicates that the presented adaptive procedure can achieve the desired accuracy with far less nodes.
x 10-6 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
Strain Energy
2.0885
2.088
Incorporated with an ideal error indicator and effective refinement scheme, the present adaptive procedure based on NS-PIM for 3-D problems has achieved very good performance. All the results have demonstrated that the proposed adaptive procedure is simple, robust and efficient.
2.0875
2.087
0
4000
8000
12000 DOF
16000
20000
24000
Fig. 26. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 4.
Relative Error of Strain Energy (%)
10-1 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
104 DOF
The author wish to thank the support of the China-funded Postgraduates’ Studying Aboard Program for Building Top University and the National Natural Science Foundation of China. The authors also give sincerely thanks to the support of Centre for ACES at National University of Singapore. It is also supported by the program for Changjiang Scholar and Innovative Research Team in University (5311050050037). References
10-2
103
Acknowledgments
105
Fig. 27. Comparison of convergence for the problem 4.
5. Conclusions In this paper, an effective adaptive procedure using the NS-PIM is proposed for 3-D problems. Elasticity problems involving stress concentration and solution singularities have been studied to examine the performance of the proposed adaptive procedure. From the results obtained in the numerical experiments, the following conclusions can be drawn: (1) NS-PIM works well with the simplest four-node tetrahedral elements, which is easy to be implemented for complicated shape. (2) The proposed error indicator works well for all the four numerical cases studied in this work. It effectively captures the location of steep gradient of stress and carried out the refinement in the area of stress concentration. (3) The h-type refinement scheme is performed by adding in new nodes to refine old mesh and the information of old nodes in previous step can be saved. (4) Based on new nodes configuration, an automatic 3-D mesh generator is used to produce high quality tetrahedral meshes. (5) The NS-PIM models, including the NS-PIM with polynomial PIM shape functions and NS-RPIM with radial PIM shape functions, can provide upper bound solutions in energy norm.
[1] Zienkiweicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering 1987;142:337–57. [2] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique. International Journal for Numerical Methods in Engineering 1992;33:1331–64. [3] Lee CK, Lo SH. Automatic adaptive 3D finite element refinement using different-order tetrahedral element. International Journal for Numerical Methods in Engineering 1997;40:2195–226. [4] Lee CK, Lo SH. A full 3D finite element analysis using adaptive refinement and PCG solver with back interpolation. International Journal for Numerical Methods in Engineering 1999;170:39–64. [5] Liu GR, Liu MB. Smoothed particle hydrodynamics-A meshfree practical method. World Scientific; 2003. [6] Liu GR. Meshfree methods: moving beyond the finite element method. 2nd ed. CRC Press; 2009. [7] Liu GR, Zhang GY, Dai KY, Wang YY, Huang HT, Zhong ZH, et al. A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems. International Journal of Computational Methods 2005;2:645–65. [8] Liu GR. A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods. International Journal of Computational Methods 2008;5:199–236. [9] Zhang GY, Liu GR, Wang YY, Huang HT, Zhong ZH, Li GY, et al. A linearly conforming point interpolation method (LC-PIM) for three-dimensional elasticity problems. International Journal for Numerical Methods in Engineering 2007;72:1524–43. [10] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering 2001;50:435–66. [11] Zhang GY, Liu GR, Nguyen TT, Song CX, Han X, Zhong ZH, et al. The upper bound property for solid mechanics of the linearly conforming radial point interpolation method (LC-RPIM). International Journal of Computational Methods 2007;4:521–41. [12] Liu GR, Zhang GY. Upper bound solution to elasticity problems: a unique property of the linearly conforming point interpolation method (LC-PIM). International Journal for Numerical Methods in Engineering 2007;74:1128–61. [13] Zhang GY, Liu GR. An efficient adaptive analysis procedure for certified solutions with exact bounds of strain energy for elasticity problems. Finite Elements in Analysis and Design 2008;44:831–41. [14] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 2: error estimates and adaptivity. International Journal for Numerical Methods in Engineering 1992;33:1365–82. [15] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery (SPR) and adaptive finite element refinement. Computer Methods in Applied Mechanics and Engineering 1992;101:207–24. [16] Babuska I, Rheinboldt WC. A-posteriori error estimates for the finite element method. International Journal for Numerical Methods in Engineering 1978;12:1597–615. [17] Babuska I, Rheinboldt WC. Error estimates for adaptive finite element computations. SIAM Journal on Numerical Analysis 1989;15:739–54.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
[18] Oden JT, Demkowicz L, Rachowicz W, Westermann TA. Toward a universal h-p adaptive finite element strategy. Part 2: a posteriori error estimation. Computer Methods in Applied Mechanics and Engineering 1989;77:113–80. [19] Ainsworth M, Oden JT. A unified approach to a posteriori error estimation using element residual methods. Numerische Mathematik 1993;65:23–50. [20] Bernard BT, Liu GR, Zhang GY, Lu C. A residual based error estimator using radial basis functions. Finite Elements in Analysis and Design 2008;44: 631–45. [21] Babuska I, Suri M. The p- and h-p version of the finite element method, an overview. Computer Methods in Applied Mechanics and Engineering 1990;80:5–20. [22] Zienkiewicz OC, Taylor RL, editors. Oxford: Butterworth-Heinemann; 2000. [V1: The basis]. [23] Chung HJ, Belytschko T. An error estimate in the EFG method. Computational Mechanics 1998;21:91–100. [24] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: part I: stress recovery and a posteriori error estimation. Computers & Structures 2004;82:413–28. [25] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: part II: adaptive refinement. Computers & Structures 2004;82:429–43. [26] Durate CA, Oden JT. An h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 1996;139:237–62. [27] Belytschko T, Liu WK, Singer M. On adaptivity and error criteria for meshfree methods. Advances in adaptive computational methods in mechanics. Elsevier Science 1998;47:217–28. [28] Gan NF, Li GY, Long SY. 3D adaptive RKPM method for contact problems with elastic–plastic dynamic large deformation. Engineering Analysis with Boundary Elements 2009;33:1211–22. [29] Liu GR, Tu ZH. An adaptive procedure based on background cells for meshless methods. Computer Methods in Applied Mechanics and Engineering 2002;191:1923–43. [30] Rabczuk T, Belytschko T. Adaptive for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering 2005;63:1559–82. [31] Park SH, Kwon KC. A posteriori error estimates and an adaptive scheme of least-square meshfree method. International Journal for Numerical Methods in Engineering 2003;58:1213–50. [32] Angulo A, Pozo L, Perazzo F. A posteriori error estimator and an adaptive technique in meshless finite points method. Engineering Analysis with Boundary Elements 2009;33:1322–38.
1135
[33] Timoshenko SP, Goodier JN. Theory of elasticity.3rd ed.. NewYork: McGrawhill; 1970. [34] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering 2001;50:937–51. [35] Li Y, Liu GR, Luan MT, Dai KY, Zhong ZH, Li GY, et al. Contact analysis for solids based on linearly conforming radial point interpolation method. Computational Mechanics 2007;39(4):537–54. [36] Liu GR, Li Y, Dai KY, Luan MT, Xue W. A linearly conforming radial point interpolation method for solid mechanics problems. International Journal of Computational Methods 2006;3:401–28. [37] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. International Journal for Numerical Methods in Engineering 2002;54:1623–48. [38] Liu GR, Zhang GY. Edge-based smoothed point interpolation methods. International Journal of Computational Methods 2008;5(4):621–46. [39] Liu GR. A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: part I theory. International Journal for Numerical Methods in Engineering 2010;81. 1093-26. [40] Hang S. Adaptive tetrahedral mesh generation by constrained Delaunay refinement. International Journal for Numerical Methods in Engineering 2008;75:856–80. [41] Frey PJ, George Pl. Mesh generation: application to finite elements. Oxford, U.k: Hermes Science; 2000. [42] Weatherill NP, Hassan O. Efficient three-dimentional Delaunay triangulation with automatic point creation and imposed boundary constraints. International Journal for Numerical Methods in Engineering 1994;37:2005–39. [43] Chew PL. Guaranteed-quality Delaunay meshing in 3D. In: Proceedings of the 13th annual ACM symposium on computational geometry. Nice, France;1997. p. 391–93. [44] Miller GL, Talmor D, Teng SH, Walkington N, Wang H. Control volume meshes using sphere packing: generation, refinement and coarsening. In: Proceeding of 5th international meshing roundtable; 1996. p. 47–61. [45] Shewchuk JR. Delaunay refinement mesh generation. PhD thesis. Department of Computer Science, Carnegie Mellon University, Pittsburgh, PA. Available as Technical report CMU-CS-97-137; 1997. [46] Shewchuk JR. Triangle: engineering a 2d quality mesh generator and delaunay triangulator. Applied Computational Geometry: Towards Geometric Engineering 1996;1148:203–22. [47] ABAQUS Theory Manual and User Manual, version 6.4. Hibbit, Karlsson & Sorensen Inc., USA; 2004.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A three-dimensional adaptive analysis using the meshfree node-based smoothed point interpolation method (NS-PIM) Q. Tang a,b, G.Y. Zhang c,n, G.R. Liu d, Z.H. Zhong a, Z.C. He a,b a
State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, China Centre for Advanced Computations in Engineering Science (ACES), Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117576, Singapore c School of Mechanical and Chemical Engineering, The University of Western Australia. 35 Stirling Highway, WA 6009, Australia d School of Aerospace Systems, University of Cincinnati, Cincinnati, OH 45221-0070, USA b
a r t i c l e i n f o
abstract
Article history: Received 20 March 2009 Accepted 9 May 2011 Available online 7 June 2011
In this paper, a three-dimensional (3-D) adaptive analysis procedure is proposed using the meshfree node-based smoothed point interpolation method (NS-PIM). Previous study has shown that the NS-PIM works well with the simplest four-node tetrahedral mesh, which is easy to be implemented for complicated geometry. In contrast to the displacement-based FEM providing lower bound solutions, the NS-PIM possesses the attractive property of providing upper bound solutions in strain energy norm. In the present adaptive procedure, a novel error indicator is devised for NS-PIM settings, which evaluates the maximum difference of strain energy values among four nodes in each of the tetrahedral cells. A simple h-type local refinement scheme is adopted and coupled with 3-D mesh automatic generator based on Delaunay technology. Numerical results indicate that the adaptive refinement procedure can effectively capture the stress concentration and solution singularities, and carry out local refinement automatically. The present adaptive procedure achieves much higher convergence in strain energy solution compared to the uniform refinement, and obtains upper bound solution in strain energy efficiently for force driven problems. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Numerical methods Meshfree methods Error indicator Adaptive analysis Point interpolation Tetrahedral mesh generation
1. Introduction It is well known that adaptive analysis is important in the computational mechanics and engineering structural designs. The ultimate objective of adaptive analysis is to achieve desired high accuracy with minimum computational cost. Adaptive mesh refinement techniques along with proper error analysis have been well studied in the traditional FEM [1–4], but much less welldeveloped for meshfree method [5,6]. The meshfree node-based point interpolation method (NS-PIM or LC-PIM) [7] has been recently developed using the generalized smoothed Galerkin (GS- Galerkin) weak form [8] and PIM shape functions with a set of small number of nodes located in a local support domain. The PIM (polynomial PIM or radial PIM) shape functions possess Delta function property, which allows straightforward imposition of point essential boundary conditions. The use of node-based smoothing domains [7,9–11] provides the so-called ‘‘soft’’ effect to the discrete model and gives a number of good features. Compared with the linear FEM, NS-PIM can
n
Corresponding author. E-mail address: [email protected] (G.Y. Zhang).
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.05.019
obtain better accuracy and higher convergence rate, especially for stress results. Recently, Liu and Zhang [12] made a thorough theoretical study on the NS-PIM and obtained some important conclusions including the upper bound property. NS-PIM is easy to be implemented and suitable for adaptive analysis. An adaptive analysis using the NS-PIM has been conducted to certify solutions with exact bounds of strain energy for 2-D problems [13]. A reliable and efficient error indicator and the associated refinement strategy are crucial issues for the adaptive procedure, especially for 3-D problems. In general, two distinct types of procedures are currently available for deriving error indicators: the recovery based error indicator and the residual based error indicator. The recovery based error indicator was first introduced by Zienkiewicz and Zhu [1] in 1987, by constructing locally an improved solution from the approximation. The recovery processes play a critical role in the computation of this error indicator and there are many papers published on this recovery methodology [2,14,15]. Residual based error indicators make use of the residual of the numerical approximation, either explicitly or implicitly, which offers a very effective alternative [16–20]. Refinement schemes can also be classified into three categories: h-refinement, p-refinement and r-refinement. The h-refinement scheme changes the size of element in a localized fashion based on
1124
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
the error indicator. The p-refinement scheme is to increase the order of the polynomial, and the r-refinement keeps the total number of nodes unchanged but to adjust their position to obtain an optimal approximation [21,22]. Because the NS-PIM works well with lower order interpolation of displacement field, this work chooses to use the h-refinement in our adaptive analysis. In the framework of element-free Galerkin (EFG) method, Chung and Belytschko [23] proposed local and global error indicator for adaptive analysis. Furthermore, Lee and Zhou [24,25] have done intensive works for refinement scheme in EFG. Duarte and Oden [26] derived an error indicator that involves the computation of the interior residuals and the residuals for Neumann boundary conditions for hp-cloud method. Belytschko et al. [27] have used the adaptive analysis based on the residuals for the reproducing kernel particle method (RKPM). Gan et al. [28] presented the adaptive procedure of RKPM for 3D contact problems with elastic–plastic dynamic large deformation. Liu and Tu [29] have introduced an adaptive procedure for meshfree methods using an error indicator of energy error computed via different order of sampling in Gauss integration based on background cells. Rabczuk and Belytschko [30] proposed adaptive analysis for structured meshfree particle methods in 2-D and 3-D problems. Park and Kwon [31] used a least square mesh free method with Voronoi cells to refine interesting regions. Angulo et al. [32] implemented adaptive produce of meshfree finite point method for solving boundary value problems. In this paper, we propose an efficient error indicator and the associated refinement scheme within the framework of the NS-PIM for adaptive analysis of three-dimensional problems. The new error indicator is defined based on maximum differences of strain energy between the four of each tetrahedral cells connected to the node. A simple h-type refinement scheme is then implemented with an effective strategy for adding in nodes into the regions identified by the error indicator. An automatic three-dimensional mesh generator based on Delaunay technology is next coded to generate high quality tetrahedral meshes for each step in the adaptive process. Adaptive analysis is finally performed for a number of 3-D problems, including ones with stress concentration and singularities. The results demonstrate that the present adaptive procedure performed very well for the NS-PIM to obtain solutions of desired accuracy and with upper bounds to the exact solution. The paper is organized as follows. In Section 2, a brief description and the basic equations of NS-PIM are given. Section 3 describes the effective adaptive procedure, including the definition of error indicator, the calculation of local critical value, the strategy of refinement and the three-dimensional mesh automatic generation. In Section 4, analyses of some 3-D numerical problems are carried out to assess the capabilities of the proposed adaptive procedure. Conclusions are stated in Section 5.
2.1. Basic equations Consider a solid mechanics problem in three-dimensional domain
O bounded by G(G ¼ Gt þ Gu). The standard strong form governing equations can be expressed by the following equations [33]: Equilibrium equation in O,
where L is a differential operator in the following form, 2@ @ @ 3 0 0 @y 0 @z @x 6 7 @ @ @ 0 @y 0 @x 0 7, LT ¼ 6 @z 4 5 @ @ @ 0 0 @z 0 @y @x
o
is the vector containing stress component, u o¼ u v w is the displacement vector, n T and b ¼ bx by bz is the external body force vector. Essential boundary conditions: u ¼ up on Gu ,
ð1Þ
ð2Þ
ð3Þ
where up is the prescribed displacement on the essential boundaries. Natural boundary conditions:
r U n ¼ tp on Gt ,
ð4Þ
where tp is the prescribed traction on the natural boundaries, and n is the vector of unit outward normal on Gt. 2.2. Construction of PIM shape functions PIM shape functions are constructed using a set of small number of nodes located in a local support domain. There are two types of PIM shape functions, which have been developed with different basis functions, i.e. polynomial basis functions [7,9,12,34] and radial basis functions [11,35–38]. In this work we use both the simplest linear polynomial basis functions and effective radial basis functions to construct PIM shape functions. For the polynomial PIM, the formulations start with the following assumption: uðxÞ ¼
n X
Pi ðxÞai ¼ PT ðxÞa,
ð5Þ
i¼1
where u(x) is a field variable function defined in the Cartesian coordinate space xT ¼ x y z , Pi(x) is the basis function of monomials, which is usually built utilizing Pascal’s triangles, ai is the corresponding coefficient, and n is the number of nodes in the local support domain. For the radial PIM, using radial basis functions augmented with polynomials, a field variable function u(x) can be approximated as follows: uðxÞ ¼
n X
Ri ðxÞai þ
i¼1
m X
Pj ðxÞbj ¼ RT ðxÞa þ PT ðxÞb,
ð6Þ
j¼1
where Ri(x) and Pj(x) are radial basis functions and polynomial basis functions, respectively, ai and bj are corresponding coefficient, n is the number of nodes in the local support domain used for constructing RPIM shape functions and m is the number of polynomial terms augmented to the radial basis functions [6]. The coefficients in Eqs. (5) and (6) can be determined by enforcing the field function to be satisfied at the n nodes within the local support domain. Finally, the PIM shape functions can be obtained and a displacement component can be interpolated as uðxÞ ¼ UT ðxÞUs ,
2. Briefing on the NS-PIM
LT r þb ¼ 0
n
rT ¼ sxx syy szz sxy syz szx T
ð7Þ
where U(x) is the row-matrix of PIM shape functions and Us is the vector of nodal displacement [6]. When linear polynomial PIM shape functions are used, four vertexes of the background four-node tetrahedron cell are taken to perform the interpolation of the interest points located inside the cell, which is same as the linear FEM does. This can be easily implemented and can always ensure the invertibility of the moment matrix, as long as the four vertexes of the tetrahedron are not on a plane. When radial PIM shape functions are used, we use eight nodes to perform the approximation: four vertexes of the cell plus four nodes located at the remote vertices of the four neighboring cells. For a boundary cell with three neighboring cells, four vertexes of the cell, three vertexes of its neighboring cells and another closet to the point of interest will be selected for constructing the RPIM shape functions.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
2.3. Node-based strain smoothing operation In the formulation of NS-PIM, the compatible strains are replaced by the generalized smoothed strains which are constructed over each node-based smoothing domain. Based on the four-node tetrahedron background cells, the problem domain O is partitioned into small smoothing domains O ¼ O1[O2[y[ON with corresponding nodes, where N is the total number of field nodes. The smoothing domain for each field node is centered by the node and constructed based on the background cells of fournode tetrahedrons. As shown in Fig. 1, a typical 3D node-based smoothing domain is a polyhedron formed by connecting several hexahedral sub-smoothing-domains which are one quarter of the surrounding cells of the node. To make the illustration easier to follow, Fig. 2 shows the construction of one hexahedral subsmoothing-domain, which is constructed by connecting the midedge-points, the centroids of the surface triangles and the centroid of a tetrahedron cell. The sub-domain of the smoothing domain for node k located in the particular cell j can be obtained by connecting the mid-edge-points, the centroids of the surface triangles and the centroid of cell j. Finding out other sub-domains located in cells, which contain node k and the smoothing domain for node k can be constructed uniting all the sub-domains.
k6
Because the PIM shape functions can be discontinuous, we introduce the ‘‘smoothed’’ boundary flux-approximation strain as follow [8,39]: Z 1 _ e ij ðxk Þ ¼ ðu n þ uj ni ÞdG, ð8Þ 2Vk Gk i j where Gk is the boundary surface of the smoothing domain for node k, i.e. Ok. Substituting Eq. (7) into Eq. (8), we can obtain the smoothed strain as follows: X_ _ e ðxk Þ ¼ ð9Þ B i ðxk ÞUi , i A Gk
where Gk contains a number of nodes whose shape functions support cover node k. In three dimensional space, n_ _ _ _ _ o _ _ e ¼ e xx e yy e zz 2e xy 2e yz 2e zx n o ð10Þ UTi ¼ uxi uyi uzi , 2_ 3 b ix ðxk Þ 0 0 6 7 _ 6 0 0 7 b iy ðxk Þ 6 7 6 7 _ 6 0 0 b iz ðxk Þ 7 _ 6 7 B i ðxk Þ ¼ 6_ 7, _ 6 b iy ðxk Þ b ix ðxk Þ 7 0 6 7 6 7 _ _ 6 0 b iz ðxk Þ b iy ðxk Þ 7 4 5 _ _ b iz ðxk Þ 0 b ix ðxk Þ
k4
_ 1 b il ¼ Vk Field node k1
Mid−edge−point k k3
Centroid of the surface triangle Centroid of the tetrahedron
k2 k5 Fig. 1. A typical 3D node-based smoothing domain for node k, which is a polyhedron formed by connecting several hexahedral sub-smoothing-domains belonging to the surrounding cells of node k.
Cell j
Field node Mid-edge-point Centroid of a surface triangle
Node k Centroid of the tetrahedron cell j
Fig. 2. A background four-node tetrahedral cell j and one of the hexahedral subsmoothing-domain for node k in cell j. The sub-smoothing-domain constructed by connecting the mid-edge-points, the centroids of the surface triangles and the centroid of the tetrahedron.
1125
Z Gk
fi ðxÞnl ðxÞdO ðl ¼ x,y,zÞ:
ð11Þ
ð12Þ
where ji is PIM shape function. Applying Gauss integration among each part of smoothing domain surface Gk, the above equation can be written in algebraic form as " N # Ns g X _ 1 X b il ¼ wn ðfi ðxmn Þnl ðxm ÞÞ , ð13Þ Vk m ¼ 1 n ¼ 1 where Ns is the number of surface areas (quadrangles) of smoothing domain Ok, as shown in Fig. 1, the Ns ¼12; Ng is the number of Gauss points distributed in each area, and wn is the corresponding weight number of Gauss integration scheme. In the present method, Ng ¼4 is adopted which means that 2 2 Gauss points are used for integration on each quadrangular surface area of the smoothing domain. 2.4. Discretized system equations In the NS-PIM, it has been proved that the generalized smoothed Galerkin weak form that allows the use of discontinuous displacement functions can be used [8,39]. Z Z Z dð_ e ðuÞÞT Dð_ e ðuÞÞdO duT bdO duT tp dG ¼ 0, ð14Þ O
O
Gt
where d means the virtual change of displacement and D is the matrix of material constants. Substituting the assumed displacements (Eq. (7)) and the smoothed strains (Eq. (9)) into Eq. (14), the discretized system equation can be obtained in the following matrix form: _ Ku ¼ f, ð15Þ where f is the force vector that can be obtained as Z Z fi ¼ fi tp dG þ fi bdO, Gt
ð16Þ
O
_ and the stiffness matrix K is assembled from the sub-stiffness matrix for all the integration cells, which are exactly the node-based
1126
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
smoothing domains for the present method. N N X X _ ðkÞ _T _ _ Vk B i DB j , K ij ¼ K ij ¼ k¼1
ð17Þ
k¼1
_ where the smoothed strain matrix B can be obtained using Eq. (11).
3. Adaptive scheme In this section, a novel error indicator and refinement strategy are introduced and implemented based on the NS-PIM. The proposed adaptive analysis procedure is illustrated by the flow chart shown in Fig. 3. Adaptive analysis is an iterative process that terminate when the error in all cells falls below the specified critical value or when the number of step meets the maximum number. The proposed adaptive procedure involves the following four main points: 1. Definition of a reliable error indicator for adaptation. 2. Achievement of a local critical value to decide whether a cell should be subdivided. 3. A h-type refinement scheme for adding in nodes into the region where is identified by the error indicator. 4. Generation high quality tetrahedral meshes for three-dimensional domain based on Delaunay technology.
Therefore, the stresses in these four parts of any cell generally are different and a big difference will occur in these cells located in the volume of steep gradient stress. Making use of this solution feature of the NS-PIM, the proposed error indicator uses the maximum difference of the strain energy values among the four parts (nodes) in each tetrahedral cell as the basic error measure. For each four-node tetrahedral cell, we can compute the differences of stress components between two different nodes, namely node 1 and node 2 by the following equation: ð1Þ ð2Þ Wsð12Þ xx ¼ sxx sxx , ð12Þ ð1Þ Wsyy ¼ syy sð2Þ yy , ð1Þ ð2Þ Wsð12Þ zz ¼ szz szz , ð1Þ ð2Þ Wsð12Þ xy ¼ sxy sxy , ð1Þ ð2Þ Wsð12Þ yz ¼ syz syz , ð1Þ ð2Þ Wsð12Þ xz ¼ sxz sxz :
ð18Þ
In a similar way to Eq. (18), we can obtain the difference of strain components between node 1 and node 2 of the corresponding cell. ð1Þ ð2Þ Weð12Þ xx ¼ exx exx , ð12Þ ð1Þ Weyy ¼ eyy eð2Þ yy , ð1Þ ð2Þ Weð12Þ zz ¼ ezz ezz , ð1Þ ð2Þ Weð12Þ xy ¼ exy exy , ð1Þ ð2Þ Weð12Þ yz ¼ eyz eyz ,
3.1. Error indicator
ð1Þ ð2Þ Weð12Þ xz ¼ exz exz :
A reliable error indicator is required in numerical simulations to identify the regions for refinement, and it plays a very crucial role for adaptive analysis to effectively improve solution accuracy. In the NS-PIM, the stress (or strain) obtained is constant in each of the node-based smoothing domains. Any cell in a NS-PIM model contributes four parts, respectively, to four nodes of the cell (see Fig. 2). START
Invoke NS-PIM solver
Obtain the maximum difference in strain energy between nodes for each background cell
YES
If the adaptive steps exceed the maximum steps? for each cell
ð19Þ
We then define the energy for these differences in stresses and strains for each cell as follows: E¼
1 ðWrÞT ðWeÞ, 2
ð20Þ
where Wr ¼ ðWsxx , Wsyy , Wszz , Wsxy , Wsyz , Wsxz ÞT and We ¼ ðWexx , Weyy , Wezz , Wexy , Weyz , Wexz ÞT , substituting Eqs. (18) and (19) into Eq. (20), we can obtain for a pair of two nodes in a cell. 1 ð12Þ ð12Þ ð12Þ ð12Þ ð12Þ Wsð12Þ E12 ¼ xx Wexx þ Wsyy Weyy þ Wsxy Wexy 2 ð12Þ ð12Þ ð12Þ ð12Þ ð12Þ : ð21Þ þ Wsð12Þ zz Wezz þ Wsyz Weyz þ Wsxz Wexz In the same way, we can calculate the E13, E14, E23, E24 and E34 for the cell. Finally, we use the maximum energy for the difference in strains among four nodes as the error indicator. Thus the error indicator Ei for the ith cell can be obtained as follow :Ei : ¼ MaxðE12 ,E13 ,E14 ,E23 ,E24 ,E34 Þ:
NO
ð22Þ
Estimate the error for a cell and calculate the local critical value
2 If error indicator for the
NO
6
cell exceeds Eth ?
7 YES
5
Performing the h-type refinement scheme by adding in new nodes
4 9
Regenerate the new mesh by automatic three-dimensional mesh generator based on Delaunay technique
END Fig. 3. The flow chart of the adaptive procedure for NS-PIM.
8 10
1 Initial nodes
Added new nodes
Fig. 4. Illustration of the h-type refinement strategy.
3
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
that the cell with the sequence number of (NenZ) is chosen and its error indicator value is used as the local critical value. It is clear that Z is the ratio of the cells to be refined to all the cells at a particular step. To control the adaptation and obtain results with desired accuracy, users can change the value of Z. For each adaptive step, a bigger value of Z leads to adding in more nodes; while a smaller Z leads to adding in fewer nodes.
3.2. Local refinement criteria and refinement rate Generally the refinement in each adaptive step is always performed for a small part of the background cells. In our adaptive scheme, we simply use a refinement rate to select the cells for refinement. First, we obtain the error indicator values using Eq. (22) for each background cell, and sort out in descending order according to these values. Then the local critical value for a particular adaptive step (Eth) is determined as Eth ¼ EbNe nZc ,
1127
3.3. Refinement strategy
ð23Þ
In each of adaptive steps, we first compute the error indicator value for each background cell using Eq. (22), and then sort all these values in descending order. A cell (i) will be refined by
where Ne is the number of background cells for this step and ZA[0,1] is a refinement rate prescribed by users. Eq. (23) shows
6
2 6
6
4
7 5
9
7 8
5
7
9
5 8
9
8 10
10
3
10
1 6
6 7
5
5
9 8
9
5
7 5
9
9
8 10 10
Fig. 5. Formation of new cells splitting of a tetrahedron into four tetrahedrons and an octahedron, then splitting the octahedron into four tetrahedrons.
z
H a
a
Face z = 3,
I
Ty = 1.0
G K b
A
E a
Face x = 3, D
Tx = −1.0
F y
B
c
x C
Face z = 0, u = v = w = 0.0
Fig. 6. Problem 1: a biaxial bending of a column footing.
1128
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
adding in nodes if the corresponding value of error indicator (Ei) meets the following equation:
40
Ei 4 Eth ,
ð24Þ
36
where Eth is the local critical value. In three-dimensions, several refinement techniques have been developed in recent years for tetrahedral meshes by means of bisection. Our simple and effective h-type refinement strategy consists of two main steps: firstly, add a new vertex at the midpoint of each edge to the targeted tetrahedron, which should be refined, and then connect the midpoints of the edges, as depicted in Fig. 4. Secondly, as shown in Fig. 5, the targeted tetrahedron is divided into eight sub-tetrahedrons by cutting off the four corners by the midpoints of the six edges, the remaining octahedron is subdivided further into four tetrahedrons by adding a diagonal. If a cell and its neighbor cell are both refined, there will be some duplicated nodes, which are added twice in the sharing face of these two cells. After refinement, the duplicated nodes should be deleted and the boundary conditions must be reset for the new field nodes.
34
Step 1: 819 Nodes 3355 Tetrahedrons
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (NS-RPIM) Uniform Mesh (NS-RPIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
Strain Energy
38
32 30 28 26 24 22 20 18 0
4000
8000
12000
Step 2: 1458 Nodes 6748 Tetrahedrons
Step 3: 2787 Nodes 13847 Tetrahedrons
Step 2: 1466 Nodes 6797 Tetrahedrons
Step 3: 2801 Nodes 13987 Tetrahedrons
Fig. 8. Meshes of each adaptive step using NS-RPIM for the problem 1.
Step 1: 819 Nodes 3355 Tetrahedrons
Step 2: 1388 Nodes 6095 Tetrahedrons
20000
24000
Fig. 10. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 1.
Fig. 7. Meshes of each adaptive step using NS-PIM for the problem 1.
Step 1: 819 Nodes 3355 Tetrahedrons
16000
DOF
Step 3: 4332 Nodes 21606 Tetrahedrons
Fig. 9. Uniform refinement of the problem 1.
Step 4: 7405 Nodes 38227 Tetrahedrons
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
3.4. Mesh generation
40 Relative Error of Strain Energy (%)
1129
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (NS-RPIM) Uniform Mesh (NS-RPIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
30
20
10 9 8 7 6 5 103
104 DOF
105
Fig. 11. Comparison of convergence for the problem 1.
Unit uniform tension
4. Numerical examples
a
c
a
c b
Crack
Crack
Unit uniform tension Face z = 0.5, Tz = 1.0
z E K
a
Face y = 0, v=0
Face x = 0, u=0
B
a
J
A D
x
y b
Face, z = 0, x< = 0.25, w = 0.0
Fig. 12. Problem 2: crack under uniform tensile force. (a) whole model and (b) one-eighth of model.
Step 1: 562 Nodes 2384 Tetrahedrons
In this section, four numerical examples involving stress concentration and singularities have been employed to investigate the performance of the present adaptive procedure. In order to demonstrate the effectiveness of the adaptive refinement procedure, both uniform refinement and the proposed adaptive refinement were implemented in the framework of both NS-PIM and FEM. FEM analyses using the commercial software ABAQUS [47] are also carried out for these problems to provide reference solutions. The units used are based on international standard unit system unless specially mentioned.
4.1. Problem 1
F
a H
C
Crack
The present adaptive procedure now needs an efficient automatic three-dimensional mesh generator. While there are a large number of 2-D mesh generators available, 3-D mesh generators are scare. Due to the complexities associated with generating 3-D meshes, much work has been devoted to the automation of the 3-D mesh generation. In this paper an outstanding automatic 3-D mesh generator called Tetgen [40] is used. Tetgen can generate high quality tetrahedral meshes from three-dimensional domain based on Delaunay technology [41–43]. The geometry of three-dimensional computational domain is described using a set of vertices, segments and facets, which is called piecewise linear complex [44]. To generate the best possible shaped mesh for each adaptive step, Tetgen uses an incremental point insertion approach extended from a classical Delaunay refinement scheme proposed by Shewchuk [45–46]. Tetgen generates the new mesh based on the new nodes configuration, which preserves the information of old nodes in the previous step. Therefore, Tetgen is coded in our adaptive code.
Step 2: 2032 Nodes 9933 Tetrahedrons
In this first example, a problem involving biaxial bending of a column footing is considered. Fig. 6 shows the dimensions of the problem domain, the loading and boundary conditions applied. The column footing is fixed on the bottom face BCFE. Traction boundary conditions are a uniform shear Ty at the top face GHIK and a uniform pressure Tx on the face ABCD. The parameters are taken as Young’s modulus E¼1.0, Poisson’s ratio u ¼0.3, a¼1, b¼2, c ¼3, Ty ¼1.0 and Tx ¼ 1.0. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 564,215 nodes and 3,266,207 elements) with refinement around the region of stress concentration.
Step 3: 4163 Nodes 21846 Tetrahedrons
Fig. 13. Uniform refinement of the problem 2.
Step 4: 7384 Nodes 40163 Tetrahedrons
1130
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
Step 2: 1377 Nodes 6848 Tetrahedrons
Step 1: 562 Nodes 2384 Tetrahedrons
Step 3: 3578 Nodes 19193 Tetrahedrons
Fig. 14. Meshes of each adaptive step using NS-PIM for the problem 2.
Fig. 15. Nodes distribution of each adaptive step using NS-PIM for the problem 2.
Relative Error of Strain Energy (%)
0.106 0.104 0.102 Strain Energy
0.1 0.098 0.096 0.094 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
0.092 0.09 0.088
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
10 9 8 7 6 5 4 3 2
1
0.086 0
4,000
8000
12000 DOF
16000
20000
24000
103
104
105
DOF
Fig. 16. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 2.
Fig. 17. Comparison of convergence for the problem 2.
Adaptive analysis is conducted for this problem in the schemes of both FEM and NS-PIM (including both polynomial PIM and radial PIM). The entire adaptive analysis takes three steps to complete, ended up with 2941 nodes (FEM), 2787 nodes (NS-PIM) and 2801 nodes (NS-RPIM) from 819 nodes initially, with refinement rate Z ¼0.1. The adaptive meshes with NS-PIM and NS-RPIM for each step are shown in Figs. 7 and 8. The figures clearly indicate that the present error indicator effectively implements the refinement of nodes, which are added into the regions for stress concentration. For comparison, four uniform refinement
models with 819, 1388, 4332 and 7405 nodes are also studied, as presented in Fig. 9. Fig. 10 shows the comparisons on convergence of strain energy between the results of uniform and adaptive meshes. It clearly illustrates that both the NS-PIM and NS-RPIM provide upper bound solutions in energy norm, and FEM provides lower bounds to the exact strain energy. The proposed adaptive strategy for FEM, NS-PIM and NS-RPIM obtain much higher convergence rate than the uniformly refined models. The figure shows that the NSRPIM achieves better accuracy and tighter upper bound of
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
z Face z = 1,w = 0.0
Face z = 2, Tz = 100
A
D
B
b
J
C
Face y = 2, Ty = 100
E
Face x = 2, Tx = 100
G
1131
solutions for both uniform and adaptive models than the NS-PIM does, which is consistent with previous studies [11]. Plotting the relative errors of strain energy with respect to the degree of freedom (DOF) in Fig. 11, we can find that the proposed adaptive strategy converges much faster than the uniform model for FEM, NS-PIM and NS-RPIM. Based on the previous study [6], more support nodes are required to create RPIM shape functions and will lead to the increase of bandwidth of system matrix. Considering the time consuming problem and parameters involved in the construction of radial PIM shape functions, we use the simplest linear NS-PIM to study following cases.
a
4.2. Problem 2
F
H I M
K
b
x Face y = 1, v = 0.0
a L
The second problem is a 3-D solid with a crack under uniform tensile force. The dimensions of the problem domain, the loading and boundary conditions applied are presented in the Fig. 12a. Due to symmetry, only one-eighth of model is considered, together with the appropriate boundary and loading conditions are shown in Fig. 12b. Essential boundary conditions are imposed on the faces ABCD, ABEF, BEKJ: on the face ABCD the displacement along the z direction is set
y
Face x = 1, u = 0.0
Fig. 18. Problem 3: a rigid cube without a 1/8 part.
Step 1: 434 Nodes 1662 Tetrahedrons
Step 2: 2011 Nodes 9612 Tetrahedrons
Step 3: 5189 Nodes 27055 Tetrahedrons
Step 4: 7551 Nodes 40485 Tetrahedrons
Fig. 19. Uniform refinement of the problem 3.
Step 1: 434 Nodes 1662 Tetrahedrons
Step 2: 655 Nodes 2674 Tetrahedrons
Step 4: 1683 Nodes 7394 Tetrahedrons
Step 5: 2743 Nodes 12293 Tetrahedrons
Step 3: 1033 Nodes 4422 Tetrahedrons
Step 6: 4503 Nodes 20483 Tetrahedrons
Fig. 20. Meshes of each adaptive step using NS-PIM for the problem 3.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
to zero, on the face ABEF the displacement along the x direction is set to zero and on the face BEKJ the displacement along the y direction is set to zero. Traction boundary condition is uniform tensile force Tz imposed on the face EFHK. The parameters are taken as Young’s modulus E¼1.0, Poisson’s ratio u¼0.3, a¼ 0.5, b¼0.25, c¼1 and Tz ¼1.0. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 153,915 nodes and 852,457 elements) with refinement around the region of crack. This problem is studied by using linear NS-PIM and FEM with both uniform and adaptive meshes. Four uniform refinement models with 562, 2032, 4163 and 7384 nodes are presented in Fig. 13. The adaptive analysis started from uniform model of 562 nodes and is performed for three steps with refinement rate Z ¼0.05. The adaptive mesh and nodal distribution obtained using NS-PIM for each step are plotted respectively in Figs. 14 and 15. These figures illustrate that our error indicator is able to accurately identify the steep gradient of stresses and successfully leads to the refinement of nodes concentrated along the singular line CD, where the crack tip is located. Fig. 16 shows the comparisons on convergence of strain energy between the results of uniform meshes and adaptive meshes. The results show that by using NS-PIM and FEM with adaptive analysis procedure, a much tighter upper and lower bound to the reference energy can be obtained efficiently. It can be found that the NS-PIM has better accuracy than FEM in both uniform and adaptive meshes. As shown in Fig. 17 the relative errors of strain energy are plotted against the DOF for both uniform and adaptive refinements. One can also see that the adaptive models obtain much faster convergence than the uniform models, no matter NS-PIM or FEM.
mesh of 257,404 nodes and 1,318,451 elements) with refinement around the region of stress concentration. An adaptive analysis of this problem is conducted using both NS-PIM and FEM. Four uniform refinement models with 434, 2011, 5189 and 7551 nodes are studied, as presented in Fig. 19. The adaptive analysis started from uniform model of 434 nodes and is performed for six steps with refinement rate Z ¼ 0.05. During the adaptation, problem domain is refined with newly added nodes based on our adaptive strategy, and the adaptive mesh obtained using NS-PIM for each step is plotted in Fig. 20. The figures show that our error indicator can accurately capture the steep gradient of stresses and lead to the refinement of nodes concentrated along the three singular lines GD, GH and GF.
102 Relative Error of Strain Energy (%)
1132
Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
101
103
4.3. Problem 3 The third problem is a rigid cube with 1/8 part removed and subjected to uniform tensile forces, as shown in Fig. 18. Essential boundary conditions are imposed on the faces DEFG, CDGH and HIFG: on the face DEFG the displacement along the x direction is set to zero, on the face CDGH the displacement along the y direction is set to zero and on the face HIFG the displacement along the z direction is set to zero. Traction boundary conditions are uniform forces Tx,Ty,Tz imposed on the faces BCHILM, EFILKJ and ABCDEJ, respectively. The parameters are taken as Young’s modulus E¼ 3.0 107, Poisson’s ratio u ¼0.3, a¼2, b¼ 1, Tx ¼100, Ty ¼100 and Tz ¼100. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral
105
104 DOF Fig. 22. Comparison of convergence for the problem 3.
z
Unit uniform tension y
K
G
b E
F x
a b a D
x 10-3 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
6
Strain Energy
5.5
C
A
Unit uniform tension
B b
5 4.5 J
z
K
4
Face z = 5, Tz = 1
3.5 Face x = 0, u=0
3
I
W
2.5
F
0
4000
8000
12000
16000
20000
24000
DOF Fig. 21. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 3.
Face y = 0, v = 0
y
G
D P
E A
C Face z = 0, w = 0
H R B
x
Fig. 23. Problem 4: a cube with a cube hole in the middle. (a) Whole model and (b) one-eighth of model.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
Fig. 21 shows the comparisons on convergence of strain energy between the results of uniform meshes and adaptive meshes. It can be seen that using all the models, no matter using uniform meshes or adaptive meshes, NS-PIM gives upper bound solutions (in strain energy norm) to the reference one; while FEM gives lower bound. This founding is consistent with our previous study [11,12]. The computed upper and lower bound strain energy obtained using the present adaptive strategy converges much faster than that computed using uniform mesh. In Fig. 22, it is also found that the relative error of strain energy of the adaptive model converges much faster than that of the uniform refinement for both NS-PIM and FEM. 4.4. Problem 4 The last example is the cube domain with a cube cavity in the middle, as illustrated in Fig. 23a. Traction boundary conditions are two uniform tensions Tz on the top face EFGK and on the bottom face ABCD. Due to symmetry, only one-eighth of model is considered, together with the appropriate boundary and loading conditions are
Step 1: 529 Nodes 2131 Tetrahedrons
Step 2: 2156 Nodes 10318 Tetrahedrons
shown in Fig. 23b. The parameters are taken as Young’s modulus E¼3.0 107, Poisson’s ratio u¼0.3, a¼1, b¼10, Tz ¼1. A reference value of strain energy is obtained by ABAQUS using very fine mesh (four-node tetrahedral mesh of 425,931 nodes and 2,430,905 elements) with refinement around the region of cube hole. Adaptive analyses using the NS-PIM and FEM started from the mesh of 529 nodes, and five adaptive steps are performed for both these two methods with refinement rate Z ¼0.05. For comparison, four uniform meshes with 529, 2156, 4165 and 7294 nodes are also computed, as presented in Fig. 24. The adaptive meshes for each step using NS-PIM are plotted in Fig. 25. The figure illustrates that the proposed error indicator is able to accurately capture the steep gradient of stresses and successfully leads to the refinement concentrated along the singular area, i.e. the cube hole. Fig. 26 shows that NS-PIM and FEM provide upper and lower bounds to the exact strain energy, respectively, and the adaptive models provide much tighter bounds than the uniform models. In Fig. 27, one can see that the adaptive models obtain much faster convergence rate than the uniform models for both NS-PIM and FEM.
Step 3: 4165 Nodes 21727 Tetrahedrons
Step 4: 7294 Nodes 39609 Tetrahedrons
Fig. 24. Uniform refinement of the problem 4.
Step 3: 1055 Nodes 4683 Tetrahedrons
1133
Step 1: 529 Nodes
Step 2: 738 Nodes
2131 Tetrahedrons
3110 Tetrahedrons
Step 4: 1586 Nodes 7432 Tetrahedrons
Step 5: 2484 Nodes 11945 Tetrahedrons
Fig. 25. Meshes of each adaptive step using NS-PIM for the problem 4.
1134
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
NS-RPIM achieves better accuracy and tighter upper bound for both uniform and adaptive models than the NS-PIM. (6) Compared with uniform refinements, the proposed adaptive strategy based on NS-PIM improves convergence significantly. It clearly indicates that the presented adaptive procedure can achieve the desired accuracy with far less nodes.
x 10-6 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM) Reference Solution
Strain Energy
2.0885
2.088
Incorporated with an ideal error indicator and effective refinement scheme, the present adaptive procedure based on NS-PIM for 3-D problems has achieved very good performance. All the results have demonstrated that the proposed adaptive procedure is simple, robust and efficient.
2.0875
2.087
0
4000
8000
12000 DOF
16000
20000
24000
Fig. 26. Comparison of convergence of strain energy solutions between adaptive refinement and uniform refinement for the problem 4.
Relative Error of Strain Energy (%)
10-1 Adaptive Mesh (NS-PIM) Uniform Mesh (NS-PIM) Adaptive Mesh (FEM) Uniform Mesh (FEM)
104 DOF
The author wish to thank the support of the China-funded Postgraduates’ Studying Aboard Program for Building Top University and the National Natural Science Foundation of China. The authors also give sincerely thanks to the support of Centre for ACES at National University of Singapore. It is also supported by the program for Changjiang Scholar and Innovative Research Team in University (5311050050037). References
10-2
103
Acknowledgments
105
Fig. 27. Comparison of convergence for the problem 4.
5. Conclusions In this paper, an effective adaptive procedure using the NS-PIM is proposed for 3-D problems. Elasticity problems involving stress concentration and solution singularities have been studied to examine the performance of the proposed adaptive procedure. From the results obtained in the numerical experiments, the following conclusions can be drawn: (1) NS-PIM works well with the simplest four-node tetrahedral elements, which is easy to be implemented for complicated shape. (2) The proposed error indicator works well for all the four numerical cases studied in this work. It effectively captures the location of steep gradient of stress and carried out the refinement in the area of stress concentration. (3) The h-type refinement scheme is performed by adding in new nodes to refine old mesh and the information of old nodes in previous step can be saved. (4) Based on new nodes configuration, an automatic 3-D mesh generator is used to produce high quality tetrahedral meshes. (5) The NS-PIM models, including the NS-PIM with polynomial PIM shape functions and NS-RPIM with radial PIM shape functions, can provide upper bound solutions in energy norm.
[1] Zienkiweicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering 1987;142:337–57. [2] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique. International Journal for Numerical Methods in Engineering 1992;33:1331–64. [3] Lee CK, Lo SH. Automatic adaptive 3D finite element refinement using different-order tetrahedral element. International Journal for Numerical Methods in Engineering 1997;40:2195–226. [4] Lee CK, Lo SH. A full 3D finite element analysis using adaptive refinement and PCG solver with back interpolation. International Journal for Numerical Methods in Engineering 1999;170:39–64. [5] Liu GR, Liu MB. Smoothed particle hydrodynamics-A meshfree practical method. World Scientific; 2003. [6] Liu GR. Meshfree methods: moving beyond the finite element method. 2nd ed. CRC Press; 2009. [7] Liu GR, Zhang GY, Dai KY, Wang YY, Huang HT, Zhong ZH, et al. A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems. International Journal of Computational Methods 2005;2:645–65. [8] Liu GR. A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods. International Journal of Computational Methods 2008;5:199–236. [9] Zhang GY, Liu GR, Wang YY, Huang HT, Zhong ZH, Li GY, et al. A linearly conforming point interpolation method (LC-PIM) for three-dimensional elasticity problems. International Journal for Numerical Methods in Engineering 2007;72:1524–43. [10] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering 2001;50:435–66. [11] Zhang GY, Liu GR, Nguyen TT, Song CX, Han X, Zhong ZH, et al. The upper bound property for solid mechanics of the linearly conforming radial point interpolation method (LC-RPIM). International Journal of Computational Methods 2007;4:521–41. [12] Liu GR, Zhang GY. Upper bound solution to elasticity problems: a unique property of the linearly conforming point interpolation method (LC-PIM). International Journal for Numerical Methods in Engineering 2007;74:1128–61. [13] Zhang GY, Liu GR. An efficient adaptive analysis procedure for certified solutions with exact bounds of strain energy for elasticity problems. Finite Elements in Analysis and Design 2008;44:831–41. [14] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 2: error estimates and adaptivity. International Journal for Numerical Methods in Engineering 1992;33:1365–82. [15] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery (SPR) and adaptive finite element refinement. Computer Methods in Applied Mechanics and Engineering 1992;101:207–24. [16] Babuska I, Rheinboldt WC. A-posteriori error estimates for the finite element method. International Journal for Numerical Methods in Engineering 1978;12:1597–615. [17] Babuska I, Rheinboldt WC. Error estimates for adaptive finite element computations. SIAM Journal on Numerical Analysis 1989;15:739–54.
Q. Tang et al. / Engineering Analysis with Boundary Elements 35 (2011) 1123–1135
[18] Oden JT, Demkowicz L, Rachowicz W, Westermann TA. Toward a universal h-p adaptive finite element strategy. Part 2: a posteriori error estimation. Computer Methods in Applied Mechanics and Engineering 1989;77:113–80. [19] Ainsworth M, Oden JT. A unified approach to a posteriori error estimation using element residual methods. Numerische Mathematik 1993;65:23–50. [20] Bernard BT, Liu GR, Zhang GY, Lu C. A residual based error estimator using radial basis functions. Finite Elements in Analysis and Design 2008;44: 631–45. [21] Babuska I, Suri M. The p- and h-p version of the finite element method, an overview. Computer Methods in Applied Mechanics and Engineering 1990;80:5–20. [22] Zienkiewicz OC, Taylor RL, editors. Oxford: Butterworth-Heinemann; 2000. [V1: The basis]. [23] Chung HJ, Belytschko T. An error estimate in the EFG method. Computational Mechanics 1998;21:91–100. [24] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: part I: stress recovery and a posteriori error estimation. Computers & Structures 2004;82:413–28. [25] Lee CK, Zhou CE. On error estimation and adaptive refinement for element free Galerkin method: part II: adaptive refinement. Computers & Structures 2004;82:429–43. [26] Durate CA, Oden JT. An h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 1996;139:237–62. [27] Belytschko T, Liu WK, Singer M. On adaptivity and error criteria for meshfree methods. Advances in adaptive computational methods in mechanics. Elsevier Science 1998;47:217–28. [28] Gan NF, Li GY, Long SY. 3D adaptive RKPM method for contact problems with elastic–plastic dynamic large deformation. Engineering Analysis with Boundary Elements 2009;33:1211–22. [29] Liu GR, Tu ZH. An adaptive procedure based on background cells for meshless methods. Computer Methods in Applied Mechanics and Engineering 2002;191:1923–43. [30] Rabczuk T, Belytschko T. Adaptive for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering 2005;63:1559–82. [31] Park SH, Kwon KC. A posteriori error estimates and an adaptive scheme of least-square meshfree method. International Journal for Numerical Methods in Engineering 2003;58:1213–50. [32] Angulo A, Pozo L, Perazzo F. A posteriori error estimator and an adaptive technique in meshless finite points method. Engineering Analysis with Boundary Elements 2009;33:1322–38.
1135
[33] Timoshenko SP, Goodier JN. Theory of elasticity.3rd ed.. NewYork: McGrawhill; 1970. [34] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering 2001;50:937–51. [35] Li Y, Liu GR, Luan MT, Dai KY, Zhong ZH, Li GY, et al. Contact analysis for solids based on linearly conforming radial point interpolation method. Computational Mechanics 2007;39(4):537–54. [36] Liu GR, Li Y, Dai KY, Luan MT, Xue W. A linearly conforming radial point interpolation method for solid mechanics problems. International Journal of Computational Methods 2006;3:401–28. [37] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. International Journal for Numerical Methods in Engineering 2002;54:1623–48. [38] Liu GR, Zhang GY. Edge-based smoothed point interpolation methods. International Journal of Computational Methods 2008;5(4):621–46. [39] Liu GR. A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: part I theory. International Journal for Numerical Methods in Engineering 2010;81. 1093-26. [40] Hang S. Adaptive tetrahedral mesh generation by constrained Delaunay refinement. International Journal for Numerical Methods in Engineering 2008;75:856–80. [41] Frey PJ, George Pl. Mesh generation: application to finite elements. Oxford, U.k: Hermes Science; 2000. [42] Weatherill NP, Hassan O. Efficient three-dimentional Delaunay triangulation with automatic point creation and imposed boundary constraints. International Journal for Numerical Methods in Engineering 1994;37:2005–39. [43] Chew PL. Guaranteed-quality Delaunay meshing in 3D. In: Proceedings of the 13th annual ACM symposium on computational geometry. Nice, France;1997. p. 391–93. [44] Miller GL, Talmor D, Teng SH, Walkington N, Wang H. Control volume meshes using sphere packing: generation, refinement and coarsening. In: Proceeding of 5th international meshing roundtable; 1996. p. 47–61. [45] Shewchuk JR. Delaunay refinement mesh generation. PhD thesis. Department of Computer Science, Carnegie Mellon University, Pittsburgh, PA. Available as Technical report CMU-CS-97-137; 1997. [46] Shewchuk JR. Triangle: engineering a 2d quality mesh generator and delaunay triangulator. Applied Computational Geometry: Towards Geometric Engineering 1996;1148:203–22. [47] ABAQUS Theory Manual and User Manual, version 6.4. Hibbit, Karlsson & Sorensen Inc., USA; 2004.