Computational Mechanics 30 (2003) 196–211
Springer-Verlag 2003 DOI 10.1007/s00466-002-0379-y

The least-squares meshfree method for solving linear elastic problems K.-C. Kwon, S.-H. Park, B.-N. Jiang, S.-K. Youn

196 Abstract A meshfree method based on the first-order least-squares formulation for linear elasticity is presented. In the authors’ previous work, the least-squares meshfree method has been shown to be highly robust to integration errors with the numerical examples of Poisson equation. In the present work, conventional formulation and compatibility-imposed formulation for linear elastic problems are studied on the convergence behavior of the solution and the robustness to the inaccurate integration using simply constructed background cells. In the least-squares formulation, both primal and dual variables can be approximated by the same function space. This can lead to higher rate of convergence for dual variables than Galerkin formulation. In general, the incompressible locking can be alleviated using mixed formulations. However, in meshfree framework these approaches involve an additional use of background grids to implement lower approximation space for dual variables. This difficulty is avoided in the present method and numerical examples show the uniform convergence performance in the incompressible limit. Therefore the present method has little burden of the requirement of background cells for the purposes of integration and alleviating the incompressible locking. Keywords Least-squares, Meshfree method, Integration error, Incompressible locking, Linear elasticity

1 Introduction Meshfree methods have been focused on with the increasing demands of large-scale numerical simulations in diverse fields for which frequent reconstruction of analysis model is required (Nayroles et al., 1992; Belytschko et al., 1994; Liu et al., 1995; Duarte and Oden, 1995; Melenk and Babusˆka, 1996). While no mesh is required in the construction of shape functions, most meshfree methods use background cells with Gaussian quadrature for integration

Received: 16 December 2001 / Accepted: 4 November 2002

K.-C. Kwon, S.-H. Park, S.-K. Youn (&) Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, 373-1, Gusung-dong, Yusung-ku, Taejon, 305-701, Korea e-mail: [email protected] B.-N. Jiang Department of Mathematics and Statistics, Oakland University, Rochester, Michigan, USA

and the construction of background cells for accurate integration is also burdensome and can make meshfree methods less effective. Therefore several attempts have been made to avoid the use of integration cells. Beissel and Belytschko (1996) employed a nodal integration instead of using background cells, showed that spurious singular modes can be induced by the nodal integration as in SPH, and suggested a stabilization scheme. Nagashima (1999) suggested a stabilization scheme for nodal integration based on the Taylor series expansion and also proposed a bucket algorithm for computing nodal volumes. On˜ate et al. (1996) used the collocation method instead of Galerkin method. The collocation method, however, can suffer from the stability problem as in the nodal integration and also it requires higher-order derivatives of shape functions and results in non-symmetric stiffness matrix. Breitkopf et al. (2000) introduced evaluation points in addition to nodes and employed a two-stage approach in order to avoid the computation of the second-order derivatives of shape functions in the collocation method. The evaluation points can also serve as stabilization terms. Chen et al. (2001) suggested a stabilized conforming nodal integration scheme. They pointed out that the inaccuracy in solution by the nodal integration is due to the disobedience of the integral identity given by the divergence theorem. In their scheme, linear fields are exactly integrated, where the integration is performed at the vertices of the Voronoi cells constructed from the nodes. Atluri and Zhu (1998) suggested the local Petrov–Galerkin approach in which integration is performed locally. Recently Park and Youn (2001) proposed a meshfree method based on leastsquares formulations, in which integral identities are not involved. Through numerical examples of Poisson equation it is also shown in the same paper that the leastsquares meshfree method is robust to integration errors. The difficulty of the incompressible locking can be resolved by various methods in the finite element method. In meshfree methods, the volumetric locking can be alleviated using large domains of influence (Belytschko et al., 1994). However, for small domains of influence, the direct application of displacement-based Galerkin formulation exhibits volumetric locking (Dolbow and Belytschko, 1999; Chen et al., 2000a, b; De and Bathe, 2001). In many practical problems, small domains of influence are preferred for improving the local resolution and computational efficiency. Therefore some meshfree methods based on mixed formulations have been proposed to alleviate the incompressible locking. However, the implementation of

mixed formulations in meshfree framework is not straightforward. In general, mesh (or cells) is involved in the implementation of lower order approximation for the dual variables to satisfy the so-called LBB condition. This can diminish the effectiveness of meshfree methods. Based on a mixed formulation, Dolbow and Belytschko (1999) proposed a selective reduced integration procedure. They used background cells with full quadrature to integrate the deviatoric part of the weak form and nodal integration strategy by Voronoi cells is used for the volumetric part. As a similar approach the pressure projection technique for meshfree method (Chen et al., 2000a, b) has been developed, in which reduced quadrature, usually one-point quadrature, in a background cell is used for the volumetric part. De and Bathe (2001) have presented a numerical inf-sup test to analyze the stability and optimality of several mixed discretization schemes using the meshfree approximations. From the mathematical point of view, the first-order least-squares formulation can be regarded as another type of mixed formulation. However it differs from the mixed formulations based on symmetric weak formulations in that the shape functions of same order can be employed for both primal and dual variables. This leads to more accurate computation of dual variables. Moreover the system matrix is always positive-definite. Therefore the uniform convergence performance in the incompressible limit can be achieved (Cai et al., 1997). This is a very important feature in meshfree methodology since the additional use of cells can be avoided. For more details of the first-order least-squares formulations for linear elasticity and their regularity requirements for the optimal convergence, one can refer to Cai et al. (1997, 1998, 2000a, b). The aim of the present work is to propose a truly meshfree method based on first-order least-squares formulations for linear elasticity, i.e. with little requirement of carefully constructed background cells. It is accomplished by the outstanding features of the first-order leastsquares formulation, the robustness to integration errors and higher order approximation for dual variables without the incompressible locking. For this purpose, first-order least-squares meshfree methods employing the conventional and the compatibility-imposed formulations are proposed. In the presented method, a simple integration algorithm can be used and no additional treatment for the incompressible locking is required. In the following sections, after the moving least-squares approximation and general first-order least-squares formulation are briefly reviewed, the conventional formulation and compatibilityimposed formulation are presented and then numerical examples are solved for the comparison of least-squares meshfree method and Galerkin meshfree method. Concluding remarks are given in the last section.

