- Email: [email protected]
Smoothed Galerkin methods using cell-wise strain smoothing technique
Engineering Analysis with Boundary Elements 36 (2012) 825–835
Contents lists available at SciVerse ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Smoothed Galerkin methods using cell-wise strain smoothing technique X.Y. Cui, G.Y. Li n State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha, 410082, PR China
a r t i c l e i n f o
abstract
Article history: Received 14 October 2011 Accepted 21 December 2011 Available online 11 January 2012
A smoothed Galerkin method (SGM) using cell-wise strain smoothing operation is formulated in this paper. In present method, the field nodes can be divided into two types: boundary field nodes and interior field nodes. The background cells are divided into SC smoothing cells, and the strains in each smoothing cell are obtained using a gradient smoothing technique which can avoid evaluating derivatives of shape functions at integration point. The field variables of boundary points are approximated using linear interpolation of neighbour boundary field nodes, and the shape functions possess the Kronecker Delta property and facilitate the impositions of essential boundary conditions. The field variables of interior points are approximated using moving least-squares approximation using the support field nodes around them. A number of numerical examples are studied and confirm the significant features of the present methods: (1) can pass the standard patch test; (2) can easily impose essential boundary conditions as those in finite element method; (3) can avoid evaluating derivatives of shape functions; (4) no numerical parameter is required. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Meshfree Smoothed Galerkin method Gradient smoothing operation Cell-wise strain smoothing
1. Introduction In the past decades, various meshfree methods including smooth particle hydrodynamics [1], diffuse element method (DEM) [2], element free Galerkin (EFG) method [3], reproducing kernel particle method (RKPM) [4], finite point method (FPM) [5], H–P clouds [6], meshless local Petrov–Galerkin (MLPG) method [7], point interpolation method (PIM)[8], radial point interpolation method (RPIM) [9], among others, have been proposed and applied in more and more fields of particular engineering and scientific problems [10]. One of the best known meshfree methods is the EFG method, which uses the moving least-squares (MLS) approximation [11] for the interpolation of the displacement functions. The EFG method is a very promising method for the treatment of the partial different equations and has been applied and improved in a large variety of problems. However, there are some major problems in using EFG such as the imposition of essential boundary conditions, the complex computation in evaluating derivatives of meshfree shape functions, the uncertain numerical parameters like the nodal influence radius, the minimization of the computational cost and the implementation of an appropriate numerical integration scheme. Some works have been developed for the alleviation of the above-mentioned problems. The essential boundary conditions can not be imposed directly, because the MLS is not a true
n
Corresponding author. Tel.: þ86 731 88821717; fax: þ 86 731 88822051. E-mail address: [email protected] (G.Y. Li).
0955-7997/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.12.013
interpolation and the shape function constructed by it lacks the Kronecker delta property. Belytschko et al. [3] used the Lagrange multipliers to impose essential boundary conditions. Collocation methods [12] and penalty method [13] were also used to solve this problem. Other methods have been developed to overcome this problem including: the modification of the variational principle [14], the usage of the singular weighting functions [15], and the combination with finite elements [16]. Some coupling schemes were also investigated in the recent years. A coupled FE–EFG method was proposed by Belytschko et al. [17]. Gu and Liu proposed a coupled EFG–BEM for the analysis of two-dimensional solids [18], a natural neighbour-based MLS approach is developed coupled EFG and NEM method [19]. In this paper, the shape functions of field node local on the boundary of the problem domain are enforced to satisfy the Kronecker delta property, and the field variables of the point on the boundary are interpolated using the neighbour field nodes on the boundary, the essential boundary conditions can be imposed directly as those in FEM. Due to the complexity involved using Gauss integration in EFG methods, Beissel and Belytschko [20] proposed a nodal integration procedure of EFG. The unstabilized EFG with nodal integration displays spurious near-singular modes in some problems. To overcome this problem, Chen et al. [21] proposed a stabilized conforming nodal integration using a strain smoothing technique, and this method can reproduce a field exactly. Based on the strain smoothing technique, Liu et al. [22] proposed a smoothed finite element method (SFEM). The SFEM further divides the elements into some smoothing cells, computes the integrals along the edge of the smoothing cells, and has been proven to have good
826
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
properties. Liu et al. gave detailed theoretical aspects including stability, bound property and convergence about SFEM [23]. In present work, the idea of dividing the elements into some smoothing cells was further studied and used in the smoothed Galerkin method. In this paper, SGM is formulated using cell-wise strain smoothing operation in each background integration cells. In SGM, the field nodes are divided into two types: boundary field nodes and interior field nodes. The node shape functions of the boundary field nodes possess the Kronecker Delta property and the interior nodal variables are approximated using MLS approximation. The support nodes of each node for MLS approximation are selected based on the background cells. Depending on the requirement of the accuracy and stability, the background cells can be divided into several smoothing cells, and the strains in each smoothing cell are obtained using a gradient smoothing technique which can avoid evaluating derivatives of meshfree shape functions at integration point. The shape functions are directly created using linear interpolation for point on the boundary of the problem domain and MLS approximation for inner points. The support nodes selections of these points are also base on the background cells. For integration points in the background cells and integration points on the edges of the background cells, different support-node selections will be adopted to guard the continuity of the displacement field. The strains in the smoothing cells can be obtained as constants, hence there need not any other numerical integration in forming the stiffness matrix. To show the performance of the proposed method, a series of benchmark examples were presented, and excellent results have been obtained illustrating the efficiency and accuracy of the present methods. The paper is outlined as follows. Section 2 briefly introduces the basic equations of linear elasticity. Gradient smoothing method is introduced in Section 3. In Section 4, the shape function construction is presented in detail. Section 5 gives the presentation of the stiffness matrix and numerical examples are presented in Section 6. Conclusions and remarks are made in Section 7.
In Eq. (3), eT ¼ {exx, eyy, 2exy} is the vector of strains that relates to the displacements by the following compatibility equation
e ¼ Lu where u¼ {ux, uy}T is the displacement vector. Boundary conditions are given as follows: u ¼ uG on Gu
ð6Þ
nT r ¼ tG on Gt
ð7Þ
where uG is the specified displacement on the essential boundary
Gu, tG is the given traction on the natural boundary Gt, and n is the unit outward normal matrix expressed as 2 3 nx 0 6 7 n ¼ 4 0 ny 5 ny nx
O
@ @y
O
Gt
where de ¼L(du).
3. Gradient smoothing method The gradient smoothing method was proposed by Chen et al. [21], it can avoid evaluating derivatives of meshfree shape functions at point of interest and thus eliminates spurious modes. In present work, each background cell is further divided into some non-overlapping smoothing cells Ok ¼ [SC C ¼ 1 OkðCÞ , as shown in Fig. 1. The strain field in the Cth smoothing cell Ok(C) can be treated as a constant using following equation: Z Z ekðCÞ ¼ eðxÞfC ðxÞdO ¼ LuðxÞfC ðxÞdO ð10Þ Ok
Ok
ð1Þ
where fC ðxÞ is a given smoothing function that satisfies at least unity property Z fC ðxÞdO ¼ 1 ð11Þ
ð2Þ
A constant smoothing function is adopted as follows: ( 1=AC x A OkðCÞ fC ðxÞ ¼ 0 x= 2OkðCÞ
A solid mechanics problem of static elasticity can be described by equilibrium equation in the domain O which can be given by
where L is differential operator matrix defined as 2@ 3 0 @x 60 @ 7 L¼6 @y 7 4 5
ð8Þ
By applying principle of virtual work, a weak form is obtained for the equilibrium equation and the boundary conditions as follows: Z Z Z deT rdO duT bdO duT tG dG ¼ 0 ð9Þ
2. Brief of basic equation
LT r þ b ¼ 0
ð5Þ
@ @x
OkðCÞ
ð12Þ
where AC is the area of the smoothing cell Ok(C).
rT ¼{sxx, syy, sxy} is the stress vector and bT ¼{bx, by } is the body force vector. The stresses relate the strain via the constitutive equation as follows:
r ¼ De in which D is follows: 2 1 E 6n D¼ 4 1n2 0
ð3Þ the matrix of material constants that defined as
n 1 0
0
3
0 7 5
ð4Þ
1v 2
where E is Young’s modulus and n is Poisson’s ratio.
Fig. 1. The background cell Ok is divided into SC smoothing cells, Ok(C) is the Cth smoothing cell with the boundary Gk(C).
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
minimizing the weighted residual as follows:
Substituting Eq. (12) into Eq. (10), the smoothed strain field can be rewritten as Z Z 1 1 ekðCÞ ¼ LuðxÞdO ¼ nUuðxÞdG ð13Þ AC OkðCÞ AC GkðCÞ
JðxÞ ¼
NP X
NI ðxÞuI
ð14Þ
where NP is the number of the support nodes of point x, uI ¼{uI, vI}T is the SGM approximation for the displacement vector at node I, and NI(xi) is the SGM shape function introduced in the next section.
aðxÞ ¼ A1 ðxÞ BðxÞ u
ð19Þ
where u is the vector that collects the nodal parameters of displacement field in the support domain, and A(x) is called the weighted moment matrix defined by
4.1. Moving least-squares approximation
pj ðxÞaj ðxÞ ¼ pT ðxÞaðxÞ
n X
AðxÞmm ¼
The moving least-squares (MLS) approximation of a function u(x) denoted by uh(x) [3] is expressed as follows:
wðxxI ÞpðxI ÞpT ðxI Þ
BðxÞmn ¼ wðxx1 Þpðx1 Þ
ð15Þ
wðxx2 Þpðx2 Þ
wðxxn Þpðxn
ð21Þ
where pðxÞ ¼ 1
ð20Þ
I¼1
j¼1
ð17Þ
where r is the radius of the support of the weight function w(r). The stationary of functional J with respect to a(x) leads to:
4. Shape function construction for SGM
m X
wðr I Þ½pT ðxI ÞaðxÞuI 2
in which n is the number of nodes in the support domain of the interested point x, uI is the nodal value of u at node I, rI is the distance between two point x and xI, and w(r) is the weight function, for example, 8 2 3 > 2 > 4 rr þ 4 rr for rr r 12 > 3 > > < 2 3 ð18Þ wðrÞ ¼ 4 4 r þ 4 r 4 r for 12 o rr r 1 >3 3 r r r > > > > :0 for rr 41
I¼1
uh ðxÞ ¼
n X I¼1
where Gk(C) is the boundary of the smoothing cell Ok(C), n is the outward normal matrix containing the components of the outward normal vector on the boundary Gk. In the present method, the displacements of point x can be interpolated using the nodal displacements at the support nodes as follows: uh ðxÞ ¼
827
x
y
T
Substituting Eq. (19) into Eq. (15), the MLS approximation can be obtained
ð16Þ
uh ðxÞ ¼
is the basis vector in 2D, and a(x) is the corresponding coefficient vector. The coefficients a(x) are obtained at any point x by
n X
jI ðxÞuI ¼ uT ðxÞu
I¼1
Boundary
n1 SC3
xe 4
xb1
SC1
SC1
xi1 xb 2
n2
SC2
xi 3
SC2
SC4
SC3
xi 2
x e1
xe3
xe2 Boundary field node
SC1
Interior field node
Integration point in cell
Integration point on the edge of the cell
SC2 SC3
Integration point on the boundary
Fig. 2. Distribution of Integration points for a background divided into smoothing cells.
ð22Þ
828
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
where jI(x) is the MLS shape function defined by
jI ðxÞ ¼
m X
5. Discretized system equations
pj ðxÞ½A1 ðxÞBðxÞjI
ð23Þ
j¼1
4.2. SGM shape function In the present work, the problem domain is first discretized with background cells (triangular, quadrangular or n-side polygon), then further divided into SC smoothing cells. Though the smoothing operation, only shape functions of the integration points on the boundary of the smoothing cells will be evaluated. The shape functions of SGM are created using linear interpolation for points on the boundary of the problem domain, and using MLS approximation for the point in the problem domain. As shown in Fig. 2, for the integration point on the boundary of the problem domain, such as xb1and xb2, shape functions of them can be created by linear interpolation of nodes n1 and n2. For the integration point in the problem domain, such as xe1 xe4 and xi1 xi3, MLS shape functions are established, and the support-node selections are base on background cells. For the integration point in the background cells and the integration point on the edge of the background cell, different support-node selections will be adopted to guard the continuity of the displacement field. For the point of interest xi located in a cell i, the support node contains all nodes of the background cell which has common node of the cell i. As shown in Fig. 3, i1–i16 are support nodes for point xi. For the point of interest xe local on an edge (only for interior edge) k, support nodes contain all nodes of the background cells which have common node of that edge k. For example, e1–e12 are support nodes for point xe. For the point of interest xn is a field node (only for interior field node), support nodes contain all nodes of the background cells which contain the node n. For example, n1–n9 are support nodes for point xn. This scheme can usually select enough nodes and leads to less time consuming than the traditional scheme, in which circular support domain is used. This also scheme can be used to create MLS shape functions with high order of consistence and for extremely irregularly distributed nodes.
As shown in Fig. 4, Gk(C) is the boundary of smoothing cell Ok(C) divided into NG segments, and nJ is the outward normal matrix of the Jth segment. Substituting Eq. (14) into Eq. (13), the smoothed strain field can be rewritten by ! ! Z NG NP NP X 1 X 1 X ekðCÞ ¼ nUNI ðxÞdG uI ¼ nJ UNI ðxJ ÞlJ uI AC I ¼ 1 AC I ¼ 1 J ¼ 1 GkðCÞ I
¼ BkðCÞ uI
ð24Þ
where lJ is the length of Jth segment, xJ is the coordinate vector of I the gauss point on the Jth segment, and BkðCÞ is the smoothed strain displacement matrix corresponding to field node I. 2 3 " # nxJ 0 NG NG X X 0 1 1 I 6 0 n 7 NI ðxJ Þ yJ 5U BkðCÞ ¼ nJ UNI ðxJ ÞlJ ¼ l 4 0 NI ðxJ Þ J AC J ¼ 1 AC J ¼ 1 nyJ nxJ 2 3 N I ðxJ ÞnxJ 0 NG X 1 6 0 NI ðxJ ÞnyJ 7 ¼ ð25Þ 4 5lJ AC J ¼ 1 NI ðxJ ÞnyJ NI xJ nxJ where nJx and nJy are the components of the outward unit normal to the Jth boundary segment. One assumes that the material parameters E, u are same over a smoothing cell, the smoothed stress in smoothing cell Ok(C) can be easily obtained by
rkðCÞ ¼ DekðCÞ
ð26Þ
We now seek for a weak form solution of displacement field u that satisfies the following smoothed Galerkin weak form Z Z Z deT rdO duT bdO duT tG dG ¼ 0 ð27Þ O
O
Substituting Eqs. (14), (24) and (26) into Eq. (27), a set of discretized algebraic system equations can be obtained in the following matrix form: Kdf ¼ 0
ð28Þ
where f is the force vector defined as Z Z f ¼ UT ðxÞbdO UT ðxÞtG dG O
e
e
e
e
e
edge k e
x
e
e
i
i
i
i
e
e
i
i
i
i
e
cell i i
i
i
i
i
i
i
i
e
n
n
n
x
Gt
in which UT(x) is the shape function vector to build force vector. In Eq. (28), K is the (global) smoothed stiffness matrix assembled in the form of KIJ ¼
N cell X
KIJðkÞ
ð30Þ
k¼1
where the summation means an assembly process same as the practice in the FEM, Ncell is the number of the number of the background cells of the whole problem domain O, and KIJðkÞ is
node n
n
n
n2
n
l2 2
x n
n
n
n3 3 smoothing
Boundary field nodes
Fig. 3. Support-node selection in SGM for shape function construction. For interest point xi in a cell i, the support nodes are i1 i16. For interest point xe on an edge k, nodes e1 e12 are selected as the support nodes. For interest point xn as field node n, nodes n1 n9 are selected as the support nodes.
Γ
l1
domain Ω
l3 Interior field nodes
ð29Þ
Gt
1
l4
n1
4
n4
gauss point
Fig. 4. A smoothing cell Ok(C) with the boundary Gk(C), which is consist of four segments. lJ is the length of the Jth segment, nJ is the outward normal vector of the Jth segment.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
the stiffness matrix associated with Ok that is computed using KIJðkÞ ¼
SC Z X C¼1
OkðCÞ
I
J
ðBkðCÞ ÞT DðBkðCÞ ÞdO ¼
SC X
I
J
ðBkðCÞ ÞT DðBkðCÞ ÞAC
ð31Þ
C¼1
in which SC is number of the smoothing cells, AC and BkðCÞ are the area and smoothed curvature-deflection matrix of the Cth smoothing cell, respectively.
6. Numerical examples In this section, examples for 2D elasto-statics analyses are presented. For MLS approximation, linear basis vector (m¼ 3 in Eq. (15)) and weight function w given in Eq. (18) are used for constructing the shape functions. For all numerical examples, quadrilateral cells with one and four smoothing cells, triangular cells with one and three smoothing cells are adopted. For convenience, triangular cells with one smoothing cell for SGM will be named as SGM-TSC1; quadrilateral cells with four smoothing cells for SGM will be named as SGM-QSC4; etc.
6.1. Patch test The first numerical example is the standard displacement patch test; the mesh of a square parch is depicted in Fig. 5. A prescribed linear displacement field is imposed on the boundary Gu of the patch are computed using u ¼ 0:001ðx þyÞ v ¼ 0:001ðxyÞ
829
displacement is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uNnode uP exact unum ÞT ðuexact unum Þ u I I I u I ¼ 1 ðuI u Ed ¼ u NP node t ðuexact ÞT ðuexact Þ I I
ð33Þ
I¼1
where the superscript num denotes the displacement vector obtained using numerical methods, exact denotes the exact or analytical solution, and Nnode is the number of total field nodes. Table 1 lists the displacement norm errors of the numerical results for standard displacement patch test using both regularly and irregularly distributed field nodes shown in Fig. 5. All the numerical models can pass of the patch test to the machine accuracy and hence are capable to reproduce linear displacement field ‘‘exactly’’. 6.2. Cantilever beam problem A cantilever beam with length L and height D is studied as benchmark problem to test the convergence of the method, which is subjected to a parabolic traction at the free end as shown in Fig. 6. The beam is assumed to have a unit thickness so that plane stress condition is valid. The analytical solution from Timoshenko and Goodier [24] can be given by " !# 2 Py 2 D ux ðx,yÞ ¼ ð6L3xÞx þ ð2þ nÞ y 6EI 4
Table 1 Displacement error norm of numerical results for the standard patch test.
ð32Þ
The material properties of patch are E¼1.0 and n ¼0.25. To satisfy the patch test, the linear displacement field approximated exactly, and the numerical solution at any interior nodes should be in exact agreement the analytic ones given in Eq. (32). To examine the numerical error precisely, an error norm in
SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Regularly distributed nodes
Irregularly distributed nodes
1.48431617E 14 9.81195622E 15 1.79482393E 14 1.17421298E 14
1.33104257E 14 1.49256754E 14 1.00570811E 14 1.08551760E 14
Fig. 5. Displacement patch test: (a), (b) and (c) uniform grid; (d), (e) and (f) irregular grid.
830
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
" # P D2 x 2 2 3ny ðLxÞ þ ð4 þ5nÞ þð3LxÞx uy ðx,yÞ ¼ 6EI 4
ð34Þ
D2 y2 4
!
D x
L Fig. 6. Cantilever beam subjected to a parabolic traction on the right tip.
-1
-0.5 SGM-TSC1 SGM-TSC3 FEM, T3
-1
SGM-QSC1 SGM-QSC4 FEM, Q4
-1.5
-2
Log10(Ed)
Log10(Ed)
-1.5
-2
-2.5
-3
-3.5 -2.5 -4
-3 -0.4
-0.2
0 0.2 Log10(Δh)
0.4
-4.5 -0.4
0.6
-0.2
0 0.2 Log10(Δh)
0.4
0.6
Fig. 7. Convergence of the numerical results in displacement error norm for the problem of cantilever beam solved using different method.
-3
-3.4 SGM-QSC1
SGM-TSC1 -3.2
-3.6
SGM-TSC3
SGM-QSC4 FEM, Q4
FEM, T3
-3.8
-3.4
-4 Log10(Ee)
Log10(Ee)
-3.6 -3.8
-4.2 -4.4
-4 -4.6 -4.2
-4.8
-4.4 -4.6 -0.4
-5
-0.2
0 0.2 Log10(Δh)
0.4
ð35Þ
where I¼D3/12 is the moment of inertia for the beam. The related parameters for the problem are: E¼3.0 107 Pa, n ¼0.3, L¼48 m, D ¼12 m and P¼ 1000 N. Using the same set of distributed field nodes, the cantilever beam is studied using the proposed methods with triangular and quadrilateral background cells, respectively. Further, this problem is also studied using FEM triangular element (T3) and quadrilateral (Q4) with the same distributed field nodes. In the numerical computations, the displacement boundary on the left uses the exact displacements obtained from Eq. (34) and the loading boundary uses the distributed shear stresses in Eq. (35).
y
P
Py P sxx ðx,yÞ ¼ ðLxÞ, syy ðx,yÞ ¼ 0, sxy ðx,yÞ ¼ I 2I
0.6
-5.2 -0.4
-0.2
0 0.2 Log10(Δh)
0.4
0.6
Fig. 8. Convergence of the numerical results in energy error norm for the problem of cantilever beam solved using different method.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
Fig. 7 plots the convergence of the solutions in displacement error norm for the cantilever beam solved using different methods. The mesh size parameter Dh is taken to be the average nodal spacing. For the left figure of Fig. 7, all methods use the triangular cells. It can be found that SGM-TSC1 and SGM-TSC3 are more accurate than FEM (T3). For the right figure of Fig. 7, quadrilateral cells are used for presented methods. Ones find that SGM-QSC4 gives better displacement results than FEM (Q4). SGM-QSC4 not only gives the best result, but also has the highest rate of convergence 2.43, which is much larger than the theoretical value 2.0. Fig. 8 shows the convergence of solutions in energy norm for the cantilever beam solved using different methods. The energy
error norm is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uR ðeexact enum ÞT Dðeexact enum ÞdO 1u Ee ¼ t O R exact ÞT Dðeexact ÞdO A O ðe
ð36Þ
where A is the area of the problem domain. It is found that SGMTSC1 and SGM-TSC3 have better accuracy and converge faster compared to the FEM (T3). SGM-TSC1 provides good accuracy in energy and a high convergence rate (1.193). Form the right figure of Fig. 8, it can be found that all of the SGM have better accuracy and higher rate of convergence than FEM (Q4). It is observed that SGM-QSC1 presents superconvergent property and the
5
4.8 SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
4.8
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
4.75 4.7
4.6
4.65 Strain energy
Strain energy
831
4.4
4.2
4.6 4.55 4.5 4.45
4
4.4 3.8 4.35 3.6 0
1000
2000 3000 DOF
4000
4.3
5000
0
1000
2000 3000 DOF
4000
5000
Fig. 9. Solutions (strain energy) converging to the exact solution for the problem of cantilever beam solved using different method.
x 10-3
-8.6
Vertical displacement of middle tip node
-7.6 -7.8 -8 -8.2 -8.4 -8.6 -8.8 -9 -9.2
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
-8.7 Vertical displacement of middle tip node
SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
-7.4
x 10-3
-8.8 -8.9 -9 -9.1 -9.2 -9.3 -9.4
-9.4 0
1000
2000 3000 DOF
4000
5000
-9.5
0
1000
2000 3000 DOF
4000
5000
Fig. 10. Solutions (vertical displacement of middle tip node) converging to the exact solution for the problem of cantilever beam solved using different method.
832
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
convergence rates of them are 1.65, which are much higher than theoretical value 1.0. To study the convergence property, the strain energy of the cantilever beam is computed using different methods and plotted
1000 800 600
Analytical solution SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Normal stress σx
400 200 0 -200 -400 -600 -800 -1000 -6
-4
-2
0 y
2
4
6
Fig. 11. Normal stress sx distribution along the cross-section of the cantilever at x¼ L/2 obtained using different methods and the same distributed 637 nodes.
Analytical solution SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Shear stress σxy
-40 -60 -80 -100 -120 -140 -6
-4
-2
0 y
6.3. Infinite plate with a circular hole An infinite plate with a circular hole (a) subjected to a unidirectional tensile load (p) in the x direct is studied. Owing to the symmetry, only one quarter is model as shown in Fig. 13. Symmetry conditions are imposed on the left and bottom edges, and the inner boundary of the hole is traction free. The analytical solution for the stress can be given in [24]
a2 3 3a4 cos 2y þcos 4y þ 4 cos 4y sxx ¼ p 1 2 2 r 2r 2
a 1 3a4 cos 2ycos 4y þ 4 cos 4y syy ¼ p 2 2 r 2r 2
a 1 3a4 sin 2ysin 4y 4 sin 4y sxy ¼ p 2 ð37Þ 2 r 2r
0 -20
in Fig.9. In the left figure, it is easy to see that SGM-TSC1 has the upper bound solutions and SGM-TSC3 has lower bound solution. From the results, ones find that the numerical model is softer than the exact model when 1 smoothing cell is used for SGM, the strain energy of SGM-TSC1 is a litter larger than the exact solution. The model of SGM-TSC3 is stiffer than the exact model but softer than FEM (T3), and it gives more accurate results than FEM (T3) but still a litter lower than exact solution. Form the right figure of Fig. 9, it can be found that all methods present lower and upper bound property. For the present SGM, the numerical model becomes stiff when the number of smoothing cell increases. SGM-QSC4 presents best accurate than other schemes plotted in the figure. It is very close to exact solution and has the lower bound property. The vertical displacements of middle tip node of the cantilever beam computed using different methods are plotted in Fig. 10. The results again show that SGM-TSC1 and SGM-TSC3 can give better results than FEM (T3) with the same mesh. We find that the convergence trend of displacements is same as it of strain energy. SGM-TSC4 gives the best results than other schemes plotted in the figure. Figs. 11 and 12 illustrate the distribution of the normal stress sx and shear stress sxy on the cross-section x¼ L/2 of the beam. It can be found that the results obtained using present methods agree well with the analytical solutions. SGM-QSC4 exhibits best results compared with other presented numerical methods.
2
4
6
Fig. 12. Shear stress sxy distribution along the cross-section of the cantilever at x¼ L/2 obtained using different methods and the same distributed 637 nodes.
where r and y are the polar coordinates and y is measured counterclockwise from the positive x-axis. The displacement fields can be calculated as follows:
a4 p kþ1 a2 ur ¼ r þcos 2y þ 1 þðk þ1Þcos 2y 3 cos 2y 4G 2 r r 2 4 p a a rð1kÞ þ 3 sin 2y ð38Þ uy ¼ 4G r r
y y
p
p
b
r a
θ
b
x
x Fig. 13. Infinite plate with a circular hole subjected to uniform tensile and its quadrant model.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
833
Fig. 14. Domain discretization for the infinite plate with a hole (9 9 nodes).
-0.5
-1
-1.5
SGM-TSC1 SGM-TSC3 FEM, T3
-1
SGM-QSC1 SGM-QSC4 FEM, Q4
-1.5
Log10(Ed)
Log10(Ed)
-2
-2.5
-2
-2.5
-3 -3 -3.5
-4 -1.2
-3.5
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
-4 -1.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
Fig. 15. Convergence of the numerical results in displacement error norm for the problem of infinite plate with a circular hole solved using different method.
where G¼
E , 2ð1þ nÞ
(
k¼
34n n
1n
for plane stress for plane strain
ð39Þ
E is Young’s modulus and n is Poisson’s ratio. We studied the problem under plane strain conditions and traction boundary conditions are imposed on the upper and right edges with the analytical stresses obtained using Eq. (37). The parameters for this problem are: a¼1 m, b¼5 m, p ¼1 Pa, E¼1 103 Pa and n ¼0.3. The 9 9, 17 17, 25 25, 33 33 and 49 49 distributed nodes are used for computation and a 9 9 node distribution is shown in Fig. 14. Fig. 15 plots the convergence of the solutions in displacement error norm for the problem of infinite plate with a circular hole solved using different methods. Form the left figure of Fig. 15, we find that SGM-TSC3 exhibits the best accuracy compared with other methods using the triangular cells and SGM-TSC1 has the highest convergence rate (2.135). For the methods using quadrilateral background cells, the solutions of all the methods converge well and exhibit different accuracy and convergence rate. In terms of accuracy, SGM-QSC1 is found worse than the FEM (Q4), and all other schemes are at least better than the FEM (Q4). SGM-QSC4 still gives good accuracy as the previous case. It must be pointed that the convergence rate of SGM-Q1 is close to 3.264, but the accuracy of it is not good enough.
Fig. 16 plots the convergence of the solutions in energy error norm for the problem of infinite plate with a circular hole solved using different methods. The methods using both the triangular and quadrilateral background cells are found converging well. SGM-TSC1 again gives the highest convergence rate (1.211) in energy norm as in displacement norm, and it is higher than the theoretical value 1.0. Form the right figure of Fig. 16, it can be found that all of the SGM have better accuracy and higher rate of convergence than FEM (Q4). SGM-QSC1 still present superconvergent property and the convergence rates of them are 2.824, which are more than twice higher than the theoretical value 1.0. The convergence trends of strain energy of the infinite plate with a circular hole computed using different methods are plotted in Fig. 17. SGM-TSC1 still has the upper bound solutions and SGM-TSC3 has lower bound solution. For quadrilateral background cells, SGM-QSC4 presents lower bound property, and SGM-QSC1 presents upper bound property.
7. Conclusion In this work, SGM is proposed using cell-wise strain smoothing operation. The background cells are divided into several smoothing cells and the strain fields in them are projected to a set of constant fields. The strain smoothing stabilization technique is
834
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
-3.6
-2.5 SGM-TSC1 SGM-TSC3 FEM, T3
-3.8
SGM-QSC1 SGM-QSC4 FEM, Q4
-3
-4
-4.2
Log10(Ee)
Log10(Ee)
-3.5
-4.4
-4
-4.5 -4.6 -5
-4.8
-5 -1.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-5.5 -1.2
-0.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
0.01183
0.01185
0.01182
0.01184
0.01181
0.01183
0.01180
0.01182
0.01179
0.01181
Strain energy
Strain energy
Fig. 16. Convergence of the numerical results in energy error norm for the problem of infinite plate with a circular hole solved using different method.
0.01178 0.01177 0.01176
0.01179 0.01178
0.01175
SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
0.01174 0.01173
0.01180
0
1000
2000 3000 DOF
4000
5000
0.01177
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
0.01176 0.01175
0
1000
2000 3000 DOF
4000
5000
Fig. 17. Solutions (strain energy) converging to the exact solution for the problem of infinite plate with a circular hole solved using different method.
employed to restore the conformability and to improve the accuracy and the rate of convergence. Only line integrations along the smoothing cells are needed and no derivative of the shape functions is involved in computing the field gradients to form the stiffness matrix. The support-node-selection for the MLS approximation is based on the background cells and no numerical parameter like the nodal influence radius is required. The results of investigated numerical examples illustrate the effectivity of the present method. Besides the good accuracy, stabilization and high convergence rate, the following features of the present method are noteworthy: (1) can pass the standard patch test; (2) can directly impose essential boundary conditions as those in finite element method; (3) can avoid evaluating derivatives of meshfree shape functions; (4) no uncertain numerical parameter like the nodal influence radius is required.
Acknowledgments This work was supported by National Science Foundation of China (11002053), National Basic Research Program of China (2010CB328005), Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body (51075005). The authors also give sincere thanks to Prof. Liu G.R. at University of Cincinnati for his help. References [1] Monaghan JJ. An Introduction of SPH. Comput. Phys. Commun. 1982;48:89–96. [2] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput. Mech. 1992;10(5):307–18.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
[3] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int. J. Numer. Methods Eng. 1994;37(2):229–56. [4] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. Int. . Numer. Methods Eng. 1995;20:1081–106. [5] Onate Idelsohn ES, Zienkiewicz OC, Tylor RL. A finite point method in computational mechanics, applications to convective transport and fluid flow. Int. J. Numer. Methods Fluids 1996;39:3839–66. [6] Duarte CA, Oden JT. An H–P adaptive methods using clouds. Comput. Methods Appl. Mech. Eng. 1996;139:237–62. [7] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput. Mech. 1998;22(2):117–27. [8] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. Int. Numer. Methods Eng. 2001;50:937–51. [9] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int. J. Numer. Methods Eng. 2002;54:1623–48. [10] Liu GR. Meshfree methods: Moving Beyond the Finite Element Method. Boca Raton, FL: CRC Press; 2002. [11] Lancaster P, Salkauskas K. Surface generated by moving least squares methods. Math. Comput 1981;37:141–58. [12] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin method. Int. . Numer. Methods . Eng. 1996;39:141–58. [13] Zhu T, Atluri N. Modified collocation method and a penalty function for enforcing essential boundary conditions in the element free Galerkin method. Comput. Mech. 1998;21:211–22.
835
[14] Lu YY, Belytschko T, Gu L. A new implementation of the element-free Galerkin method. Comput. Methods in Appl. Mech. Eng. 1994;113:397–414. [15] Kaljevic I, Saigal S. An improved element free Galerkin formulation. Int. J. Numer. Methods Eng. 1997;40:2953–74. [16] Krongauz Y, Belytschko T. Enforcement of essential boundary conditions in meshless approximations using finite elements. Comput. Methods Appl Mech Eng. 1996;131:133–45. [17] Belytschko T, Organ D, Krongauz Y. A coupled finite element-element-free Galerkin method. Comput. Mech. 1995;17:186–95. [18] Gu YT, Liu GR. A coupled element free Galerkin/boundary element method for stress analysis of two-dimensional solids. Comput. Methods Appl. Mech. Eng. 2001;190:4405–19. [19] Most T. A natural neighbour-based moving least-squares approach for the element-free Galerkin method. Int. J. Numer Methods Eng. 2007;71:224–52. [20] Beissel S, Belytschko T. Nodal integration of the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 1996;139:49–74. [21] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin meshfree methods. Int. J. Numer. Methods Eng. 2001;50: 435–66. [22] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Comput. Mech. 2007;39:859–77. [23] Liu GR, Nguyen TT, Dai KY, Lam KY. Theoretical aspects of the smoothed finite element method (SFEM). Int. J. Numer. Methods Eng. 2007;71:902–30. [24] Timoshenko SP, Goodier JN. Theory of Elasticity.3rd edn. New York: McGrawHill; 1970.
Contents lists available at SciVerse ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Smoothed Galerkin methods using cell-wise strain smoothing technique X.Y. Cui, G.Y. Li n State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha, 410082, PR China
a r t i c l e i n f o
abstract
Article history: Received 14 October 2011 Accepted 21 December 2011 Available online 11 January 2012
A smoothed Galerkin method (SGM) using cell-wise strain smoothing operation is formulated in this paper. In present method, the field nodes can be divided into two types: boundary field nodes and interior field nodes. The background cells are divided into SC smoothing cells, and the strains in each smoothing cell are obtained using a gradient smoothing technique which can avoid evaluating derivatives of shape functions at integration point. The field variables of boundary points are approximated using linear interpolation of neighbour boundary field nodes, and the shape functions possess the Kronecker Delta property and facilitate the impositions of essential boundary conditions. The field variables of interior points are approximated using moving least-squares approximation using the support field nodes around them. A number of numerical examples are studied and confirm the significant features of the present methods: (1) can pass the standard patch test; (2) can easily impose essential boundary conditions as those in finite element method; (3) can avoid evaluating derivatives of shape functions; (4) no numerical parameter is required. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Meshfree Smoothed Galerkin method Gradient smoothing operation Cell-wise strain smoothing
1. Introduction In the past decades, various meshfree methods including smooth particle hydrodynamics [1], diffuse element method (DEM) [2], element free Galerkin (EFG) method [3], reproducing kernel particle method (RKPM) [4], finite point method (FPM) [5], H–P clouds [6], meshless local Petrov–Galerkin (MLPG) method [7], point interpolation method (PIM)[8], radial point interpolation method (RPIM) [9], among others, have been proposed and applied in more and more fields of particular engineering and scientific problems [10]. One of the best known meshfree methods is the EFG method, which uses the moving least-squares (MLS) approximation [11] for the interpolation of the displacement functions. The EFG method is a very promising method for the treatment of the partial different equations and has been applied and improved in a large variety of problems. However, there are some major problems in using EFG such as the imposition of essential boundary conditions, the complex computation in evaluating derivatives of meshfree shape functions, the uncertain numerical parameters like the nodal influence radius, the minimization of the computational cost and the implementation of an appropriate numerical integration scheme. Some works have been developed for the alleviation of the above-mentioned problems. The essential boundary conditions can not be imposed directly, because the MLS is not a true
n
Corresponding author. Tel.: þ86 731 88821717; fax: þ 86 731 88822051. E-mail address: [email protected] (G.Y. Li).
0955-7997/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.12.013
interpolation and the shape function constructed by it lacks the Kronecker delta property. Belytschko et al. [3] used the Lagrange multipliers to impose essential boundary conditions. Collocation methods [12] and penalty method [13] were also used to solve this problem. Other methods have been developed to overcome this problem including: the modification of the variational principle [14], the usage of the singular weighting functions [15], and the combination with finite elements [16]. Some coupling schemes were also investigated in the recent years. A coupled FE–EFG method was proposed by Belytschko et al. [17]. Gu and Liu proposed a coupled EFG–BEM for the analysis of two-dimensional solids [18], a natural neighbour-based MLS approach is developed coupled EFG and NEM method [19]. In this paper, the shape functions of field node local on the boundary of the problem domain are enforced to satisfy the Kronecker delta property, and the field variables of the point on the boundary are interpolated using the neighbour field nodes on the boundary, the essential boundary conditions can be imposed directly as those in FEM. Due to the complexity involved using Gauss integration in EFG methods, Beissel and Belytschko [20] proposed a nodal integration procedure of EFG. The unstabilized EFG with nodal integration displays spurious near-singular modes in some problems. To overcome this problem, Chen et al. [21] proposed a stabilized conforming nodal integration using a strain smoothing technique, and this method can reproduce a field exactly. Based on the strain smoothing technique, Liu et al. [22] proposed a smoothed finite element method (SFEM). The SFEM further divides the elements into some smoothing cells, computes the integrals along the edge of the smoothing cells, and has been proven to have good
826
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
properties. Liu et al. gave detailed theoretical aspects including stability, bound property and convergence about SFEM [23]. In present work, the idea of dividing the elements into some smoothing cells was further studied and used in the smoothed Galerkin method. In this paper, SGM is formulated using cell-wise strain smoothing operation in each background integration cells. In SGM, the field nodes are divided into two types: boundary field nodes and interior field nodes. The node shape functions of the boundary field nodes possess the Kronecker Delta property and the interior nodal variables are approximated using MLS approximation. The support nodes of each node for MLS approximation are selected based on the background cells. Depending on the requirement of the accuracy and stability, the background cells can be divided into several smoothing cells, and the strains in each smoothing cell are obtained using a gradient smoothing technique which can avoid evaluating derivatives of meshfree shape functions at integration point. The shape functions are directly created using linear interpolation for point on the boundary of the problem domain and MLS approximation for inner points. The support nodes selections of these points are also base on the background cells. For integration points in the background cells and integration points on the edges of the background cells, different support-node selections will be adopted to guard the continuity of the displacement field. The strains in the smoothing cells can be obtained as constants, hence there need not any other numerical integration in forming the stiffness matrix. To show the performance of the proposed method, a series of benchmark examples were presented, and excellent results have been obtained illustrating the efficiency and accuracy of the present methods. The paper is outlined as follows. Section 2 briefly introduces the basic equations of linear elasticity. Gradient smoothing method is introduced in Section 3. In Section 4, the shape function construction is presented in detail. Section 5 gives the presentation of the stiffness matrix and numerical examples are presented in Section 6. Conclusions and remarks are made in Section 7.
In Eq. (3), eT ¼ {exx, eyy, 2exy} is the vector of strains that relates to the displacements by the following compatibility equation
e ¼ Lu where u¼ {ux, uy}T is the displacement vector. Boundary conditions are given as follows: u ¼ uG on Gu
ð6Þ
nT r ¼ tG on Gt
ð7Þ
where uG is the specified displacement on the essential boundary
Gu, tG is the given traction on the natural boundary Gt, and n is the unit outward normal matrix expressed as 2 3 nx 0 6 7 n ¼ 4 0 ny 5 ny nx
O
@ @y
O
Gt
where de ¼L(du).
3. Gradient smoothing method The gradient smoothing method was proposed by Chen et al. [21], it can avoid evaluating derivatives of meshfree shape functions at point of interest and thus eliminates spurious modes. In present work, each background cell is further divided into some non-overlapping smoothing cells Ok ¼ [SC C ¼ 1 OkðCÞ , as shown in Fig. 1. The strain field in the Cth smoothing cell Ok(C) can be treated as a constant using following equation: Z Z ekðCÞ ¼ eðxÞfC ðxÞdO ¼ LuðxÞfC ðxÞdO ð10Þ Ok
Ok
ð1Þ
where fC ðxÞ is a given smoothing function that satisfies at least unity property Z fC ðxÞdO ¼ 1 ð11Þ
ð2Þ
A constant smoothing function is adopted as follows: ( 1=AC x A OkðCÞ fC ðxÞ ¼ 0 x= 2OkðCÞ
A solid mechanics problem of static elasticity can be described by equilibrium equation in the domain O which can be given by
where L is differential operator matrix defined as 2@ 3 0 @x 60 @ 7 L¼6 @y 7 4 5
ð8Þ
By applying principle of virtual work, a weak form is obtained for the equilibrium equation and the boundary conditions as follows: Z Z Z deT rdO duT bdO duT tG dG ¼ 0 ð9Þ
2. Brief of basic equation
LT r þ b ¼ 0
ð5Þ
@ @x
OkðCÞ
ð12Þ
where AC is the area of the smoothing cell Ok(C).
rT ¼{sxx, syy, sxy} is the stress vector and bT ¼{bx, by } is the body force vector. The stresses relate the strain via the constitutive equation as follows:
r ¼ De in which D is follows: 2 1 E 6n D¼ 4 1n2 0
ð3Þ the matrix of material constants that defined as
n 1 0
0
3
0 7 5
ð4Þ
1v 2
where E is Young’s modulus and n is Poisson’s ratio.
Fig. 1. The background cell Ok is divided into SC smoothing cells, Ok(C) is the Cth smoothing cell with the boundary Gk(C).
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
minimizing the weighted residual as follows:
Substituting Eq. (12) into Eq. (10), the smoothed strain field can be rewritten as Z Z 1 1 ekðCÞ ¼ LuðxÞdO ¼ nUuðxÞdG ð13Þ AC OkðCÞ AC GkðCÞ
JðxÞ ¼
NP X
NI ðxÞuI
ð14Þ
where NP is the number of the support nodes of point x, uI ¼{uI, vI}T is the SGM approximation for the displacement vector at node I, and NI(xi) is the SGM shape function introduced in the next section.
aðxÞ ¼ A1 ðxÞ BðxÞ u
ð19Þ
where u is the vector that collects the nodal parameters of displacement field in the support domain, and A(x) is called the weighted moment matrix defined by
4.1. Moving least-squares approximation
pj ðxÞaj ðxÞ ¼ pT ðxÞaðxÞ
n X
AðxÞmm ¼
The moving least-squares (MLS) approximation of a function u(x) denoted by uh(x) [3] is expressed as follows:
wðxxI ÞpðxI ÞpT ðxI Þ
BðxÞmn ¼ wðxx1 Þpðx1 Þ
ð15Þ
wðxx2 Þpðx2 Þ
wðxxn Þpðxn
ð21Þ
where pðxÞ ¼ 1
ð20Þ
I¼1
j¼1
ð17Þ
where r is the radius of the support of the weight function w(r). The stationary of functional J with respect to a(x) leads to:
4. Shape function construction for SGM
m X
wðr I Þ½pT ðxI ÞaðxÞuI 2
in which n is the number of nodes in the support domain of the interested point x, uI is the nodal value of u at node I, rI is the distance between two point x and xI, and w(r) is the weight function, for example, 8 2 3 > 2 > 4 rr þ 4 rr for rr r 12 > 3 > > < 2 3 ð18Þ wðrÞ ¼ 4 4 r þ 4 r 4 r for 12 o rr r 1 >3 3 r r r > > > > :0 for rr 41
I¼1
uh ðxÞ ¼
n X I¼1
where Gk(C) is the boundary of the smoothing cell Ok(C), n is the outward normal matrix containing the components of the outward normal vector on the boundary Gk. In the present method, the displacements of point x can be interpolated using the nodal displacements at the support nodes as follows: uh ðxÞ ¼
827
x
y
T
Substituting Eq. (19) into Eq. (15), the MLS approximation can be obtained
ð16Þ
uh ðxÞ ¼
is the basis vector in 2D, and a(x) is the corresponding coefficient vector. The coefficients a(x) are obtained at any point x by
n X
jI ðxÞuI ¼ uT ðxÞu
I¼1
Boundary
n1 SC3
xe 4
xb1
SC1
SC1
xi1 xb 2
n2
SC2
xi 3
SC2
SC4
SC3
xi 2
x e1
xe3
xe2 Boundary field node
SC1
Interior field node
Integration point in cell
Integration point on the edge of the cell
SC2 SC3
Integration point on the boundary
Fig. 2. Distribution of Integration points for a background divided into smoothing cells.
ð22Þ
828
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
where jI(x) is the MLS shape function defined by
jI ðxÞ ¼
m X
5. Discretized system equations
pj ðxÞ½A1 ðxÞBðxÞjI
ð23Þ
j¼1
4.2. SGM shape function In the present work, the problem domain is first discretized with background cells (triangular, quadrangular or n-side polygon), then further divided into SC smoothing cells. Though the smoothing operation, only shape functions of the integration points on the boundary of the smoothing cells will be evaluated. The shape functions of SGM are created using linear interpolation for points on the boundary of the problem domain, and using MLS approximation for the point in the problem domain. As shown in Fig. 2, for the integration point on the boundary of the problem domain, such as xb1and xb2, shape functions of them can be created by linear interpolation of nodes n1 and n2. For the integration point in the problem domain, such as xe1 xe4 and xi1 xi3, MLS shape functions are established, and the support-node selections are base on background cells. For the integration point in the background cells and the integration point on the edge of the background cell, different support-node selections will be adopted to guard the continuity of the displacement field. For the point of interest xi located in a cell i, the support node contains all nodes of the background cell which has common node of the cell i. As shown in Fig. 3, i1–i16 are support nodes for point xi. For the point of interest xe local on an edge (only for interior edge) k, support nodes contain all nodes of the background cells which have common node of that edge k. For example, e1–e12 are support nodes for point xe. For the point of interest xn is a field node (only for interior field node), support nodes contain all nodes of the background cells which contain the node n. For example, n1–n9 are support nodes for point xn. This scheme can usually select enough nodes and leads to less time consuming than the traditional scheme, in which circular support domain is used. This also scheme can be used to create MLS shape functions with high order of consistence and for extremely irregularly distributed nodes.
As shown in Fig. 4, Gk(C) is the boundary of smoothing cell Ok(C) divided into NG segments, and nJ is the outward normal matrix of the Jth segment. Substituting Eq. (14) into Eq. (13), the smoothed strain field can be rewritten by ! ! Z NG NP NP X 1 X 1 X ekðCÞ ¼ nUNI ðxÞdG uI ¼ nJ UNI ðxJ ÞlJ uI AC I ¼ 1 AC I ¼ 1 J ¼ 1 GkðCÞ I
¼ BkðCÞ uI
ð24Þ
where lJ is the length of Jth segment, xJ is the coordinate vector of I the gauss point on the Jth segment, and BkðCÞ is the smoothed strain displacement matrix corresponding to field node I. 2 3 " # nxJ 0 NG NG X X 0 1 1 I 6 0 n 7 NI ðxJ Þ yJ 5U BkðCÞ ¼ nJ UNI ðxJ ÞlJ ¼ l 4 0 NI ðxJ Þ J AC J ¼ 1 AC J ¼ 1 nyJ nxJ 2 3 N I ðxJ ÞnxJ 0 NG X 1 6 0 NI ðxJ ÞnyJ 7 ¼ ð25Þ 4 5lJ AC J ¼ 1 NI ðxJ ÞnyJ NI xJ nxJ where nJx and nJy are the components of the outward unit normal to the Jth boundary segment. One assumes that the material parameters E, u are same over a smoothing cell, the smoothed stress in smoothing cell Ok(C) can be easily obtained by
rkðCÞ ¼ DekðCÞ
ð26Þ
We now seek for a weak form solution of displacement field u that satisfies the following smoothed Galerkin weak form Z Z Z deT rdO duT bdO duT tG dG ¼ 0 ð27Þ O
O
Substituting Eqs. (14), (24) and (26) into Eq. (27), a set of discretized algebraic system equations can be obtained in the following matrix form: Kdf ¼ 0
ð28Þ
where f is the force vector defined as Z Z f ¼ UT ðxÞbdO UT ðxÞtG dG O
e
e
e
e
e
edge k e
x
e
e
i
i
i
i
e
e
i
i
i
i
e
cell i i
i
i
i
i
i
i
i
e
n
n
n
x
Gt
in which UT(x) is the shape function vector to build force vector. In Eq. (28), K is the (global) smoothed stiffness matrix assembled in the form of KIJ ¼
N cell X
KIJðkÞ
ð30Þ
k¼1
where the summation means an assembly process same as the practice in the FEM, Ncell is the number of the number of the background cells of the whole problem domain O, and KIJðkÞ is
node n
n
n
n2
n
l2 2
x n
n
n
n3 3 smoothing
Boundary field nodes
Fig. 3. Support-node selection in SGM for shape function construction. For interest point xi in a cell i, the support nodes are i1 i16. For interest point xe on an edge k, nodes e1 e12 are selected as the support nodes. For interest point xn as field node n, nodes n1 n9 are selected as the support nodes.
Γ
l1
domain Ω
l3 Interior field nodes
ð29Þ
Gt
1
l4
n1
4
n4
gauss point
Fig. 4. A smoothing cell Ok(C) with the boundary Gk(C), which is consist of four segments. lJ is the length of the Jth segment, nJ is the outward normal vector of the Jth segment.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
the stiffness matrix associated with Ok that is computed using KIJðkÞ ¼
SC Z X C¼1
OkðCÞ
I
J
ðBkðCÞ ÞT DðBkðCÞ ÞdO ¼
SC X
I
J
ðBkðCÞ ÞT DðBkðCÞ ÞAC
ð31Þ
C¼1
in which SC is number of the smoothing cells, AC and BkðCÞ are the area and smoothed curvature-deflection matrix of the Cth smoothing cell, respectively.
6. Numerical examples In this section, examples for 2D elasto-statics analyses are presented. For MLS approximation, linear basis vector (m¼ 3 in Eq. (15)) and weight function w given in Eq. (18) are used for constructing the shape functions. For all numerical examples, quadrilateral cells with one and four smoothing cells, triangular cells with one and three smoothing cells are adopted. For convenience, triangular cells with one smoothing cell for SGM will be named as SGM-TSC1; quadrilateral cells with four smoothing cells for SGM will be named as SGM-QSC4; etc.
6.1. Patch test The first numerical example is the standard displacement patch test; the mesh of a square parch is depicted in Fig. 5. A prescribed linear displacement field is imposed on the boundary Gu of the patch are computed using u ¼ 0:001ðx þyÞ v ¼ 0:001ðxyÞ
829
displacement is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uNnode uP exact unum ÞT ðuexact unum Þ u I I I u I ¼ 1 ðuI u Ed ¼ u NP node t ðuexact ÞT ðuexact Þ I I
ð33Þ
I¼1
where the superscript num denotes the displacement vector obtained using numerical methods, exact denotes the exact or analytical solution, and Nnode is the number of total field nodes. Table 1 lists the displacement norm errors of the numerical results for standard displacement patch test using both regularly and irregularly distributed field nodes shown in Fig. 5. All the numerical models can pass of the patch test to the machine accuracy and hence are capable to reproduce linear displacement field ‘‘exactly’’. 6.2. Cantilever beam problem A cantilever beam with length L and height D is studied as benchmark problem to test the convergence of the method, which is subjected to a parabolic traction at the free end as shown in Fig. 6. The beam is assumed to have a unit thickness so that plane stress condition is valid. The analytical solution from Timoshenko and Goodier [24] can be given by " !# 2 Py 2 D ux ðx,yÞ ¼ ð6L3xÞx þ ð2þ nÞ y 6EI 4
Table 1 Displacement error norm of numerical results for the standard patch test.
ð32Þ
The material properties of patch are E¼1.0 and n ¼0.25. To satisfy the patch test, the linear displacement field approximated exactly, and the numerical solution at any interior nodes should be in exact agreement the analytic ones given in Eq. (32). To examine the numerical error precisely, an error norm in
SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Regularly distributed nodes
Irregularly distributed nodes
1.48431617E 14 9.81195622E 15 1.79482393E 14 1.17421298E 14
1.33104257E 14 1.49256754E 14 1.00570811E 14 1.08551760E 14
Fig. 5. Displacement patch test: (a), (b) and (c) uniform grid; (d), (e) and (f) irregular grid.
830
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
" # P D2 x 2 2 3ny ðLxÞ þ ð4 þ5nÞ þð3LxÞx uy ðx,yÞ ¼ 6EI 4
ð34Þ
D2 y2 4
!
D x
L Fig. 6. Cantilever beam subjected to a parabolic traction on the right tip.
-1
-0.5 SGM-TSC1 SGM-TSC3 FEM, T3
-1
SGM-QSC1 SGM-QSC4 FEM, Q4
-1.5
-2
Log10(Ed)
Log10(Ed)
-1.5
-2
-2.5
-3
-3.5 -2.5 -4
-3 -0.4
-0.2
0 0.2 Log10(Δh)
0.4
-4.5 -0.4
0.6
-0.2
0 0.2 Log10(Δh)
0.4
0.6
Fig. 7. Convergence of the numerical results in displacement error norm for the problem of cantilever beam solved using different method.
-3
-3.4 SGM-QSC1
SGM-TSC1 -3.2
-3.6
SGM-TSC3
SGM-QSC4 FEM, Q4
FEM, T3
-3.8
-3.4
-4 Log10(Ee)
Log10(Ee)
-3.6 -3.8
-4.2 -4.4
-4 -4.6 -4.2
-4.8
-4.4 -4.6 -0.4
-5
-0.2
0 0.2 Log10(Δh)
0.4
ð35Þ
where I¼D3/12 is the moment of inertia for the beam. The related parameters for the problem are: E¼3.0 107 Pa, n ¼0.3, L¼48 m, D ¼12 m and P¼ 1000 N. Using the same set of distributed field nodes, the cantilever beam is studied using the proposed methods with triangular and quadrilateral background cells, respectively. Further, this problem is also studied using FEM triangular element (T3) and quadrilateral (Q4) with the same distributed field nodes. In the numerical computations, the displacement boundary on the left uses the exact displacements obtained from Eq. (34) and the loading boundary uses the distributed shear stresses in Eq. (35).
y
P
Py P sxx ðx,yÞ ¼ ðLxÞ, syy ðx,yÞ ¼ 0, sxy ðx,yÞ ¼ I 2I
0.6
-5.2 -0.4
-0.2
0 0.2 Log10(Δh)
0.4
0.6
Fig. 8. Convergence of the numerical results in energy error norm for the problem of cantilever beam solved using different method.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
Fig. 7 plots the convergence of the solutions in displacement error norm for the cantilever beam solved using different methods. The mesh size parameter Dh is taken to be the average nodal spacing. For the left figure of Fig. 7, all methods use the triangular cells. It can be found that SGM-TSC1 and SGM-TSC3 are more accurate than FEM (T3). For the right figure of Fig. 7, quadrilateral cells are used for presented methods. Ones find that SGM-QSC4 gives better displacement results than FEM (Q4). SGM-QSC4 not only gives the best result, but also has the highest rate of convergence 2.43, which is much larger than the theoretical value 2.0. Fig. 8 shows the convergence of solutions in energy norm for the cantilever beam solved using different methods. The energy
error norm is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uR ðeexact enum ÞT Dðeexact enum ÞdO 1u Ee ¼ t O R exact ÞT Dðeexact ÞdO A O ðe
ð36Þ
where A is the area of the problem domain. It is found that SGMTSC1 and SGM-TSC3 have better accuracy and converge faster compared to the FEM (T3). SGM-TSC1 provides good accuracy in energy and a high convergence rate (1.193). Form the right figure of Fig. 8, it can be found that all of the SGM have better accuracy and higher rate of convergence than FEM (Q4). It is observed that SGM-QSC1 presents superconvergent property and the
5
4.8 SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
4.8
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
4.75 4.7
4.6
4.65 Strain energy
Strain energy
831
4.4
4.2
4.6 4.55 4.5 4.45
4
4.4 3.8 4.35 3.6 0
1000
2000 3000 DOF
4000
4.3
5000
0
1000
2000 3000 DOF
4000
5000
Fig. 9. Solutions (strain energy) converging to the exact solution for the problem of cantilever beam solved using different method.
x 10-3
-8.6
Vertical displacement of middle tip node
-7.6 -7.8 -8 -8.2 -8.4 -8.6 -8.8 -9 -9.2
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
-8.7 Vertical displacement of middle tip node
SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
-7.4
x 10-3
-8.8 -8.9 -9 -9.1 -9.2 -9.3 -9.4
-9.4 0
1000
2000 3000 DOF
4000
5000
-9.5
0
1000
2000 3000 DOF
4000
5000
Fig. 10. Solutions (vertical displacement of middle tip node) converging to the exact solution for the problem of cantilever beam solved using different method.
832
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
convergence rates of them are 1.65, which are much higher than theoretical value 1.0. To study the convergence property, the strain energy of the cantilever beam is computed using different methods and plotted
1000 800 600
Analytical solution SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Normal stress σx
400 200 0 -200 -400 -600 -800 -1000 -6
-4
-2
0 y
2
4
6
Fig. 11. Normal stress sx distribution along the cross-section of the cantilever at x¼ L/2 obtained using different methods and the same distributed 637 nodes.
Analytical solution SGM-TSC1 SGM-TSC3 SGM-QSC1 SGM-QSC4
Shear stress σxy
-40 -60 -80 -100 -120 -140 -6
-4
-2
0 y
6.3. Infinite plate with a circular hole An infinite plate with a circular hole (a) subjected to a unidirectional tensile load (p) in the x direct is studied. Owing to the symmetry, only one quarter is model as shown in Fig. 13. Symmetry conditions are imposed on the left and bottom edges, and the inner boundary of the hole is traction free. The analytical solution for the stress can be given in [24]
a2 3 3a4 cos 2y þcos 4y þ 4 cos 4y sxx ¼ p 1 2 2 r 2r 2
a 1 3a4 cos 2ycos 4y þ 4 cos 4y syy ¼ p 2 2 r 2r 2
a 1 3a4 sin 2ysin 4y 4 sin 4y sxy ¼ p 2 ð37Þ 2 r 2r
0 -20
in Fig.9. In the left figure, it is easy to see that SGM-TSC1 has the upper bound solutions and SGM-TSC3 has lower bound solution. From the results, ones find that the numerical model is softer than the exact model when 1 smoothing cell is used for SGM, the strain energy of SGM-TSC1 is a litter larger than the exact solution. The model of SGM-TSC3 is stiffer than the exact model but softer than FEM (T3), and it gives more accurate results than FEM (T3) but still a litter lower than exact solution. Form the right figure of Fig. 9, it can be found that all methods present lower and upper bound property. For the present SGM, the numerical model becomes stiff when the number of smoothing cell increases. SGM-QSC4 presents best accurate than other schemes plotted in the figure. It is very close to exact solution and has the lower bound property. The vertical displacements of middle tip node of the cantilever beam computed using different methods are plotted in Fig. 10. The results again show that SGM-TSC1 and SGM-TSC3 can give better results than FEM (T3) with the same mesh. We find that the convergence trend of displacements is same as it of strain energy. SGM-TSC4 gives the best results than other schemes plotted in the figure. Figs. 11 and 12 illustrate the distribution of the normal stress sx and shear stress sxy on the cross-section x¼ L/2 of the beam. It can be found that the results obtained using present methods agree well with the analytical solutions. SGM-QSC4 exhibits best results compared with other presented numerical methods.
2
4
6
Fig. 12. Shear stress sxy distribution along the cross-section of the cantilever at x¼ L/2 obtained using different methods and the same distributed 637 nodes.
where r and y are the polar coordinates and y is measured counterclockwise from the positive x-axis. The displacement fields can be calculated as follows:
a4 p kþ1 a2 ur ¼ r þcos 2y þ 1 þðk þ1Þcos 2y 3 cos 2y 4G 2 r r 2 4 p a a rð1kÞ þ 3 sin 2y ð38Þ uy ¼ 4G r r
y y
p
p
b
r a
θ
b
x
x Fig. 13. Infinite plate with a circular hole subjected to uniform tensile and its quadrant model.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
833
Fig. 14. Domain discretization for the infinite plate with a hole (9 9 nodes).
-0.5
-1
-1.5
SGM-TSC1 SGM-TSC3 FEM, T3
-1
SGM-QSC1 SGM-QSC4 FEM, Q4
-1.5
Log10(Ed)
Log10(Ed)
-2
-2.5
-2
-2.5
-3 -3 -3.5
-4 -1.2
-3.5
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
-4 -1.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
Fig. 15. Convergence of the numerical results in displacement error norm for the problem of infinite plate with a circular hole solved using different method.
where G¼
E , 2ð1þ nÞ
(
k¼
34n n
1n
for plane stress for plane strain
ð39Þ
E is Young’s modulus and n is Poisson’s ratio. We studied the problem under plane strain conditions and traction boundary conditions are imposed on the upper and right edges with the analytical stresses obtained using Eq. (37). The parameters for this problem are: a¼1 m, b¼5 m, p ¼1 Pa, E¼1 103 Pa and n ¼0.3. The 9 9, 17 17, 25 25, 33 33 and 49 49 distributed nodes are used for computation and a 9 9 node distribution is shown in Fig. 14. Fig. 15 plots the convergence of the solutions in displacement error norm for the problem of infinite plate with a circular hole solved using different methods. Form the left figure of Fig. 15, we find that SGM-TSC3 exhibits the best accuracy compared with other methods using the triangular cells and SGM-TSC1 has the highest convergence rate (2.135). For the methods using quadrilateral background cells, the solutions of all the methods converge well and exhibit different accuracy and convergence rate. In terms of accuracy, SGM-QSC1 is found worse than the FEM (Q4), and all other schemes are at least better than the FEM (Q4). SGM-QSC4 still gives good accuracy as the previous case. It must be pointed that the convergence rate of SGM-Q1 is close to 3.264, but the accuracy of it is not good enough.
Fig. 16 plots the convergence of the solutions in energy error norm for the problem of infinite plate with a circular hole solved using different methods. The methods using both the triangular and quadrilateral background cells are found converging well. SGM-TSC1 again gives the highest convergence rate (1.211) in energy norm as in displacement norm, and it is higher than the theoretical value 1.0. Form the right figure of Fig. 16, it can be found that all of the SGM have better accuracy and higher rate of convergence than FEM (Q4). SGM-QSC1 still present superconvergent property and the convergence rates of them are 2.824, which are more than twice higher than the theoretical value 1.0. The convergence trends of strain energy of the infinite plate with a circular hole computed using different methods are plotted in Fig. 17. SGM-TSC1 still has the upper bound solutions and SGM-TSC3 has lower bound solution. For quadrilateral background cells, SGM-QSC4 presents lower bound property, and SGM-QSC1 presents upper bound property.
7. Conclusion In this work, SGM is proposed using cell-wise strain smoothing operation. The background cells are divided into several smoothing cells and the strain fields in them are projected to a set of constant fields. The strain smoothing stabilization technique is
834
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
-3.6
-2.5 SGM-TSC1 SGM-TSC3 FEM, T3
-3.8
SGM-QSC1 SGM-QSC4 FEM, Q4
-3
-4
-4.2
Log10(Ee)
Log10(Ee)
-3.5
-4.4
-4
-4.5 -4.6 -5
-4.8
-5 -1.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-5.5 -1.2
-0.2
-1
-0.8 -0.6 Log10(Δh)
-0.4
-0.2
0.01183
0.01185
0.01182
0.01184
0.01181
0.01183
0.01180
0.01182
0.01179
0.01181
Strain energy
Strain energy
Fig. 16. Convergence of the numerical results in energy error norm for the problem of infinite plate with a circular hole solved using different method.
0.01178 0.01177 0.01176
0.01179 0.01178
0.01175
SGM-TSC1 SGM-TSC3 FEM, T3 Analytical
0.01174 0.01173
0.01180
0
1000
2000 3000 DOF
4000
5000
0.01177
SGM-QSC1 SGM-QSC4 FEM, Q4 Analytical
0.01176 0.01175
0
1000
2000 3000 DOF
4000
5000
Fig. 17. Solutions (strain energy) converging to the exact solution for the problem of infinite plate with a circular hole solved using different method.
employed to restore the conformability and to improve the accuracy and the rate of convergence. Only line integrations along the smoothing cells are needed and no derivative of the shape functions is involved in computing the field gradients to form the stiffness matrix. The support-node-selection for the MLS approximation is based on the background cells and no numerical parameter like the nodal influence radius is required. The results of investigated numerical examples illustrate the effectivity of the present method. Besides the good accuracy, stabilization and high convergence rate, the following features of the present method are noteworthy: (1) can pass the standard patch test; (2) can directly impose essential boundary conditions as those in finite element method; (3) can avoid evaluating derivatives of meshfree shape functions; (4) no uncertain numerical parameter like the nodal influence radius is required.
Acknowledgments This work was supported by National Science Foundation of China (11002053), National Basic Research Program of China (2010CB328005), Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body (51075005). The authors also give sincere thanks to Prof. Liu G.R. at University of Cincinnati for his help. References [1] Monaghan JJ. An Introduction of SPH. Comput. Phys. Commun. 1982;48:89–96. [2] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput. Mech. 1992;10(5):307–18.
X.Y. Cui, G.Y. Li / Engineering Analysis with Boundary Elements 36 (2012) 825–835
[3] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int. J. Numer. Methods Eng. 1994;37(2):229–56. [4] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. Int. . Numer. Methods Eng. 1995;20:1081–106. [5] Onate Idelsohn ES, Zienkiewicz OC, Tylor RL. A finite point method in computational mechanics, applications to convective transport and fluid flow. Int. J. Numer. Methods Fluids 1996;39:3839–66. [6] Duarte CA, Oden JT. An H–P adaptive methods using clouds. Comput. Methods Appl. Mech. Eng. 1996;139:237–62. [7] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput. Mech. 1998;22(2):117–27. [8] Liu GR, Gu YT. A point interpolation method for two-dimensional solids. Int. Numer. Methods Eng. 2001;50:937–51. [9] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int. J. Numer. Methods Eng. 2002;54:1623–48. [10] Liu GR. Meshfree methods: Moving Beyond the Finite Element Method. Boca Raton, FL: CRC Press; 2002. [11] Lancaster P, Salkauskas K. Surface generated by moving least squares methods. Math. Comput 1981;37:141–58. [12] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin method. Int. . Numer. Methods . Eng. 1996;39:141–58. [13] Zhu T, Atluri N. Modified collocation method and a penalty function for enforcing essential boundary conditions in the element free Galerkin method. Comput. Mech. 1998;21:211–22.
835
[14] Lu YY, Belytschko T, Gu L. A new implementation of the element-free Galerkin method. Comput. Methods in Appl. Mech. Eng. 1994;113:397–414. [15] Kaljevic I, Saigal S. An improved element free Galerkin formulation. Int. J. Numer. Methods Eng. 1997;40:2953–74. [16] Krongauz Y, Belytschko T. Enforcement of essential boundary conditions in meshless approximations using finite elements. Comput. Methods Appl Mech Eng. 1996;131:133–45. [17] Belytschko T, Organ D, Krongauz Y. A coupled finite element-element-free Galerkin method. Comput. Mech. 1995;17:186–95. [18] Gu YT, Liu GR. A coupled element free Galerkin/boundary element method for stress analysis of two-dimensional solids. Comput. Methods Appl. Mech. Eng. 2001;190:4405–19. [19] Most T. A natural neighbour-based moving least-squares approach for the element-free Galerkin method. Int. J. Numer Methods Eng. 2007;71:224–52. [20] Beissel S, Belytschko T. Nodal integration of the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 1996;139:49–74. [21] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin meshfree methods. Int. J. Numer. Methods Eng. 2001;50: 435–66. [22] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Comput. Mech. 2007;39:859–77. [23] Liu GR, Nguyen TT, Dai KY, Lam KY. Theoretical aspects of the smoothed finite element method (SFEM). Int. J. Numer. Methods Eng. 2007;71:902–30. [24] Timoshenko SP, Goodier JN. Theory of Elasticity.3rd edn. New York: McGrawHill; 1970.