more details, please refer to other literatures (Lancaster and Salkauskas, 1981; Duarte and Oden, 1995).
, Let X be an open domain of Rd , d = 1, 2, or 3, xI 2 X n I ¼ 1; 2; . . . ; N, the given nodal points, fpi gi¼1 ; n  N, some basis for the approximation and wI ðxÞ non-negative weight function associated with the nodal point xI . Then the MLS shape function for the I-th node wI can be written as

wI ðxÞ ¼ pT ðxÞA1 ðxÞBI ðxÞ

ð1Þ

where

197

pT ðxÞ ¼ ½p1 ðxÞ; p2 ðxÞ; . . . ; pn ðxÞ

ð2aÞ

AðxÞ ¼ PT WðxÞP

ð2bÞ

BI ðxÞ ¼ wI ðxÞpðxI Þ 2 p1 ðx1 Þ p2 ðx1 Þ 6 p1 ðx2 Þ p2 ðx2 Þ 6 P¼6 . .. 4 .. .

  .. .

pn ðx1 Þ p1 ðx1 Þ 7 7 .. 7 . 5

p1 ðxN Þ p2 ðxN Þ    pn ðxN Þ 2 w1 ðxÞ 0  0 6 0 w ðxÞ    0 2 6 WðxÞ ¼ 6 . .. .. .. 4 .. . . . 0

0

ð2cÞ

3

ð2dÞ

3 7 7 7 5

ð2eÞ

   wN ðxÞ

Derivatives of the MLS shape functions are then given as follows:

wI ;i ¼ pT ;i A1 BI þ pT A1 BI ;i pT A1 A;i A1 BI

ð3Þ

In computation of the MLS shape functions, no explicit mesh, i.e., element connectivity, is required. This fact enables the meshfree strategies for the solution of differential equations. Throughout the present work, the linear basis, pT ¼ ½1; x; y , and the following weight function are used, which is known to be integrated easily (Duarte and Oden, 1995).

pffiffiffiffiffiffiffiffi 4=pð1  kx  xI k2 =h2I Þ4 wI ðxÞ ¼ 0

if kx  xI k < hI if kx  xI k  hI ð4Þ

where hI is the influence radius of the nodal point xI . If not stated otherwise, hI has been selected so that hI ¼ 1:5d, where d is the nodal spacing.

3 First-order least-squares formulation The Galerkin method, which is based on the weighted residual form, has been successfully applied to various 2 fields of engineering and applied sciences. However, Moving least-squares approximation serious difficulties such as oscillation, instabilities of the The moving least-squares method (Lancaster and Salkauskas, 1981) is the most widely used approximation solution and inaccuracy of its derivatives have been scheme in recent meshfree methods. Also in the present encountered during the attempts to apply the Galerkin work the MLS approximation is employed. In the follow- finite element method to the non-self-adjoint equations in ing, the construction of the MLS shape functions and their the fields of fluid dynamics and electro-magnetics. The least-squares finite element method (LSFEM) (Jiang, 1998) derivatives is briefly reviewed with their properties. For

198

has been used to overcome these difficulties by a simple idea of minimizing the squared residuals. The Galerkin formulation based meshfree methods are extremely sensitive to the integration errors. Both the Galerkin formulation and the least-squares formulation are based on the concept of averaging by integration. One of the essential ingredients in the Galerkin formulation, is the integration by parts and the application of the divergence theorem. This procedure assumes exact integration and thus inaccuracy in the integration is directly reflected to the solution accuracy. In the least-squares formulation, these integral identities are not employed and integration may be regarded only as a tool of averaging. This implies that simply generated integration points, which do not provide accurate integration, can be used without seriously degrading the solution accuracy in the least-squares formulation. Therefore the least-squares formulation can be employed effectively in devising a truly meshfree method (Park and Youn, 2001; Park et al., 2002). In the present section, a general first-order least-squares formulation based on L2 -norm is reviewed with convergence characteristics. For details and the improved schemes for various problems, please refer to the work of Jiang (1998) and the references therein.

3.1 General first-order least-squares formulation All the higher-order differential equations can be transformed to the first order differential equations. Many least-squares formulations are based on the first-order differential equations. One of the main reasons for using the first-order equations is to make use of C0 -elements. Let X be an open domain of Rd ; d ¼ 1, 2 or 3 and C its boundary. Let us suppose the linear first-order boundary value problem. Au ¼ f

in X

ð5aÞ

Bu ¼ g on C

ð5bÞ

where A is a first-order linear differential operator:

Au ¼ A0 u þ

d X i¼1

Ai

ou oxi

ð6Þ

IðvÞ ¼ kRðvÞk20 ¼ kAv  fk20 ¼ ðAv  f; Av  fÞ0

ð9Þ

2

where the inner product ð; Þ0 in L ðXÞ is defined by

ðu;vÞ0 ¼

Z

uT v dX

ð10Þ

X

The least-squares method is to minimize the quadratic functional IðvÞ. Applying the stationary condition on IðvÞ, the following variational formulation is obtained: Find u 2 V such that

Bðu;vÞ ¼ FðvÞ

8v2V

ð11Þ

where

Bðu;vÞ ¼ ðAu; AvÞ0

ð12aÞ

FðvÞ ¼ ðf; AvÞ0

ð12bÞ

It is noted that the bilinear form (12a) is symmetric and only essential boundary conditions are present in the first-order least-squares formulation. In the present work, penalty method is used to enforce the boundary conditions. For a well-posed problem, the operator A is bounded below. Thus the bilinear form (12a) leads to a symmetric positive-definite matrix.

3.2 Convergence of least-squares formulation The convergence characteristics of the least-squares method is somewhat different according to the satisfaction of the H1 -elliptic condition of the differential equation (5). The following two theorems (Jiang, 1998; Park et al., 2002) present general guidelines on the convergence of the first-order least-squares formulation. Theorem 1. If the linear first-order problem (5) is bounded below by L2 -norm (i.e., kAukL2  akukL2 8 u 2 V) and its solution u is sufficiently smooth, then for its V-conforming approximation by piecewise polynomials of degree r on quasi-uniform meshes with mesh size h, there exists a constant C > 0 independent of h and u such that the least-squares approximated solution uh satisfies

ku  uh kL2 ðXÞ  Chr jujHrþ1 ðXÞ

ð13Þ

T

in which u ¼ ½u1 ; u2 ; . . . ; um is a vector of m unknown functions of xT ¼ ½x1 ; x2 ; . . . ; xd , Ai the Neq  m coefficient matrices which characterize the differential operator A, f a given function in the domain, B a boundary algebraic operator, and g a given function on the boundary. It should be mentioned that the number of equations Neq must be greater or equal to the number of unknowns m. Suppose that f 2 L2 ðXÞ. Let V be a subspace of L2 ðXÞ such that the boundary condition (5b) is satisfied:

Theorem 2. If the linear first-order problem (5) is H1 elliptic (i.e., kAukL2  akukH1 8 u 2 V) and its solution u is sufficiently smooth, then for its V-conforming approximation by piecewise polynomials of degree r on quasiuniform meshes with mesh size h, there exist constants C1 > 0 and C2 > 0 independent of h and u such that the least-squares approximated solution uh satisfies

Bv ¼ g on C 8 v 2 V

ku  uh kL2 ðXÞ  C1 hrþ1 jujHrþ1 ðXÞ

ð7Þ

ð14aÞ

It is noted that the first-order differential operator A maps ku  uh k 1  C2 hr juj rþ1 ð14bÞ H ðXÞ H ðXÞ the subspace V into L2 ðXÞ. For an arbitrary trial function v 2 V, the residual is defined at each point in X as follows: It should be noted that the dual variables (e.g., heat flux, stress, etc.) are included in the solution vector. Therefore RðvÞ ¼ Av  f in X ð8Þ the L2 -norms in (13) and (14a) are similar to the H1 -norm Then the integral of the squared residual is given in the in the displacement-based Galerkin formulation. Also it following quadratic form: can be roughly said that one-order higher regularity of the

solution is required to achieve the above convergence behavior than in the Galerkin formulation. For H1 -elliptic (coercive) problems, which are usually obtained by imposing compatibility condition, the rate of convergence is optimal, i.e., it is as high as possible for a given order of approximation. Several least-squares formulations have been developed for various kinds of differential equations so that the elliptic conditions are satisfied (Jiang, 1998). Although the rate of convergence is sub-optimal for general problems, the dual variables can be computed with at least the same rate of convergence as the Galerkin method. Even if the mixed Galerkin formulation is used, the LBB condition forces that lower-order approximants are used for the dual variables. Using leastsquares formulation, however, the LBB condition does not need to be satisfied and thus the accuracy of dual variables can be improved by changing the form of the first-order differential equations so as to satisfy the elliptic condition.

3.3 Discrete least-squares formulation Formally the least-squares formulation presented above can be stated as follows: Z Minimize ðAv  fÞT ðAv  fÞdX ð15Þ X

Using Gaussian quadrature for integration, this can be rewritten in the following discrete form.

Minimize

X ðAvI  f I ÞT ðAvI  f I ÞWI

ð16Þ

I

where WI are the weight factors at the integration points (including volume information). Considering the leastsquares curve fitting, where usually integration is not involved but summation is performed on at a given set of finite points, it can be presumed that the accurate integration is not an essential ingredient in the least-squares formulation and the discrete least-squares formulation does not deteriorate the solution with inaccurate integration. This feature makes it possible to use an inaccurate and simple integration scheme in the least-squares formulation. In other words, it may be sufficient to use in (16) some simply chosen set of integration points instead of Gaussian quadrature points which requires the construction of background cells. Although integration accuracy can be degraded with simply chosen integration points, if integration errors are in a moderate range, the solution accuracy may not be deteriorated greatly in the least-squares formulation. By this characteristics, the least-squares meshfree method can be a truly meshfree method, and this feature has been successfully shown for 2D Poisson equation (Park and Youn, 2001).

4 First-order least-squares formulation for linear elasticity In this section, the conventional formulation is introduced and the compatibility-imposed formulation is proposed for the linear elasto-static problem. For the first-order least-squares formulation second-order differential

equations of linear elastic problem are transformed into the system of first-order differential equations by introducing dual variables. It should be mentioned that other first-order least-squares formulations could be also obtained depending on what are chosen as dual variables. In the conventional formulation displacement components and stress components are used as primal and dual variables, respectively. The compatibility-imposed formulation employs other dual variables to enforce the compatibility condition. For its H1 -like ellipticity, certain assumptions on the solution’s regularity should be met. The required regularity can be mitigated employing H1 norm instead of L2 -norm. For other least-squares formulations and more details of the regularity issues, one can refer to the work of Cai et al. (1997, 1998, 2000a). If the solution is sufficiently smooth, the compatibility-imposed formulation gives the optimal rate of convergence for dual variables as shown in the next section. As reviewed in the previous section, the first-order least-squares formulation can be regarded as a sort of mixed formulation and the system matrix is always symmetric and positive-definite. Thus the volumetric locking can be avoided (Cai et al., 1997) and the convergence characteristics maintained in the nearly incompressible condition by introducing appropriate dual variables. In the presented formulations these features are achieved since no singularity is caused in the incompressible limit.

4.1 Conventional formulation Before going directly into the least-squares formulations, the governing equations of the linear elasto-static problem on the domain X bounded by the boundary C, which are second-order differential equations, are summarized. rij ;j þbj ¼ 0 in X eij ¼

1 2 ðui ;j

þuj ;i Þ

ð17aÞ in X

ð17bÞ

1þm m ð17cÞ rij  raa dij in X E E where rij and eij are the stress tensor and strain tensor, respectively, which correspond to the displacement field ui ; bi is the body force; E and m are Young’s modulus and Poisson ratio, respectively; dij is the Kronecker delta; and ðÞ;i denotes oðÞ=oxi . The governing equations are composed of the equation of equilibrium (17a), the strain– displacement relationship (17b) and the constitutive equation (17c). The corresponding boundary conditions are given as follows. eij ¼

ui ¼ u
i rij nj ¼
ti

on Cui on Cti

ð18aÞ ð18bÞ

where u
i and
ti are the prescribed displacements and tractions, respectively, on the displacement boundary Cui and on the traction boundary Cti , and ni is the unit outward normal to the boundary C. In the conventional formulation, the governing equations are transformed into first-order differential equations by introducing displacement components ui , as

199

primal variables and dimensionless stress components, sij ¼ rij =E, as dual variables. With this choice of the unknown functions fui ; sij g, the first-order differential equations and the boundary conditions in the conventional formulation are described as follows. 1 2 ðui ;j

þuj ;i Þ  ð1 þ mÞsij þ msaa dij ¼ 0 in X

ð19aÞ

bi in X ð19bÞ E ð20aÞ ui ¼ u
i on Cui
ti sij nj ¼ ð20bÞ on Cti E With this system of equations, following the procedure in the previous section is straightforward. For the two-dimensional plane strain problem the following form of the conventional formulation can be obtained after some manipulations. sij;j ¼ 

200

u; x  ð1  m2 Þp þ mð1 þ mÞq ¼ 0 in X

ð21aÞ

v; y þ mð1 þ mÞp  ð1  m2 Þq ¼ 0 in X

ð21bÞ

u;y þv;x 2ð1 þ mÞr ¼ 0 in X

ð21cÞ

p; x þ r; y ¼ bx

in X

ð21dÞ

r; x þ q; y ¼ by

in X

ð21eÞ

on Ctn

tx nx p þ ty ny q þ ðtx ny þ ty nx Þr ¼
tt =E on Ctt

It is noted that ^sij is unsymmetric and the dimensionless stress sij is given as follows.

sij ¼ ^sij þ ^sji

ð28Þ

With these dual variables, the equation of equilibrium and the constitutive equation, i.e. the relationship between primal and dual variables, can be described in the following form.

bi ¼ 0 in X E 1 sij  m^saa dij 2 ui;j ¼ ð1 þ mÞ^

^sij;j þ ^sji;j þ

ð29Þ in X

ð30Þ

ui ;jk ¼ ui ;kj

ð22bÞ

þ 2nx ny r ¼
tn =E

4.2 Compatibility-imposed formulation In order to enforce the compatibility condition to firstorder least-squares formulation, new dual variables ^sij , defined as follows, are introduced. 1 m ^sij  ð27Þ ui ;j þ ua ;a dij 2ð1 þ mÞ 2ð1 þ mÞð1  2mÞ

ð22dÞ

v ¼ v
on Cv þ

ð26dÞ

ð22cÞ

ð22aÞ

n2y q

ð26cÞ

The well-known St. Venant’s form and Beltrami–Michell form of compatibility equations, which involve secondorder differentiation of stress components and strain components, respectively, are not suitable for first-order least-squares formulation. In the present work, the following form of compatibility condition is considered.

u ¼ u
on Cu n2x p

h i B3 ¼ 0 0 n2x n2y 2nx ny

B4 ¼ 0 0 tx nx ty ny tx ny þ ty nx

ð31Þ

where p, q and r are the dual variables, sxx , syy and sxy , respectively; t i is the unit tangential vector component on the boundary C; Ctn and Ctt are the normal traction boundary and the tangential traction boundary, respectively. This can be rewritten in the matrix form as follows.

At this point it can be mentioned that the reason for introducing the new dual variables ^sij is to express (31) only with first-order differentiations of dual variables. By differentiation of (30), the left-hand side and the right-hand side of (31) are written as follows.

Au ¼ f

ð23Þ

ui ;jk ¼ 2ð1 þ mÞ^sij;k  2m^saa;k dij

ð32aÞ ð32bÞ

in X

B1 u ¼ g1

on Cu

ð24aÞ

ui;kj ¼ 2ð1 þ mÞ^sik;j  2m^saa;j dik

B2 u ¼ g2

on Cv

ð24bÞ

B3 u ¼ g3

on Ctn

ð24cÞ

From (31) and (32), the compatibility condition for dual variables is obtained as the following form.

B4 u ¼ g4

on Ctt

ð24dÞ

ð1 þ mÞ^sij;k  m^saa;k dij ¼ ð1 þ mÞ^sik;j  m^saa;j dik

ð33Þ

where

2

o=ox 0 1 þ m2 6 0 o=oy m þ m2 6 A¼6 0 6 o=oy o=ox 4 0 0 o=ox 0 0 0

in X

m þ m2 1 þ m2 0 0 o=oy

3

0 7 0 7 2ð1 þ mÞ 7 7 o=oy 5 o=ox

Removing identities and repetitions from (33), only nine equations are linearly independent. Along with this compatibility condition, the following form of boundary conditions for dual variables is added.

n  rui ¼ n  r
ui

on Cui

ð34Þ

uT ¼ ½ u v p q r



f T ¼ 0 0 0 bx by

ð25bÞ

B1 ¼ ½ 1 0 0 0 0

ð26aÞ

These additional boundary conditions are the tangential derivatives of displacement on Cui . It should be noted that rui in (34) can be expressed by dual variables ^sij and thus (34) is essential boundary condition for dual variables. The compatibility-imposed formulation for the unknown functions fui ; ^sij g, is summarized as follows.

B2 ¼ ½ 0 1 0 0 0

ð26bÞ

ui ;j 2ð1 þ mÞ^sij þ 2m^saa dij ¼ 0 in X

ð25aÞ ð25cÞ

ð35aÞ

bi in X E ð1 þ mÞ^sij;k  ð1 þ mÞ^sik;j  m^saa;k dij þ m^saa;j dik ¼ 0 in X

ui ¼ u
i

Relative error in energy norm 2 31=2 ,2 31=2 Z Z 4 eT  r dX5 ¼ 4 ðe  eh ÞT  ðr  rh ÞdX5

ð35bÞ

^sij;j þ ^sji;j ¼ 

ð35cÞ

X

X

ð36aÞ

on Cui

ð40bÞ

ui on Cui ð36bÞ n  rui ¼ n  r

ti on Cti ð36cÞ ð^sij þ ^sji Þnj ¼ E For the two-dimensional plane strain problem the following form of the compatibility-imposed formulation is obtained after some manipulations. u; x  ð1  m2 Þp þ mð1 þ mÞq ¼ 0 in X

ð37aÞ

h

In the above, u is the exact solution and u the numerical approximation. In the presented figures, GALMF, CLSMF, ELSMF are used to denote the meshfree methods with Galerkin formulation, conventional least-squares formulation and compatibility-imposed least-squares formulation, respectively; AcItg and InacItg denote accurate integration and inaccurate integration, respectively.

r;x þs;x þq;y ¼ by =E in X

ð37fÞ

2r;x ð1  mÞp;y þ mq;y ¼ 0 in X

ð37gÞ

mp;x ð1  mÞq;x þ2s;y ¼ 0 in X

ð37hÞ

u ¼ u
on Cu

ð38aÞ

5.1 Cantilever beam Consider the cantilever beam problem as shown in Fig. 1, for which the following exact solution is given (Timoshenko and Goodier, 1987).   Py 1 2 2 ð41aÞ u ¼
ð6L  3xÞx þ ð2 þ m
Þ y  D 6EI 4  P 1 2 2 2 v ¼
3
my ðL  xÞ þ ð4 þ 5
mÞD x þ ð3L  xÞx 6EI 4 ð41bÞ

v ¼ v
on Cv

ð38bÞ

The corresponding stresses are given as follows.

n  ru ¼ n  r
u on Cu

ð38cÞ

n  rv ¼ n  r
v on Cv

ð38dÞ

2

v;y þmð1 þ mÞp  ð1  m Þq ¼ 0 in X

ð37bÞ

u;y 2ð1 þ mÞr ¼ 0 in X

ð37cÞ

v; x  2ð1 þ mÞs ¼ 0 in X

ð37dÞ

p; x þ r;y þs;y ¼ bx =E

ð37eÞ

n2x p

þ

n2y q

in X

þ 2nx ny ðr þ sÞ ¼
tn =E

on Ctn

tx nx p þ ty ny q þ ðtx ny þ ty nx Þðr þ sÞ ¼
tt =E

ð38eÞ on Ctt

P rxx ¼  ðL  xÞy I ryy ¼ 0   P D2 2 y rxy ¼ 2I 4

ð42aÞ ð42bÞ ð42cÞ

ð38fÞ

where I is the bending moment of inertia of the beam, 3 where the primal variables, u and v, and the dual variables, I ¼ D =12, and for plane strain case the elastic material constants are given in the following form. p, q, r and s, are defined as follows.

½u

v

p

q

r

s ¼ ux

uy

2^sxx

2^syy

^sxy

^syx ð39Þ


E ¼ E=ð1  m2 Þ

ð43aÞ

m
¼ m=ð1  mÞ

ð43bÞ

The problem is analyzed for the plane strain case with P ¼ 104 , E ¼ 5  107 , L ¼ 10 and D ¼ 2. Due to the 5 antisymmetry about y ¼ 0, only upper half of the beam is Numerical examples modeled. For the following convergence study, uniform Through numerical examples, the convergence characteristics of the least-squares meshfree method with the pre- nodal distributions are used with 3  21, 5  41, 9  81, sent formulations is studied and compared with Galerkin 17  161 and 33  321 equally spaced nodes. meshfree method. Also the robustness of the least-squares meshfree method to integration errors and to incompressible condition is illustrated. For the error estimation and convergence studies, the following relative errors are considered.

Relative error in displacement norm 2 31=2 ,2 31=2 Z Z 4 u  u dX5 ¼ 4 ðu  uh Þ  ðu  uh ÞdX5 X

X

ð40aÞ

Fig. 1. Cantilever beam problem

201

202

5.1.1 Convergence behavior The cantilever beam problem is solved with m ¼ 0:3 and three different types of boundary conditions: only displacement prescribed case; only traction prescribed case; both displacement and traction prescribed case as shown in Fig. 2. These boundary conditions are imposed by the exact solution of (41) and (42) on the corresponding boundaries. For accurate integration, the rectangular background cells, where the vertices of the cells coincide with the uniformly distributed nodes, are used with 4  4 Gaussian quadrature rule. For each boundary condition case, the relative errors are shown in Figs. 3–5. In Fig. 3, in which only displacement is prescribed, Galerkin meshfree method and both least-squares meshfree methods give the same rate of convergence for primal variables, e.g. displacements. However for dual variables, e.g. stress or strain, leastsquares methods give higher rate of convergence than Galerkin method. A little higher rate of convergence is observed in conventional formulation. The compatibilityimposed formulation shows the same rate of convergence, which is optimal rate, for both primal and dual variables. This means that the solution satisfies the regularity requirements for H1 -ellipticity of the compatibility-imposed formulation. These results are consistent with the theorems about the convergence characteristics of leastsquares formulations in Sect. 3.3. Although there are some changes in the order of solution accuracy among the formulations, from the viewpoint of convergence similar

results are observed in other boundary condition cases as shown in Figs. 4 and 5. The convergence test was also performed with other conventional formulations and compatibility-imposed formulations with different dual variables, which gave similar behaviors with the presented results.

5.1.2 Effect of inaccurate integration For studying the effects of inaccurate integration in Galerkin meshfree method and least-squares meshfree method, the cantilever beam is inclined by 30 and is solved using the rectangular background cells of which vertices do not coincide with nodes as shown in Fig. 6. With the increase of nodes, the cells are refined so that each cell has the size of nodal distance, same as the case of accurate integration, and 2  2 Gaussian quadrature rule is used to magnify the effects of inaccurate integration. In the background cells, the integration points

Fig. 2. Boundary conditions of beam problem: a only displace- Fig. 3. Errors in displacement and energy norms for the beam ment b.c. is prescribed; b only traction b.c. is prescribed; c both problem (m ¼ 0:3) with accurate integration when only displacedisplacement b.c. and traction b.c. are prescribed ment b.c. is prescribed

203

Fig. 4. Relative errors for the beam problem (m ¼ 0:3) with accurate integration when only traction b.c. is prescribed

Fig. 5. Relative errors for the beam problem (m ¼ 0:3) with accurate integration when both displacement b.c and traction b.c are prescribed

outside domain are removed. It is noted that this integration scheme is simple but does not support accurate integration. The both displacement and traction boundary conditions are enforced (see Fig. 2c) and m ¼ 0:3 is used in this study. Figure 7 compares the relative errors using Galerkin formulation and least-squares formulations in meshfree methods. It is observed that inaccurate integration does not deteriorate the solution accuracy and the convergence behavior in the least-squares methods while the solution of Galerkin method is degraded seriously and even does not seem to converge, by the comparison with the results of accurate integration Fig. 6. Background integration cells for the beam problem with in Fig. 5. inclined domain

5.1.3 Effect of incompressibility In order to show the robustness of least-squares meshfree solution in the incompressible limit, the

cantilever beam problem with both displacement and traction boundary conditions is analyzed with different values of Poisson ratio; m = 0.3, 0.49, 0.4999 and

204

Fig. 7. Relative errors for the beam problem (m ¼ 0:3) with inaccurate integration when both displacement b.c. and traction b.c. are prescribed

0.499999. As Poisson ratio approaches to 0.5, Galerkin method gives poor solution accuracy as shown in Fig. 8, similar results can be also found in other literatures (Dolbow and Belytschko, 1999; Chen et al., 2000a, b). It is expected that the least-squares formulations presented in the previous section do not rise any singularity in the vicinity of m = 0.5. This expectation is approved in Figs. 9 and 10, which show the relative errors with m approaching to 0.5 for conventional formulation and compatibility-imposed formulation, respectively. Little change is observed in solution accuracy and convergence rates for both primal and dual variables are maintained even in the case of m = 0.5. Figures 11 and 12 show the results when using the inaccurate integration scheme described in Fig. 6, and the solution behavior is also good in this situation. For the sake of comparison, the relative errors of Galerkin and both least-squares methods are displayed for m ¼ 0:499999 in Fig. 13.

Fig. 8. Relative errors for the beam problem using Galerkin meshfree method (accurate integration) with different values of m

5.2 Infinite plate with a circular hole The exact solution for an infinite plate with a hole of radius a centered at origin, subjected to a unit tension in the x-direction at infinity (see Fig. 14), is given in Timoshenko and Goodier (1987) as follows.   1 j1 a2 r þ cosð2hÞ þ f1 þ ð1 þ jÞ cosð2hÞg ur ¼ r 4l 2  4 a ð44aÞ  3 cosð2hÞ r  1 a2 a4 ð1  jÞ  r  3 sinð2hÞ ð44bÞ uh ¼ r r 4l The corresponding stresses are written as follows.

rxx

 a2 3 3a4 cosð2hÞ þ cosð4hÞ þ 4 cosð4hÞ ¼ 2 r 2 2r ð45aÞ

205

Fig. 9. Relative errors for the beam problem using conventional least-squares meshfree method (accurate integration) with different values of m

ryy ¼ 

 a2 1 3a4 cosð2hÞ  cosð4hÞ  4 cosð4hÞ 2 r 2 2r

 a 1 3a4 sinð2hÞ þ sinð4hÞ þ 4 sinð4hÞ ¼ 2 r 2 2r

ð45bÞ

2

rxy

ð45cÞ where the material constants for plane strain case are as follows.

l¼

E 2ð1 þ mÞ

j ¼ 3  4m

ð46aÞ ð46bÞ

Due to symmetry, only the upper quadrant of the problem, 0  x  5 and 0  y  5, is modeled with a ¼ 1 as shown in Fig. 15. This problem is analyzed

Fig. 10. Relative errors for the beam problem using compatibility-imposed least-squares meshfree method (accurate integration) with different values of m

with E ¼ 1  107 and only traction-prescribed boundary condition from the exact solution is considered with enforcing three displacements for removing rigid-body mode as shown in Fig. 15. The nodal distribution for convergence study is shown in Fig. 16: 49, 169, 625 and 2401 nodes.

5.2.1 Convergence behavior For the study of convergence behavior under accurate integration, the background cells are such that the nodes are located at the vertices of the integration cells with 4  4 Gaussian quadrature rule. Figure 17 shows the results in the case of m ¼ 0:3. The trends in convergence are similar with the results of the beam problem. The same rates of convergence are observed for primal variables in both Galerkin method and least-squares methods, and the optimal rate of convergence is obtained for dual variables

206

Fig. 12. Relative errors for the beam problem using compatibilFig. 11. Relative errors for the beam problem using conventional ity-imposed least-squares meshfree method (inaccurate integraleast-squares meshfree method (inaccurate integration) with tion) with different values of m different values of m

using compatibility-imposed formulation, though in this case conventional formulation does not exhibit higher convergence rate than Galerkin formulation.

5.2.2 Effect of inaccurate integration Figure 17 also includes the results under inaccurate integration, for which the quadtree algorithm is used to generate the integration points as follows (see Fig. 18): (i)

Let Nmax be a properly chosen integer, which is not too large. In the present work Nmax ¼ 4 was used. (ii) Take a square which covers the domain of analysis. (iii) If the square has more than Nmax nodes in it, divide the square into four squares with equal size until each square has nodes less than or equal to Nmax . (iv) Considering the number of nodes in each square, generate integration points in the square applying an adequate Gaussian quadrature rule (in the present work, the rule in Table 1 is used).

(v) Remove the integration points which are outside the domain. It is noted that these simply constructed quadtree cells are not exactly matched with the analysis domain and the integration is not accurate especially in the vicinity of the boundaries. Using this inaccurate integration scheme, the solution of Galerkin method does not seem to converge as the number of nodes increases. Although the least-squares methods show a little degradation especially in the solutions for primal variables as compared with those using accurate integration, they still work well and give robust results.

5.2.3 Effect of incompressibility The problem of the infinite plate with a hole is solved with different values of Poisson ratio; m ¼ 0.3, 0.49 and 0.4999. It is clear in Fig. 19 that in the vicinity of m ¼ 0:5 Galerkin method exhibits poor solution accuracy due to

207

Fig. 15. Quarter model for the plate problem

Fig. 13. Relative errors for the beam problem with m ¼ 0:499999

Fig. 16. Nodal distribution of the plate problem

Fig. 14. Infinite plate with a circular hole problem

as shown in Figs. 20 and 21, respectively. For the purpose of comparison, the solutions of Galerkin method and both the least-squares methods with the refinement of nodal distribution are shown for m ¼ 0:4999 in Fig. 22.

the singularity in the nearly incompressible condition. However the least-squares methods are hardly affected and still result in good solution behavior even in the case of m ¼ 0:5 with both the accurate integration and the inaccurate integration using the qaudtree algorithm

5.3 L-shaped beam The solutions of the previous two examples are sufficiently smooth to achieve the theoretical rate of convergence demonstrated in (13) and (14). To study how the leastsquares meshfree method performs for problems with

Table 1. Gaussian quadrature rule for quadtree algorithm No. of nodes in cell

Quadrature order

0 1 2 3 4

2 3 3 4 4

· · · · ·

2 3 3 4 4

208

Fig. 17. Relative errors for the plate problem with m ¼ 0:3

Fig. 19. Relative errors for the plate problem using Galerkin meshfree method (accurate integration) with different values of m

Fig. 18. Quadtree algorithm

non-smooth solution, the problem of L-shaped domain shown in Fig. 23 is considered. Displacement is fixed on the left side of the domain where x ¼ 0, and unit

displacement in y-direction is prescribed on the right side of the domain. The plane strain is assumed with E ¼ 1  108 and m ¼ 0:3. This problem is analyzed with different values of r, which is the radius of the corner curvature. As the radius of curvature decreases to zero, singularity occurs at both ends of A–A section. In this case, the regularity requirement of (13) and (14) is not satisfied. In this study, the results of the least-squares methods are compared with those of Galerkin method. Integration

209

Fig. 20. Relative errors for the plate problem using least-squares meshfree methods (accurate integration) with different values of m Fig. 21. Relative errors for the plate problem using least-squares meshfree methods (inaccurate integration) with different values of m

is performed using the quadrilateral background cells whose vertices coincide with the distributed nodes. In each cell 4  4 Gaussian quadrature rule is used in Galerkin method and 2  2 rule in the least-squares methods. Figure 24 shows the magnitude of displacement and vonMises stress along A–A section in the case of r ¼ 1. Since both least-squares methods gave the similar results, only the result of the compatibility-imposed formulation is displayed. As the number of nodes is increased, the results of the least-squares method and Galerkin method are in good agreement. The problem is also solved with decreasing the radius of curvature; r ¼ 0.5, 0.25 and 0. As r approaches to zero, the stress increases abruptly at both ends of A–A section. The results in Fig. 25 were obtained using 12681 nodes for r ¼ 0:5 and 13985 nodes for both r ¼ 0.25 and 0. In all cases, the magnitudes of displacement are in good agreement. However, with r approaching to zero, Galerkin method gives higher von-Mises stress at the ends of A–A

section than the least-squares method. This is attributed to the fact that the stress field is computed by enforcing the stress-displacement relationship in least-squares sense. The stress field by the least-squares method is smoother but less sensitive to steep change than that by Galerkin method, which shows some oscillation as shown in Fig. 25. This implies that the least-squares method requires more number of nodes at the singular region than Galerkin method to capture the steep change of stress field.

6 Concluding remarks A truly meshfree method based on first-order leastsquares formulations for linear elasticity and the moving least-squares approximation has been proposed in the present work. In the present method, a simple integration algorithm can be effectively used and no additional structure of cells for alleviating the incompressible

210

Fig. 22. Relative errors for the plate problem with m ¼ 0:4999

Fig. 24. Magnitude of displacement and von-Mises stress along A–A section with r = 1

employs dimensionless stress components as dual variables to transform the second-order differential governing equations into the first-order differential equations. The compatibility-imposed formulation introduces similar dimensionless dual variables to enforce the modified compatibility equations, which are the first-order differential equations of the dual variables. Since the leastsquares formulations do not involve integral identities unlike Galerkin formulation, they are highly robust to integration errors so that simply constructed background integration cells can be effectively used. In the first-order least-squares formulations, both primal and dual variables can be approximated by the shape functions of same order and moreover the system matrix is always positive-definite. Therefore the present method does not Fig. 23. L-shaped beam problem require any scheme to remove the incompressible locking, and exhibits the uniform convergence as the material locking is used. The conventional and the compatibility- tends to become incompressible. The numerical results have shown that the least-squares imposed least-squares formulations for linear elasticity meshfree method works well while Galerkin meshfree have been presented. The conventional formulation

Beissel S, Belytschko T (1996) Nodal integration of the elementfree Galerkin method. Comp. Meth. Appl. Mech. Eng. 139: 49–74 Belytschko T, Lu YY, Gu L (1994) Element-free Galerkin methods. Int. J. Num. Meth. Eng. 37: 229–256 Breitkopf P, Touzot G, Villon P (2000) Double grid diffuse collocation method. Comput. Mech. 25: 199–206 Cai Z, Manteuffel TA, McCormick SF (1997) First-order system least squares for the Stokes equations, with application to linear elasticity. SIAM. J. Numer. Anal. 34: 1727–1741 Cai Z, Manteuffel TA, McCormick SF, Parter SV (1998) Firstorder system least squares (FOSLS) for planar linear elasticity: pure traction problem. SIAM. J. Numer. Anal. 35: 320–335 Cai Z, Lee CO, Manteuffel TA, McCormick SF (2000a) First-order system least squares for the Stokes and linear elasticity: numerical results. SIAM. J. Sci. Comput. 5: 1706–1727 Cai Z, Lee CO, Manteuffel TA, McCormick SF (2000b) First-order system least squares for the Stokes and linear elasticity equations: further results. SIAM. J. Sci. Comput. 5: 1728–1739 Chen JS, Yoon S, Wang HP, Liu WK (2000a) An improved reproducing kernel particle method for nearly incompressible finite elasticity. Comp. Meth. Appl. Mech. Eng. 181: 117–145 Chen JS, Wang HP, Yoon S, You Y (2000b) Some recent improvement in meshfree methods for incompressible finite elasticity boundary value problems with contact. Comput. Mech. 25: 137–156 Chen JS, Wu CT, Yoon S, You Y (2001) A stabilized conforming nodal integration for Galerkin mesh-free methods. Int. J. Num. Meth. Eng. 50: 435–466 De S, Bathe K-J (2001) Displacement/pressure mixed interpolation in the method of finite spheres. Int. J. Num. Meth. Eng. 51: 275–292 Dolbow J, Belytschko T (1999) Volumetric locking in the element free Galerkin method. Int. J. Num. Meth. Eng. 46: 925–942 Duarte CA, Oden JT (1995) Hp clouds – a meshless method to solve boundary value problems, Technical Report 95-05, TICAM, University of Texas at Austin Jiang BN (1998) The Least-Squares Finite Element Method – Theory and Application. Springer-Verlag, Berlin Lancaster P, Salkauskas K (1981) Surfaces generated by moving least-squares methods. Math. Comput. 37: 141–158 Liu WK, Jun S, Zhang YF (1995) Reproducing kernel particle methods. Int. J. Num. Meth. Fluids. 20: 1081–1106 Fig. 25. Magnitude of displacement and von-Mises stress along Melenk JM, Babusˆka I (1996) The partition of unity finite element A–A section with different values of r method: basic theory and applications. Comp. Meth. Appl. Mech. Eng. 139: 289–314 Nagashima T (1999) Node-by-node meshless approach and its method does not converge with simple and inaccurate applications to structural analysis. Int. J. Num. Meth. Eng. 39: schemes. In most cases, the least-squares meshfree method 341–385 gives higher rate of convergence for dual variables than Nayroles B, Touzot G, Villon P (1992) Generalizing the finite Galerkin meshfree method, and the optimal rate of conelement method: diffuse approximation and diffuse elements. vergence is achieved in the compatibility-imposed forComput. Mech. 10: 307–318 mulation. However, for this optimal performance some On˜ate E, Zienkiewicz OC, Taylor RL (1996) A finite point method in computational mechanics. Applications to convective regularity of the solution should be met. In the present transport and fluid flow. Int. J. Num. Meth. Eng. 39: 3839– formulations no volumetric locking is observed and 3866 the convergence characteristics is valid also in the Park SH, Youn SK (2001) The least-squares meshfree method. incompressible limit. Int. J. Num. Meth. Eng. 52: 997–1012 Park SH, Kwon KC, Youn SK (2002) A study on the convergence References of least-squares meshfree method under inaccurate integraAtluri SN, Zhu T (1998) A new meshless Petrov-Galerkin (MLPG) tion. Int. J. Num. Meth. Eng. in printing approach in computational mechanics. Comput. Mech. 22: Timoshenko SP, Goodier JN (1987) Theory of Elasticity, 3rd edn. 117–127 McGraw-Hill, New York

211