World Journal of Mechanics, 2012, 2, 339-360
doi:10.4236/wjm.2012.26040 Published Online December 2012 (http://www.SciRP.org/journal/wjm)
Copyright © 2012 SciRes. WJM
Elastoplastic Large Deformation Using Meshless Integral
Method
Jianfeng Ma1, X. J. Xin2
1Department of Aerospace and Mechanical Engineering, Saint Louis University, Saint Louis, USA
2Department of Mechanical and Nuclear Engineering, Kansas State University, Manhattan, USA
Email: jma15@slu.edu, xin@ksu.edu
Received September 23, 2012; revised October 24, 2012; accepted November 4, 2012
ABSTRACT
In this paper, the meshless integral method based on the regularized boundary integral equation [1] has been extended to
analyze the large deformation of elastoplastic materials. The updated Lagrangian governing integral equation is ob-
tained from the weak form of elastoplasticity based on Green-Naghdi’s theory over a local sub-domain, and the moving
least-squares approximation is used for meshless function approximation. Green-Naghdi’s theory starts with the addi-
tive decomposition of the Green-Lagrange strain into elastic and plastic parts and considers a J2 elastoplastic constitu-
tive law that relates the Green-Lagrange strain to the second Piola-Kirchhoff stress. A simple, generalized collocation
method is proposed to enforce essential boundary conditions straightforwardly and accurately, while natural boundary
conditions are incorporated in the system governing equations and require no special handling. The solution algorithm
for large deformation analysis is discussed in detail. Numerical examples show that meshless integral method with large
deformation is accurate and robust.
Keywords: Meshless Method; Large Deformation; Local Boundary Integral Equation; Moving Least-Squares
Approximation; Subtraction Method; Singularity Removal; Elastoplasticity; Green-Naghdi’s Theory
1. Introduction
Over the past two decades the meshless methods have
attracted much attention owing to their advantages in
adaptivity, higher degree of continuity in the solution
field, and the capability to handle moving boundaries and
changing geometry. In the meshless method, the concept
of an element is eliminated. The model geometry consists
of a distribution of nodes over the problem domain, and
the approximate solution is constructed entirely based on
these nodes. Most development in meshless methods to
date has been focused mostly on linear elasticity. Re-
search in large deformation plasticity using meshless
methods is only currently gaining attention.
For the formulation of elastoplastic theory at large de-
formation, Green-Naghdi’s theory [2,3] begins with the
additive decomposition of the Green-Lagrange strain into
elastic and plastic parts as ep
K
LKLKL
EEE
. On the
other hand, Lee’s theory [3,4] considers the multiplica-
tive decomposition of deformation gradient F into an
elastic Fe and plastic part Fp as ep
F
FF. Both theo-
ries are formulated on the basis of fundamental laws of
continuum mechanics. According to Chiou et al. [3],
Green-Naghdi’s theory is more flexible because it can be
applied to either isotropic or anisotropic materials and
the computing procedure involved is relatively straight-
forward. In general, theories of large deformation plas-
ticity have to be implemented numerically since it is dif-
ficult to obtain closed form solutions to practical prob-
lems involving large deformation. Researchers such as
Chiou et al. [5], Lee [6], and Hu [7], have used FEM to
solve large deformation plasticity problems.
Some researchers have used meshless methods to ana-
lyze large deformation with elastoplasticity. Belytschko
et al. [8] proposed a 3D Element-Free Galerkin (EFG)
method intended for dynamic problems with geometric
and material nonlinearities solved with explicit time in-
tegration. Rossi and Alves [9] applied a modified Ele-
ment-Free Galerkin method to large deformation proc-
esses. The proposed EFG method enables direct enforce-
ment of essential boundary conditions, and the plasticity
model assumes a multiplicative decomposition of the
deformation gradient into an elastic and plastic part and
allows for nonlinear isotropic hardening. Chen et al. [10]
and Eskandarian et al. [11] formulated the governing
equations for rate-independent large strain plasticity in
the framework of meshless method and used the method
to simulate high-speed impact. Xiong et al. [12] used
meshless local Petrov-Galerkin method to model large
deformation of micro-electronic mechanical devices
(MEMS) by using local radial point interpolation (RPIM)
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
340
approximation to construct shape function. The system
equation is obtained based on the Total Lagrangian (TL)
approach. Hu et al. [13] used meshless local Petrov-
Galerkin method for large deformation contact analysis
of elastomeric components. A nonlinear formulation of
meshless local Petrov-Galerkin finite-volume mixed me-
thod was developed by Han et al. to analyze static and
dynamic large deformation problems [14]. Li et al. [15]
developed a coupled finite element and meshless local
Petrov-Galerkin method to analyze large deformation
problems. Gu et al. [16,17] employed meshless local
Kriging method to analyze large deformation problems.
A comparison between Total Lagrangian and Updated
Lagrangian approaches with an adaptive local meshless
formulation was given by Gu in [18]. Gun et al. [19]
used a meshless formulation of the Euler-Nernouli theory
with non-linear strain-displacement relation which covers
both axial and bending deformations to predict the large
deformation behavior of fabrics. Li and Lee [20,21] em-
ployed an adaptive meshless method to solve mechanical
contact problems which incorporated a sliding line algo-
rithm with penalty method to handle contact constrains.
In [22], Zhu et al. used meshless local natural neighbor
interpolation method to analyze 2D elastoplastic large
deformation in conjunction with Updated Lagrangian
formulation. A Galerkin SPH method was utilized by
Wang to solve large deformation fracture problems [23].
Reproducing kernel particle methods were used to ana-
lyze large deformation problems in [24-26]. In [27], Li et
al. used an element-free Galerkin method to solve die
forging problems. Quak et al. [28] applied both SPH and
EFG meshless methods to analyze extraction problems.
The authors have developed a meshless integral method
for linear elasticity [1] and later extended it to elastoplas-
ticity for small deformation [29]. This meshless integral
method is truly meshless and does not require a back-
ground mesh for integration. In the present work, we
extend the meshless integral method to large deformation
elastoplasticity. The governing integral equation is ob-
tained from the weak form based on Green-Naghdi’s
large deformation elastoplasticity theory over a local sub-
domain, and moving least-squares approximation is used
for meshless function approximation. The constitutive
law uses A J2 elastoplastic flow theory based on von
Mises yield criterion with isotropic hardening in which
the Green-Lagrange strain is related to the second Piola-
Kirchhoff stress. The fixed point iteration, because of its
simplicity in implementation, is used to solve the nonlin-
ear equations arising from large deformation elastoplas-
ticity. A number of numerical tests were carried out for
model verification. It was found that the meshless results
are in excellent or satisfactory agreement with either
closed form solutions or FEM results, indicating that the
large deformation elastoplasticmeshless integral method
based on the regularized integral equation is accurate and
robust. An application of the current method to the metal
forming process has been reported separately [30].
This paper is organized as follows: In Section 2, the
regularized local boundary integral equation is derived,
and the subtraction method is used to remove the strong
singularity that is present in the initial local boundary
integral equation. In Section 3, large deformation elasto-
plastic constitutive equations are presented. Section 4
describes the moving least-squares approximation (MLSA)
for approximating the displacement, strain, and stress
fields over the problem domain. The meshless imple-
mentation of the regularized boundary integral equation
using MLSA and the treatment of weak singularity are
presented in Section 5. Section 6 discusses the enforce-
ment of essential and natural boundary conditions. The
solution algorithm for solving the nonlinear system equa-
tions is discussed in Section 7. Numerical examples are
presented in Section 8 to assess the accuracy and effec-
tiveness of the method. Discussion and conclusions from
this study are given in Section 9.
2. Regularized Local Boundary Integral
Equation Using Subtraction Method
In this work, the updated Lagrangian coordinate system
is used in which the stress directions are referred to the
last known equilibrium state. Consider a two-dimen-
sional body for a given load increment as shown in Fig-
ure 1. The equilibrium state of the body at the beginning
of the load step, which has been solved and is denoted by
0, is referred to as the reference configuration. The ref-
erence configuration designates the configuration with
respect to which subsequent deformation is measured.
The domain of the current configuration of the body at
the end of the load increment is denoted by , which is
also called the deformed configuration. The vector X for
a given material point in the reference configuration does
not change with time, and X is called the Lagrangian
coordinates; x, which describes the material point in the
current configuration, changes with time, and x is called
the Eulerian coordinates. In this work, the rectangular
Figure 1. Undeformed (initial or reference) and deformed
(current) configuration.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
341
Eulerian coordinates, xk (k = 1, 2), and Lagrangian coor-
dinates, XK (K = 1, 2), are employed. We use upper case
letter and upper case subscript to indicate that the vari-
able of interest refers to the reference state, and use lower
case letter and lower case subscript to indicate that the
variable of interest refers to the current state. Within the
updated Lagrangian framework, the following proce-
dures are implied: 1) During each load increment, all
field variables are defined with respect to the state at the
start of this load increment (reference configuration); 2)
At the end of this load increment, field variables are up-
dated with respect to the state at the end of this load in-
crement (current configuration), and the current configu-
ration will be the reference configuration for the next
load increment.
In large deformation problems, various stress measures
can be defined. The most commonly used stresses are: 1)
Cauchy stress (true stress) tensor σ which expresses the
stress relative to the current configuration; 2) First Piola-
Kirchhoff stress (nominal stress) tensor P which relates
forces in the current configuration with areas in the
reference configuration; 3) Second Piola-Kirchhoff stress
tensor T which relates forces in the reference configuration
to areas in the reference configuration, where the force in
the reference configuration is obtained via a mapping
that preserves the relative relationship between the force
direction and the area normal in the current configuration.
The relationships between these stresses are
11T
,,JJ

PFσTFσF
where F is the deformation gradient and defined as
x
x
X
Y
yy
X
Y










F, and J is the determinant of F. The
second Piola-Kirchhoff stress tensor is symmetric, and in
the current work is linked to the Green-Lagrange strain in
the elastoplastic constitutive law.
For a large deformation elastoplastic body represented
by a reference domain 0 with boundary Γ0, the govern-
ing elastoplasticity equations are as follows:
 
  

 
,0
,,,,
0
0
1,
2
IJ JI
IJI JJ IKIK J
ecor
IJIJKL KLIJ
f
E uuuu
dTCdEdT





XX
X
XX XX
X
XXX
(1)
where τ (bold face denotes vectors or tensors) is the
stress measure defined by r
IJ rIrJ rILJ
L
x
GGT





,
TLJ is the second Piola-Kirchhoff stress, G is the trans-
formation tensor between the reference configuration and
the current configuration, r
L
x
X
is the deformation gra-
dient, E is the Green strain tensor associated with the
displacement u, ρ0 is the mass density in the reference
configuration 0, fI is the body force per unit mass,
e
I
JKL
C is the linear elastic stiffness matrix, and cor
I
J
dT is
the stress correction term because of nonlinearities in
either geometry or material. An index I following a
comma designates partial differentiation with respect to
XI, and repeated indices indicate summation over the
dimensionality of the problem. The essential and the
natural boundary conditions on the boundary Γ are re-
spectively:
on
I
Iu
uu
(2)
on
I
IJ JIt
TnT

(3)
here, u represents the prescribed displacement on u
,
T represents the prescribed traction on t
; n is the
outward unit normal to the boundary; tu
 and
tu
.
The weak form of (1) over a local sub-domain
00
a
 in the reference state is:
 


0
,0 0
,d 0
a
a
IJ JII
fg


XXXYX
(4)
where the notation
in this paper is used to denote a
node (e.g., a indicates node a, b indicates node b,
etc.) in order to reserve the usual subscripts, I, J, etc., for
denoting degree of freedom (DOF) components,
00
a
 is a spherical (circular in 2D) sub-domain
related to node a, a
Y is the position vector of node
a which is also called a source point, X is the integration
or field point which may or may not coincide with a node,
and
I
g
is the test function. In the following, the func-
tional dependence on X, i.e. “(X)”, will be dropped for
brevity when no ambiguity is caused. In this work, fol-
lowing [31], we use a special test function defined by

,,
aaa
IJIJ
gue
XYXYY (5)
where
J
e represents the J-th component of a unit force
vector,
J
I
u
is the special test function, given by [31]:

 
2
2
2
22
,
1
8π1
54
43ln 1
23 4
1
a
IJ
aa
I
J
aa
ss
aa
a
IJ
aa
s
u
rvr
v
hh
rr
r
hr




 








XY
(6)
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
342
where a
r is the distance from a
Y to X;
aa
I
II
rXY; a
n is the outward unit normal to the
boundary 0
a
 at X; a
s
h is the radius of the local
sub-domain;
and EE
for plane strain, or
1


and

2
1EE
 for plane stress; E is
the Young’s modulus,
is the Poisson’s ratio, and


21E

is the shear modulus. Note that the
special test function is defined in the reference state. The
associated traction, for linear elastic behavior, is







2
32
2
,
11
12
4π134
12
2
3
34
a
IJ
aa
a
KK
I
J
aa a
s
aa
a
aa aa
IJ
IJ JI
aa a
aa aa
IJ JI
a
s
t
rn
rv
vrn vh
rr
rv
rnrn
nr r
rnrn
vh








XY
(7)
The special test function
I
J
u
has the property that it
vanishes on the boundary of a spherical 0
a
. With
I
J
u
as the test function, application of integration by parts to
(4) twice leads to:

 

 

 

 



 

 


0
0
0
0
0
0
0
*0
,
,,,
0
,
,
,d
,d
,d
,d
,d
,d
,d
a
a
a
a
a
a
a
a
I
a
IJ J
a
IJ J
acor r
IJ KLKrJ
L
ae
IJKRJR LLKMNMN
a
IJLKJ LK
a
IJLKJLJ LK
u
uT
tu
x
uT
X
uuCu
uTun
uT un












XYX X
XYX X
XYX X
XY XX
XY X
XYXX X
X
YX XX

 
0
0
,d0
a
a
IJ J
uf

XYX X
(8)
where 0
a
 is the boundary of 0
a
, 0
J
T is the J-th
component of the traction at the reference state,
I
u
is
the displacement increment from the reference state to
the current state, 0
L
K
T is the second Piola-Kirchhoff
stress at the reference state,
L
K
T
is the second Piola-
Kirchhoff stress increment from the reference state to the
current state, and
,a
XY is the Dirac delta function.
Equation (8) in its current form cannot be used directly in
numerical calculation because when a
Y is a boundary
node, the equation contains strong singularity (1/r type in
line integral) in the traction term
I
J
t
. For the local inte-
gral equation approach to be a valid numerical method,
the strong singularity must be handled appropriately.
The subtraction technique is employed in the present
study to remove the strong singularity; the technique for
small deformation has been presented in 4[1]. For simplic-
ity in implementation, the local sub-domain is always
chosen, in the reference state, as a sphere or part of a
sphere centered on a node. If node a
Y is an interior
node, a
s
h is selected such that 0
a
stays fully inside
0
. If a
Y is a boundary node, then 0
a
is the in-
tersection of 0
and a sphere a
 of radius a
s
h,
centered at the boundary node, and the boundary 0
a
is the union of the part of a
 inside 0
, denoted
by a
C, and the part of 0
 inside a
 , denoted
by Γa, as illustrated in Figure 2. A further modifica-
tion is the exclusion of a tiny sphere
of radius Δ
(which later tends to zero) centered on a
Y whose
boundary is denoted by a
C. Figure 3 shows sche-
matically this modification. For interior nodes, Γa is
empty, and a
C is a full circle.
With the above notations, the boundary 0
a
 can be
decomposed into the following sections:
0Γ
aaaa
CC
  (9)
aaa
ut
  (10)
where a
u
is the section of a
where the displace-
ment is prescribed, and a
t
the section where the trac-
tion is prescribed.
Using subtraction method, we obtain, in the limit of
0
, the following:
Figure 2. Schematic diagram showing the sub-domain for
an interior or a boundary node a
Y.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
343
Figure 3. Exclusion of a tiny sphere Δ of radius Δ centered
at a node for removing the strong singularity.

 



 

 



 







0
0
0
0
,
,,,
0
,
,
0
,d
,d
,d
,d
,d
,d
,d
a
a
a
a
a
a
acor r
IJ KLKrJ
L
ae
IJKRJR LLKMNMN
a
IJLKJ LK
a
IJLKJLJ LK
a
IJ J
aa
IJJ J
C
aa
IJJ J
x
uT
X
uuCu
uTun
uT un
uT
tuu
tuu









XY XX
XY X
XYXXX
X
YX XX
XYX X
XY XX
XY X


 
0
0
,d0
a
a
a
IJ J
uf

X
XYX X
(11)
The integration of

,a
IJ
t
XY over a
C can be
obtained in closed form:

,dΓ
a
aa
I
JIJ
C
t

XY (12)
with
 
 
212 1
21 21
sin 2sin 2cos2cos 2
2π2π34 2π34
cos 2cos 2sin2sin2
2π34 2π2π34
a
IJ
 


















(13)
Here 21

 is the internal boundary angle sub-
tended by material at a
Y on the boundary, as shown in
Figure 4.
Figure 4. Schematic diagram showing the internal bound-
ary angle θ2θ1 = θ at node a
Y on the boundary.
Several special cases are worth noting. For an interior
node, 2π
; for a boundary node where the boundary
is smooth, π
; for a corner node, θ = corner angle.
Substituting the above expressions into (11) leads to:

 



 

 

 

 


0
0
0
0
0
,
,,,
0
,
0
0
,
,d
,d
,d
,d
,d
,d
a
a
a
a
a
a
ut
a
aaa cor r
IJJIJ KLKRJ
L
ae
IJKRJR LLKMNMN
a
IJLKJ LK
a
IJ J
a
IJ J
a
IJLKJLJ LK
IJ
x
uuT X
uuCu
uTun
uf
uT
uT un
t




 





X
YX X
XY X
XYXXX
XYX X
XYX X
X
YX XX
X

 


,d
,()d
a
a
a
ut
a
J
C
aa
IJJ J
u
tuXu


YXX
XY X
(14)
In (14), the subtraction method has been used in the
last term on the right hand side. When the field point X
approaches the source node a
Y,

a
I
I
ux u tends to
zero which removes the strong singularity and makes the
integral numerically integrable. All other terms in (14)
are regular or weakly singular for which special integra-
tion quadrature gives convergent and accurate results.
The regularized Equation (14) holds for any source node
a
Y, either inside the domain or on the boundary. In the
current meshless integral method, both boundary and
interior source nodes are used, and the moving least-
squares approximation is employed for approximating
the solution field, as presented later.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
344
3. Constitutive Equation for Elastoplasticity
with Large Deformation
In this research, the Green-Naghdi’s theory with von
Mises yield criterion, associated flow rule, and isotropic
strain hardening is used to model the material behavior.
The Green-Naghdi’s theory is attractive because it can be
applied to either isotropic or anisotropic materials and
the computing procedure involved is relatively straight-
forward.
Green-Naghdi’s theory [2,3] begins with the additive
decomposition of the Green-Lagrange strain increment
into the elastic and plastic parts as ep
K
LKLKL
dEdE dE
.
The von Mises yield criterion has the following form:

2
10
2IJ IJp
fTT

 (15)
f is the yield function constructed in the space of second
Piola-Kirchhoff stress,
I
J
T is the component of devia-
toric stress of the second Piola-Kirchhoff stress, and
p

characterizes the hardening of the material. For
linear hardening,

03,
pp
YH

 where 0
Y is
the initial yield stress,
p
H
is the hardening modulus,
and
p
is the effective plastic strain.
p
H
is given by
the following equation
1
T
p
TY
E
HEE
(16)
where T
E is elastoplastic tangent modulus and Y
E is
the elastic Young’s modulus from uniaxial tension test.
The plastic strain increment
p
I
J
dE is expressed by an
associated flow rule as:
p
IJ
I
J
IJ
dL f
dE
H
T
dL T
H
(17)
where L is the loading criteria and given by the following
equation:




T
11 11222233 331212
2
1
f
dLC dE
T
ETdE TdETdETdE
v






(18)
[C] is elastic stiffness matrix which works for both
plane strain and plane stress.
The positive scalar H is expressed as:
2d
22
33d1
e
e
p
E
Hv





(19)
Using Equations (18) and (19), we obtain the plastic
strain increment:

11
22
33
12
2
p
IJ
T
T
dL
dE T
H
T







(20)
With
p
I
J
dE , the increment of effective plastic strain is
computed as follows:
2
3
22 2
33 3
pp
p
IJIJIJ IJ
eq eq
dL
dCdEdET T
H
dL dL
HH



(21)
where
3
2
eqIJ IJ
TT
is the equivalent stress of the
current stress state.
The constitutive equations are summarized as follows:
1) Elastic region
0
IJ
fT

223
1
2
KK
IJ IJIJ
I
KJL JKIL ILJKJLIK
IJ KLKL
dT
dT dE
dE
 




  
(22)
Where
and
are Lame constants.
2) Plastic region (loading,
0, 0
IJ
fT dL
)

2
2
9
23
1
2
IJ
K
LKLIJ
IJKK IJ
pe
I
KJL JKIL ILJKJLIK
IJ KLKL
dT
TdET
dE dEH
dE



 
 
(23)
3) Plastic region (unloading,
0, 0
IJ
fT dL
)

2
1
2
IJIJKK IJ
I
KJL JKIL ILJKJLIK
IJ KLKL
dTdE dE
dE




  
(24)
4. Moving Least-Squares Approximation
In the finite element method, the coupling between the
nodes is accomplished through the use of shape functions,
defined locally over each element, which interpolate the
solution field from nodal values. For a meshless method,
the absence of elements excludes the use of such shape
functions and therefore, and a different local approxima-
tion scheme based on nodal values but independent of
any elements needs to be devised. In this work, we have
chosen to exploit the non-interpolative moving least-
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
345
squares approximation (MLSA) scheme because of its
high accuracy and the ease with which it can be extended
to n-dimensional problems.
Consider a reference domain 0
that contains n
nodes:
T
12
T
123
, in 2-D
, 1
,, in 3-D
aa
a
aaa
YY
Yan
YYY





(25)
Following [29], we define a support domain for node
a
Y, which is a sphere (3D) or disk (2D) centered on
a
Y with a radius a
w
l. A weight function a
w is a
continuous function that is positive in the support domain
and zero outside, i.e.


a
a
0if
0if
aa
w
aa
w
wl
wl


XXY
XXY
(26)
As introduced previously, the sub-d oma in 0
a
for
node a
Y, located entirely inside 0
, is a sphere or
part of a sphere centered on a
Y with a radius a
s
h.
Figure 5 illustrates the meaning of local sub-domain and
support domain.
Two other frequently used concepts are the domain of
definition and the domain of influence. The domain of
definition of a point X is the set of all nodes whose
weight functions are non-zero at X, while the domain of
influence of a node a
Y is the set of all nodes whose
weight functions are non-zero in some part of the sub-
domain of node a
Y. The domain of definition and the
domain of influence are convenient terms in the descrip-
tion of MLSA and local boundary integrals, and are il-
lustrated schematically in Figure 6.
The moving least-squares approximant h
u to a func-
tion u is defined by:
Figure 5. Schematic diagram illustrating the meaning of
local sub-domain and support domain for node a
Y and
node b
Y.
Figure 6. Schematic diagram showing the domain of influence
for node a
Y and the domain of definition for point X.
T
0
(),
h
uXp XcXX (27)
The two vectors p and c are both functions of the spa-
tial coordinates:

T
12
,XXX in 2D or

T
123
,,XXXX in 3D. p is a complete monomial
basis of m terms (e.g., in 2D, m = 3 for a linear basis, and
6 for a quadratic basis), c is a coefficient vector which is
determined by minimizing a weighted discrete L2-norm:
 


2
T
1
ˆ
x
N
aaa
a
w




JXXp YcXu (28)
where ˆa
u is the fictitious nodal displacement that ap-
proximates the value of u at node a
Y, and the upper
limit of summation,
x
N, is the total number of nodes in
the domain of definition of point X. The matrices P, W
and ˆ
u are defined by


1
T
2
T
Tx
x
N
Nm









pY
pY
P
pY
(29)



10
0x
xx
N
NN
w
w






 
X
Wx
X
(30)
1
2
ˆ
ˆ
ˆ
ˆx
N
u
u
u
u
(31)
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
346
Minimization of (28) leads to:
 
T
1
ˆˆ
x
N
aa
h
a

uX Φ
X
uXu (32)
where:
 
TT1
Φ
X
pXA XBX (33)
 
1
1
x
N
a
bba
b
p


XXAXBX (34)
 


T
T
1
x
N
aaa
a
w

AXPWX PBX P
XpYp Y
(35)







T
112 2
xx
NN
ww
w
BX
PW X
XpYXpY
XpY
(36)
The MLSA for a function exists only when
AX is
non-singular. A necessary condition for a well-defined
MLSA is that for each sample point 0
X (a node or
a quadrature point), at least m weight functions are non-
zero and the nodes in the domain of definition of X are
not arranged in a degenerate pattern (such as on a straight
line). In MLSA, the shape function related to node a
Y
is

a
X
. The size of the support domain should be
large enough to ensure the coupling between a minimum
set of nodes, but small enough to capture local variations.
The influence of the choices of the basis functions and
the weight functions on the behavior and the quality of
the shape function has been discussed in [1]. In this work,
we use linear, quadratic monomial basis, and spline
weight functions, defined as follows:

234
1683,0
0,
a
aaa
aa
w
aaa
www
aa
w
w
rrrrl
lll
rl

 



X
(37)
here a
r is the distance from node a
Y to point X,
and a
w
l is the size of support domain. The weight func-
tion has only one parameter, the size of support domain
a
w
l, which makes its use simple. It is noted that MLSA
is non-interpolative, and there is a difference between the
nodal value of the MLSA approximant h
u and the ficti-
tious nodal displacement ˆa
u. For brevity, the subscript
h in h
u will be omitted in the remainder of this paper.
5. Meshless Implementation
We now apply the MLSA to the integral Equation (14) to
establish the meshless implementation. The shape func-
tion, as we have defined it, gives:
 
1
ˆ
x
N
bb
J
J
b
uu
XX (38)

,,
1
ˆ
x
N
bb
J
KKJ
b
uu
XX (39)
 
,,
1
ˆ
x
N
bb
J
KKJ
b
uu
 
XX (40)
where
x
N is the total number of nodes in the domain of
definition of point X.
The related traction term
J
T is

JIJ I
Tn
X
X (41)
where
12
,nn is the normal to the plane passing X
over which the traction acts. For a node b
Y, we define
N and b
B
matrices as:
12
21
,1
,2
,2 ,1
0;
0
0
0
b
bb
bb
nn
nn










N
B
(42)
Combining (42) with (39) and (41), we can express the
traction in terms of the shape functions as follows:



11
1
22
ˆ
ˆ
xb
N
ba
ep
b
b
Tu
Tu







 

XNC BY
X (43)
here ep
C is the elastoplastic stiffness matrix.
With the above discretization and the boundary condi-
tions that
J
J
uu
on a
u
and
J
J
TT on a
t
,
Equation (14) becomes (there is a summation on b and J
but not on a and I):

,,
11
ˆˆ
yy
NN
ab babb
IJJIJ J
bb
aa a
IJ JI
H
uLu
uG



(44)
where ,ab
IJ
H, ,ab
IJ
L and a
I
G are:
 

 

 
,,dΓ
,dΓ
,dΓ
a
u
a
a
t
aba b
ep
IJ IKKJ
ab
IJ
C
ab
IJ
Hu
t
t

 
X
YNCBX X
XYXX
XYXX
(45)


,,dΓ
a
t
aba ba
IJ IJ
Lt
X
YYX (46)
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
347




 







 

 




0
0
0
0
0
0
,,
,,,
0
,
,
,d
,d
, d
,d
,d
,d
,
a
a
t
a
a
a
a
aa
IJIJ
a
IJ J
acor
IJKLKRJRLRL
ae
IJKRJR LLKMNMN
a
IJLKJ LK
a
IJLKJLJ LK
a
IJ J
Gbu
uT
uT u
uuCu
uTun
uT un
tu










XXY X
XYX X
XY XX
XY X
XYXX X
X
YX XX
XY X



d
a
u
a
J
u
YX
(47)
and the upper limit of summation, y
N, is the total num-
ber of nodes in the domain of influence of node a
Y.
The third to sixth terms of Equation (47) are the nonlin-
ear terms because of the large deformation and elasto-
plasticity. In all numerical examples tested in this work,
the body force term is zero.
Owing to the subtraction technique, the singularity in
the third integral on the right hand side of (45) as X ap-
proaches a
Y is cancelled by the term in (46), and
similarly the seventh integral on the right hand side of
(47) is also regularized. Even though the subtraction
technique removes the strong singularity, the integrands
in the first integral on the right hand side of (45) and the
second integral on the right hand side of (47) still contain
the weakly singular ln (r) term. The logarithmic singular-
ity is integrable, but the accuracy of ordinary Legendre-
Gauss integration is poor. We found that the special inte-
gration scheme for the logarithmic singularity [22], which
is reproduced in Table 1 for completeness, achieves ex-
cellent numerical accuracy [1]. The remaining integrals
of Equation (45) and Equation (47) are regular for which
standard quadrature can be used with good accuracy.
In our numerical examples, the numbers of integration
points were as follows: 8 integration points for any inte-
gral along a straight line, and 8 × 8 integration points for
any integral over a sub-domain. For regularized integrals,
the usual Legendre-Gauss integration was used. For inte-
grals containing logarithmic singularity, the special inte-
gration of 8 integration points [32] as listed in Table 1
was used.
6. Enforcement of Boundary Conditions
Appropriate boundary conditions need to be enforced in
order to solve the simultaneous Equations (44). In mesh-
Table 1. Special numerical integration for functions con-
taining logarithmic singularity (8 integration points).
 
1
1
0
1
ln d
N
ii
i
f
rrwfr
r



Abscissas (ri) Weights (wi)
0.0133 2024 4160 8925 0.1644 1660 4728 0030
0.0797 5042 9013 8949 0.2375 2561 0023 3060
0.1978 7102 9326 1880 0.2268 4198 4431 9190
0.3541 5399 4351 9090 0.1757 5407 9006 0700
0.5294 5857 5234 9170 0.1129 2403 0246 7590
0.7018 1452 9939 1000 0.0578 7221 0717 7821
0.8493 7932 0441 1070 0.0209 7907 3742 1330
0.9533 2645 0056 3600 0.0036 8640 7104 0276
less methods, enforcing the essential (Dirichlet) bound-
ary conditions is not as easy as in the finite element
method. Because MLSA is non-interpolative, the essen-
tial boundary condition cannot be accurately enforced by
prescribing values for the fictitious nodal displacements
(i.e., ˆaa
I
I
uu). A number of techniques for enforcing
essential boundary conditions have been developed, in-
cluding the collocation methods [33], Lagrange multi-
plier method [34], penalty method [35,36], Nitsche’s
method [37], coupled meshless-finite element method
[38], and other methods [39-43]. The advantages and
disadvantages of a variety of methods for enforcing es-
sential boundary conditions have been discussed briefly
in [1].
A number of collocation methods have been developed.
The direct collocation method [33] used the condition
ˆaa
I
I
uu (48)
to replace the row of the discretized weak form equation
corresponding to the degree of freedom with prescribed
displacement a
I
u. This is actually inconsistent with the
assumption of MLSA since the fictitious nodal displace-
ment ˆa
I
u is generally not equal to the approximated
displacement value.
A generalized collocation method uses

1
ˆ
n
ababa
I
II
b
uYuu

(49)
as the collocation condition which was shown to yield
more accurate results [44]. Wagner and Liu [45] devel-
oped a corrected collocation method which restores the
consistency of the weak form and enhances convergence.
Wu and Plesha [46] proposed a boundary flux colloca-
tion method to enforce the boundary conditions exactly.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
348
Generally, there are two types of discretization in
meshless methods: 1) Galerkin based method over the
global domain, which is used in EFG [33,42,47,48],
clouds [41], RKPM [49-52]; and 2) loca l weak for m over
multiple local domains, which is employed in [31,53-55].
For Galerkin based method over the global domain, the n
system equations are obtained by applying the weak form
over the global domain once, hence all equations must
hold simultaneously in order to maintain consistency of
the weak form. Replacing a row in the matrix equation
by (49), which contains a linear combination of DOFs
rather than dictating the single value of a constrained
DOF, sacrifices the consistency of the weak form and
compromises the solution accuracy.
For local weak form over multiple local domains, each
equation is obtained by applying the weak form over a
particular local domain, and the weak form needs to be
applied n times for a problem with n DOFs. Conse-
quently, the generalized collocation method with (49)
can be directly used to enforce essential boundary condi-
tions. Because each system equation is independent of
the rest, replacing the equation corresponding to a con-
strained DOF by (49) will not cause any inconsistency in
the weak form.
As the current meshless integral method uses local
weak form over multiple sub-domains, the generalized
collocation method (49) is applicable. For a DOF with
essential boundary condition, we simply use the condi-
tion
I
I
uu rather than applying the integral Equation
(14). In numerical implementation, this means replacing
the governing equation corresponding to the prescribed
DOF with generalized collocation condition (49). Note
that the governing Equation (14) holds for interior nodes
as well as boundary (including corner) nodes, therefore
Equation (49) is applied to nodes on smooth boundaries
as well as at corners. Our numerical tests show that the
generalized collocation method (49) for enforcing essen-
tial boundary conditions works well for the meshless
integral method.
For the natural boundary condition
I
I
tt, no special
treatment is needed. The prescribed traction is directly
used in the second integral in Equation (47).
Although theoretically essential boundary conditions
can be converted into the equivalent natural boundary
conditions and the solution remains unchanged, we ob-
served that the meshless method is always more accurate
when essential boundary conditions are used.
After the boundary conditions are enforced, the gov-
erning equations can be written as

,
1
ˆ
y
N
ab ba
IJ JI
b
K
uR
(50)
where

,,
,when is unconstrained
when
ababa ba
IJIJ IJ
a
ab I
IJ ba
IJ
aa
II
HL
u
K
uu



Y
Y
(51)
when is unconstrained
when
aa
aII
Iaaa
III
Gu
R
uuu
(52)
and the upper limit of summation, y
N, is the total num-
ber of nodes in the domain of influence of node a
Y.
Detailed expressions of various terms are given in Equa-
tions (45) to (47).
7. Solution Algorithm for Elastoplasticity
with Large Deformation
In this work, we will use the fixed point iteration to solve
the governing equations which has the following advan-
tages: the derivative of the stiffness matrix is not needed,
and the implementation is relatively easy. With this algo-
rithm, for each load increment, the stiffness matrix is
computed only once in the first iteration.
The material considered in this study is rate inde-
pendent, and therefore the actual loading rate has no ef-
fect on the solution. It is, however, convenient to intro-
duce a time parameter to describe the loading history. A
piecewise linear function is used in the program to de-
scribe the history of applied loads and is implemented in
a table form in the input data. Each linear segment is
referred to as a load step. For a loading history with N
load steps, the corresponding data block takes the fol-
lowing form:
00
11
22
,
,
,
,
N
N
N
tP
tP
tP
tP

where n indicates the number of load steps. 01
,,,
N
tt t
are the times, while 01
,, ,
N
PP P are the corresponding
load levels. Both 0
t and 0
P are set to zero if the initial
state is load free.
Since the elastoplastic integral equations are strongly
nonlinear, loads have to be applied in small increments,
and iterations are needed to achieve convergence within
each load increment. For the first load step 1
P at a level
sufficient to cause yielding in the material, the solution is
first obtained by assuming elastic behavior. The load and
solution are then scaled linearly to initial yielding by a
scaling factor r where 1r
. The remaining portion of
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
349
load from 1
rP to

1
1rP is then divided into M in-
crements. Within each load increment, the solution is
iterated until convergence. For subsequent load steps, the
load is incremented in a similar manner.
Throughout this section, we use


m
var i to denote a
variable
var (which may be displacement, stress or
strain) for m-th load increment at i-th iteration, and use

var m to denote the converged solution for m-th load
increment.
7.1. Solution Algorithm
1) Divide the load into N load steps, set the load in-
crement number M for each load step and maximum it-
eration number Imax; initialize the load step index n,
load increment index m, and iteration index i to 1.
2) For load step index n and load increment index m,
the convergent state at (m – 1)-th load increment will be
used as the reference configuration. Set the initial dis-
placement relative to the reference configuration



00
m
U, and compute the stiffness matrix [K] based
on the reference configuration. Set the second Piola-
Kirchhoff stress



0
T
, where

is the con-
vergent Cauchy stress at (m – 1)-th load increment.
3) For iteration index i, use Equations (53)-(55) to up-
date the solutions.





 



11 1ii mi
mmm
FKULDEPU
 
 (53)





1ii
m
KUR F
 (54)








1iii
mm
UU U
 (55)
where
ΔR is the load increment from (m – 1)-th load
increment to m-th load increment.


1i
m
LDEPU
are
the nonlinear terms because of large deformation and
elastoplasticity corresponding to the third to sixth terms
of the Equation (47) for m-th load increment at (i – 1)-th
iteration.
4) The following convergence criterion is used to ter-
minate the iteration for each load step:




12
2
i
f
m
RF R
 (56)
f
is a predefined tolerance, usually in the order of
47
10~ 10

.
a) If the program converges in this iteration, update the
geometry, the displacement field, and the second Piola-
Kirchhoff stress for each Gaussian point (Section 7.2)
and obtain the corresponding Cauchy stress by Equation
(57):



0
,,IJIKI KKLKLJLJ L
uT TuJ
 
 (57)
where J is the determinant of the deformation gradient F
which is defined as follows:
x
x
X
Y
yy
X
Y
F (58)
then go to step 5.
b) Otherwise, increase the iteration index by one and
check if i is greater than Imax. If yes, the program cannot
converge for the preset Imax, exit program execution;
otherwise, compute the second Piola-Kirchhoff stresses
at all Gaussian points for all nodes and go to step 3.
5) Increase the load increment index m by one and
check if m is greater than load increment number M. If
yes, exit the program; otherwise, go to step 2. The inte-
grals for next load increment are to be performed over
the deformed geometry of converged solution of current
load increment, and the converged current configuration
will be the reference configuration for the next load in-
crement, hence the updated Lagrangian formulation.
As shown in the previous section, the nonlinear terms


1i
m
LDEPU
(the third to sixth terms of the Equation
(47)) for m-th load increment at (i – 1)-th iteration are
needed in order to calculate the internal force


1i
m
F
.
The integration is usually performed using Gaussian nu-
merical integration. Therefore, the stress state, along with
the stress correction term, is computed at all Gaussian
points in each iteration step. The following section de-
scribes the algorithm to compute the stress at each Gaus-
sian point.
7.2. Procedure for Computing Stresses at Each
Gaussian Point
1) After the solution is converged for (m – 1)-th load
increment, the corresponding Cauchy stress using the
second Piola-Kirchhoff stresses which are basic stress
measure obtained from solution is computed. The Cauchy
stress obtained is then used as the initial value of the
second Piola-Kirchhoff stress



0
m
T for the m-th load
increment. Given





1
i
mm
uU U
 , the strain
increment
E
and the stress increment
e
T are
initially predicted assuming elastic behavior using (59)
and (60) shown below.
and
e
T
, needed to
compute the fourth term of Equation (47), are predicted
by Equations (61) and (62) with


cor
m
T and

e
m
T
set to zero. The parameter EPF (elastoplastic flag) indi-
cates the stress state at a Gaussian point under considera-
tion with EPF = 0 being elastic and EPF = 1 being elas-
toplastic. Initial values of EPF for all Gaussian points are
set to zero.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
350

,,,,
1
2
I
JIJJIKIKJ
Euuuu (59)
ee
TCE

 
 (60)

,,
1
2
I
JIJJI
uu
 (61)


ee
TC

 
 (62)
2) Determine the loading state
2.1) If EPF = 1, the Gaussian point is in an elastoplas-
tic state in previous load step. Compute the loading crite-
rion function L using Equation (18). Use parameter r to
denote the scaling factor such that





,0
e
m
m
fTr T
 
(63)
If L > 0, let r = 0, plastic loading occurs; if L 0, let r
= 1, EPF = 0, compute the yield function after the trial
stress increment is applied




,
e
m
m
fT T

if 0f, let r = 1 and EPF = 0, the point remains in the
elastic state; if 0f, let EPF = 1, the point enters into
elastoplastic state, determine scaling factor r using (63).
Update the stress at this point by
 


e
m
TT rT
(64)
2.2) If EPF = 0, the Gaussian point is in an elastic
state in the previous load step. Compute the yield func-
tion after the trial stress increment is applied





,
e
m
m
fT T

if 0f, let r = 1 and EPF = 0, the point remains in the
elastic state; if 0f, let EPF = 1, the point enters into
elastoplastic state. Determine scaling factor r using (63).
Update the stress at this Gaussian point using Equation
(64).
3) Compute the sub-increment of strain
E:

(1 )rE
EN

(65)
where N is an integer. Integrate numerically to compute
sub-increment of stress
ij
T with n looping from 1 to
N.


n
T is used to denote the second Piola-Kirchhoff
stress increment at the end of n-th iteration and let



00
cor
TT
:

2
2
9
23
1
2
klkk ij
ijijkkij
pe
I
KJL JKIL ILJK JLIK
IJ KLKL
TET
TEEH
E

 





(66)




1nn
TT T
 
(67)
4) Update the variables:

 
N
core e
TTrTT

(68)






i
cor corcor
mm
TT T
(69)







0iN
e
mm
TT TrT (70)






i
ee e
mm
TT T


(71)
7.3. The Computation of Nonlinear Terms
For the third term in Equation (47),
,,a
IJ K
u
XY
needs to be computed, which is the derivative of
,a
IJ
u
XY with respect to
K
X
. Equation (69) is used
to calculate
cor
T for each Gaussian point, and the
third term in Equation (47) is computed accordingly.
Equation (71) is then used to calculate
e
T which is
equal to ,
e
L
KMNMN
CU and the fourth term of Equation
(47) is computed accordingly.
As for the fifth and sixth terms of Equation (47), 0
L
K
T
is the initial value of the second Piola-Kirchhoff stress of
the current load increment.
L
K
T is the increment of the
second Piola-Kirchhoff stress of the current iteration with
respect to 0
L
K
T. These two terms can be computed easily.
8. Numerical Examples
This section presents the numerical solutions to several
large deformation elastoplastic problems using the mesh-
less integral method. The tests include the uniaxial ten-
sion tests, shear tests, rotation tests, and punch test. For
large deformation elasticity cases, closed form solutions
are used for comparison. For large deformation elasto-
plasticity cases, the finite element results obtained using
commercial software, ANSYS (http://www.ansys.com/),
is used as the basis for comparison. Since FEM results
are approximate solutions, convergence study was car-
ried out for each FEM model in which the mesh was re-
fined progressively until the global L2-norm between two
successive meshes was within at least 4
10 . With such
tight convergence, we regarded the FEM solutions as
practically exact. The FEM models used for comparison,
therefore, are much finer than the corresponding mesh-
less models, but care was taken to ensure that all nodes in
a meshless model coincide exactly with a subset of nodes
in the corresponding FEM model in order to facilitate
direct comparison.
The guidelines for the selection of the size of sub-
domain a
s
h and the size of support domain a
w
l have
been discussed in detail in [1]. To summarize briefly, the
sizes are set to be proportional to the minimum nodal
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
351
distance for a node, minDist, defined as the minimum
distance between the node and it neighbors, with propor-
tional coefficients being of the order of 1. A global error
indicator, the L2-norm error in displacement, defined by
L2-norm error =




1
22
1122
1
22
12
1
2
Nnum exnum ex
j
jjj
Nex e
j
jx
j
j
uuuu
uu







is used as a measure of the overall performance of the
numerical method, where num denotes the numerical
result, and ex denotes the exact solution.
In all numerical examples presented below, for elastic
material response the properties of AISI 1020 steel with
Young’s modulus E = 203 GPa and Poisson’s ratio
0.3
under plain strain condition were used. If the
material behavior was elastoplastic, then after yielding,
bilinear stress-strain behavior was assumed with a yield
stress of 260 MPa, and an elastoplastic tangent modulus
of 1 GPa.
8.1. Uniaxial Tension Tests
The uniaxial tension patch test geometry is a 1 × 1 m2
square plate shown in Figure 7(a). Three meshless mod-
els, shown in Figures 7(b) to (d), were simulated in
which spline weight function was used. The left edge
was constrained from moving in X1 direction and traction
free in X2 direction; the bottom edge was constrained
from moving in X2 direction and traction free in X1 direc-
tion; the top edge was traction free in both X1 and X2 di-
rections; and the right edge was traction free in X2 direc-
tion and had various prescribed displacement in X1 de-
pending on the test case. The uniaxial tension patch tests
were investigated for both elastic only material and elas-
toplastic material.
For the first elastic only case, a prescribed displace-
ment U1 of 0.1 m in X1 direction was applied to the right
edgeand the effect of load increments was investigated.
For one and two load increments, closed form solutions
were obtained for comparison. Meshless results for the
same load increments were identical with closed form
results, as shown in Tables 2 and 3. For larger load in-
crements, closed form solutions were quite laborious and
therefore not pursued. Meshless results for up to 12 load
increments are plotted in Figures 8(a) to (c). The figures
show variations of 11
, 33
, and U2, respectively, of
upper edge (denoted as b) with the number of load in-
crements which indicate clearly trend of convergence.
The larger the number of the load increments is, the more
accurate the solutions are. Meshless results for 12 load
increments were σ11 = 21428 MPa, σ33 = 6430 MPa, and
b = –0.04047 m. The unrealistically high stresses are a
(a) A square plate (b) 9 regular nodes
(c) 25 regular nodes (d) 25 irregular nodes
Figure 7. (a) A square plate for the patch test; (b) Meshless
model with 9 regular nodes; (c) Meshless model with 25
regular nodes; (d) Meshless model with 25 irregular nodes.
Table 2. The comparison between hand calculation solution
and meshless result for U1 = 0.1 (1 load increment).
σ11 σ33 b
Closed form solution 23,423 7026 –0.046061
Meshless 23,423 7026 –0.046061
Table 3. The comparison between hand calculation solution
and meshless result for U1 = 0.1 (2 load increments).
σ11 σ33 b
Closed form solution 22,297 6696 –0.042859
Meshless 22,297 6696 –0.042859
result of the large prescribed displacement corresponding
to a 10% engineering strain and the enforced elastic only
behavior. For comparison, FEM results for the same load
increments are σ11 = 21261 MPa, σ33 = 6378 MPa, and b
= –0.04027 m. From this example, a guideline for deter-
mining load increment size was adopted for subsequent
simulations: the largest engineering strain increment in a
load increment should be in the range of 0.005 to 0.01.
For the second elastic only case, a prescribed dis-
placement U1 of 0.5 m in X1 direction was applied to the
right edge using 20, 40, and 60 load increments. Differ-
ence in results using 40 and 60 load increments was
small. Figures 9 and 10 show the variation of σ11 and σ33,
respectively, with applied displacement U1 from the
meshless model using 60 load increments. Corresponding
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
352
21000
21500
22000
22500
23000
23500
24000
0510 15
Number of load increments
σ
11
(MPa)
Meshless
Closed form
(a)
6400
6500
6600
6700
6800
6900
7000
7100
05101
5
Number of load increments
σ
33
(MPa)
Meshless
Closed form
(b)
-0.047
-0.046
-0.045
-0.044
-0.043
-0.042
-0.041
-0.04
-0.039
0510 15
Number of load increments
U
2
of upper edge (m)
Meshless
Closed form
(c)
Figure 8. (a) The convergence test of meshless method for
Cauchy stress σ11 for elastic case; (b) The convergence test
of meshless method for Cauchy stress σ33 for elastic case; (c)
The convergence test of meshless method for displacement
of upper edge for elastic case.
FEM simulation was carried out for comparison. Good
agreement between the meshless results and FEM results
is evident. Nonlinear behavior, although not strong, is
0
10000
20000
30000
40000
50000
60000
70000
80000
90000
100000
00.1 0.2 0.3 0.40.5 0.6
U
1
(m)
σ
11
(MPa)
FEM
Meshless
Figure 9. FEM results versus meshless results of Cauchy
stress σ11 for uniaxial tension simulation for elastic case.
0
5000
10000
15000
20000
25000
30000
00.1 0.2 0.3 0.4 0.5 0.6
U
1
(m)
σ
33
(MPa)
FEM
Meshless
Figure 10. FEM results versus meshless results of Cauchy
stress σ11 for uniaxial tension simulation for elastic case.
noticeable.
For elastoplastic material behavior, since closed form
results are not available, FEM results were used as the
basis for comparison. For the first elastoplastic case, a
prescribed displacement U1 of 0.1 m in X1 direction was
applied to the right edge under different number of load
increments. Figures 11(a) to (c) show variations of σ11,
σ33, and U2, respectively, of the upper edge (denoted as b)
with the number of load increments, together with the
corresponding FEM results.
The figures indicate trend of convergence, and con-
verged values for FEM and meshless method are close to
each other. The meshlessresults for 12 load increments
are σ11 = 426 MPa, σ33 = 212 MPa, and b = –0.09009 m.
The corresponding results of FEM with 12 load incre-
ments are σ11 = 426 MPa, σ33 = 212 MPa, and b =
–0.08978 m. It is noted that the current meshless method
is based on Green-Naghdi’s theory, while FEM is based
on E. H. Lee’s theory [2]. The differences between mesh-
less results and FEM results show that large deformation
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
353
300
320
340
360
380
400
420
440
460
480
500
0510 15
number of load increments
σ
11
(MPa)
Meshless
FEM
(a)
100
120
140
160
180
200
220
240
260
280
300
0510 15
number of load increments
σ
33
(MPa)
Meshless
FEM
(b)
-0.16
-0.14
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
0
051015
number of load increments
U
2
of upper edge (m)
Meshless
FEM
(c)
Figure 11. (a) The convergence test of FEM results and
meshless results for Cauchy stress σ11 for elastoplastic case;
(b) The convergence test of FEM results and meshless
results for Cauchy stress σ33 for elastoplastic case; (c) The
convergence test of FEM results and meshless results for
displacement of upper edge for elastoplastic case.
elastoplasticity theories may have some influence on the
numerical results.
For the second elastoplastic case, a prescribed dis-
placement U1 of 0.5 m in X1 direction was applied to the
right edge using 60 load increments. Figures 12 and 13
show the variation of σ11 and σ33, respectively, with ap-
plied displacement U1 from the meshless model and
FEM model. Good agreement between the meshless re-
sults and FEM results is evident.
8.2. The Shear Tests
The shear patch test geometry is a 1 × 1 m2 square plate
shown in Figure 14. Three meshless models, shown in
Figures 14(b) to (d), with spline weight function were
simulated. The motion of the square was described by
112
x
XX
and 22
x
X
. Three values of
, 0.1,
0.5 and 1.0, were simulated. For all three models, pre-
scribed displacements were applied to all edges of the
plate. For the shear tests, both elastic only material and
elastoplastic material were investigated.
0
100
200
300
400
500
600
700
800
900
00.1 0.2 0.3 0.4 0.50.6
U
1
(m)
σ
11
(MPa)
FEM
Meshless
Figure 12. FEM results versus meshless results of σ11 for
uniaxial tension simulation for elastoplastic case.
0
50
100
150
200
250
300
350
400
450
00.1 0.20.30.40.5 0.6
U
1
(m)
σ
33
(MPa)
FEM
Meshless
Figure 13. FEM versus meshless results of Cauchy stress σ33
for uniaxial tension simulation for elastoplastic case.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
354
(a) A square plate (b) 9 regular nodes
(c) 25 regular nodes (d) 25 irregular nodes
Figure 14. (a) A square plate for the shear tests; (b) Meshless
model with 9 regular nodes; (c) Meshless model with 25
regular nodes; (d) Meshless model with 25 irregular nodes.
For elastic only material behavior, analytical results in
terms of the Cauchy stress using Jaumann rate for this
problem are given as [56]:

1122 1cos
 
  (72)
12 sin

(73)
The constitutive equation in terms of the Jaumann rate
is given by:
Ttrace 2
JJJ

 WW
σσ σσ
D
ID
(74)
where D is the rate-of-deformation tensor and 0.1
W is the spin tensor, and the superscript, J, denotes that
the material constants are used with the Jaumann rate.
For the first elastic only shear test with 0.1
, the
L2-norms between meshless results and analytical solu-
tions for the 9 node, 25 node regular, and 25 node ir-
regular models are 3.2 × 10–16, 2.6 × 10–18, and 3.6 ×
10–17 respectively. The stresses for all three models are
σ11 = 390 MPa, σ22 = –391.48 MPa, σ12 = 7794.88 MPa
which are practically identical with the analytical solu-
tion using Jaumann rate: σ11 = 390 MPa, σ22 = –390 MPa,
σ12 = 7794.68 MPa. Note that although it is pure shear
loading, non-zero normal stresses in the X1 and X2 direc-
tions are present due to large deformation.
For the second elastic only shear test with 0.5
,
the L2-norms between meshless results and analytical
solutions for the 9 node, 25 node regular, and 25 node
irregular models are 6.2 × 10–16, 5.5 × 10–18, and 4.7 ×
10–16 respectively. The stresses for all three models are
σ11 = 9560 MPa, σ22 = –9555 MPa, σ12 = 37431 MPa
which are practically identical with the analytical solu-
tion using Jaumann rate: σ11 = 9557 MPa, σ22 =–9557
MPa, σ12 = 37432 MPa.
For the third elastic only shear test, Figures 15(a)-(c)
show the variation of normal stresses σ11, σ22, and shear
stress σ12 with increasing
up to 3.0 from meshless
model, together with the analytical results. These three
figures reveal that the meshless results match the ana-
lytical results almost identically.
(a)
(b)
(c)
Figure 15. (a) Analytical versus meshless results of Cauchy
stress σ11 of the shear test for elastic case; (b) Analytical
versus meshless results of Cauchy stress σ22 of the shear test
for elastic case; (c) Analytical versus meshless results of
Cauchy stress σ12 of the shear test for elastic case.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
355
For elastoplastic material behavior, FEM results were
used as the basis for comparison. For 0.1
, the
L2-norms between the meshless and FEM results for the
9 node, 25 node regular, and 25 node irregular models
are 3.7 × 10–13, 4.6 × 10–13, and 6.4 × 10–14 respectively.
The shear stress for all three models is σ12 = 182.32 MPa
which is in good agreement with the FEM solution of σ12
= 182.60 MPa.
For elastoplastic behavior with 0.5
, the L2-norms
between the meshless and FEM results for the 9 node, 25
node regular, and 25 node irregular models are 5.2 ×
10–13, 1.8 × 10–13, and 9.4 × 10–11 respectively. The shear
stress for all three models is σ12 = 316.52 MPa which is
again in good agreement with the FEM solution of σ12 =
316.89 MPa. Figure 16 shows the variation of shear
stress σ12 with
from the meshless model, together
with FEM results. The meshless results match the FEM
results almost identically.
8.3. The Rigid Body Rotation Tests
The rigid body rotation test geometry is a 1 × 1 m2 square
plate shown in Figure 17(a). Figure 17(b) shows the
current configuration with rotation angle θ. Three mesh-
less models (9 regular nodes, 25 regular nodes, and 25
irregular nodes), with spline weight function were simu-
lated. Properties of AISI 1020 steel with Young’s modulus
E = 203 GPa and Poisson’s ratio 0.3
under plain
strain condition were used, although the material proper-
ties have no effect on the simulation results. For all three
models, prescribed displacements were applied to all
edges of the plate. The plate was subjected to the follow-
ing initial applied stresses 0
11 20MPa
, 0
22 0MPa
,
and 0
12 0MPa
which corotate with the plate.
The simulation was carried out with θ increasing from
0 to π2, and the results were reported in step of
π12
 as follows:
θ σ11 (MPa) σ22 (MPa) σ12 (MPa)
π12 18.66 1.3397 5
π6 15 5 8.66
π4 10 10 10
π3 5 15 8.66
5π12 1.3397 18.66 5
π2 0 20 0
All these results are identical with the analytical re-
sults given by the following equation [56]:
2
0
11
2
1
cossin 2
2
1sin 2sin
2






σ (75)
For the three meshless models with different nodal
0
50
100
150
200
250
300
350
00.1 0.2 0.3 0.4 0.5 0.6
γ
σ
12
(MPa)
FEM
Meshless
Figure 16. FEM results versus meshless results of the shear
test for elastoplastic case.
(a) Original configuration (b) Current configuration
Figure 17. Rotation of a prestressed square with no defor-
mation. (a) Original configuration; (b) Current configura-
tion.
densities, virtually identical results were obtained. Fig-
ures 18-20 show meshless and analytical results of σ11,
σ12, and σ22, respectively, versus different rotation angle θ
for the rotation tests. The figures reveal that the meshless
results are practically identical with the analytical results.
8.4. The Punch Test
The punch test example is illustrated in Figure 21. The
model geometry had 561 nodes, and the spline weight
function and linear monomial basis were used. Both the
left edge and the right edge were constrained from mov-
ing in X1 direction and were traction free in X2 direction.
The bottom edge was constrained from moving in X2
direction and was traction free in X1 direction. To simu-
late the compression from a punch on top, uniform pre-
scribed displacement was applied on the left half of the
top edgein the X2 direction, and the traction free condi-
tion was enforcedon the left half of the top edgein the X1
direction. The right half of the top edge was traction free.
For the punch test, elastoplastic material was investi-
gated.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
356
0
5
10
15
20
25
00.511.52
θ
σ
11
(MPa)
Analytical
Meshless
Figure 18. Analytical versus meshless results of Cauchy
stress σ11 for rotation test for different rotation angle.
0
2
4
6
8
10
12
00.511.52
θ
σ
12
(MPa)
Analytical
Meshless
Figure 19. Analytical versus meshless results of Cauchy
stress σ12 of the rotation test for different rotatation angle.
0
5
10
15
20
25
00.511.52
θ
σ
22
(MPa)
Analytical
Meshless
Figure 20. Analytical versus meshless results of Cauchy
stress σ22 of the rotation test for different rotation angle.
A prescribed displacement U2 of –0.05 m was first ap-
plied to the left half of the top edge. Figure 22 shows the
Figure 21. The square plate for punch test.
Figure 22. The undeformed model (solid diamonds), de-
formed meshless model (solid squares), and FEM model
(triangles) for the punch test with U2 = –0.05 m.
undeformedmeshless model (solid diamonds), deformed
meshless model (solid squares) and FEM model (empty
triangles). The meshless results match FEM results very
well and the L2-norm between meshless and FEM results
is 0.0029. Figures 23 and 24 present the distribution of
σ22 and von Mises stress, respectively, along X2 = 0, in-
dicatingthat the meshless results and FEM solution are
comparable with each other.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
357
850
750
650
550
450
350
250
00.2 0.4 0.6 0.81
X
1
σ
22
along X
2
=0 (MPa)
FEM
Meshless
Figure 23. The distribution of Cauchy stress σ22 along X2 =
0 of the 561 node model with U2 = –0.05 m.
150
200
250
300
350
400
00.2 0.4 0.6 0.81
X
1
vonMisesalongX
2
=0(MPa)
FEM
Meshless
Figure 24. The distribution of von Mises stress along X2 = 0
of the 561 node model with U2 = –0.05 m.
Figure 25 shows the undeformed model (solid dia-
monds), deformed meshless model (solid squares), and
FEM model (triangles) when the prescribed displacement
of the left half of the top edge was increased to U2 =
–0.08 m. The L2-norm error between meshless results
and FEM solution is 0.0056. Figures 26 and 27 present
the distribution of σ22 and von Mises stress, respectively,
along X2 = 0. The two figures indicate that the meshless
results are in reasonable agreement with FEM results.
9. Concluding Remarks
In this paper, the large deformation elastoplasticmeshless
integral method based on regularized boundary integral
equation [1] is presented. The updated Lagrangian gov-
erning integral equation is obtained from the weak form
of elastoplasticity based on Green-Naghdi’s theory over a
local sub-domain. Green-Naghdi’s theory starts with the
additive decomposition of the Green-Lagrange strain into
Figure 25. The undeformed model (solid diamonds), de-
formed meshless model (solid squares), and FEM model
(triangles) for the punch test with U2 = –0.08 m.
-1000
-900
-800
-700
-600
-500
-400
-300
00.2 0.4 0.6 0.811.2
X
1
σ
22
along X
2
=0 (MPa)
Meshless
FEM
Figure 26. The distribution of Cauchy stress σ22 along X2 =
0 of the 561 node model with U2 = –0.08 m.
200
250
300
350
400
450
00.2 0.4 0.6 0.81
X
1
vonMisesstressalongX
2
=0(MPa)
Meshless
FEM
Figure 27. The distribution of von Mises stress along X2 = 0
of the 561 node model with Uy = –0.08 m.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
358
the elastic part and plastic part and considers a J2 elasto-
plastic constitutive law that relates the Green-Lagrange
strain to second Piola-Kirchhoff stress. The generalized
collocation method is employed to enforce the essential
boundary conditions exactly, which is simple and com-
putationally efficient. The natural boundary conditions
are incorporated in the system governing equations and
require no special handling. Numerical results show that
this method is accurate and robust.
REFERENCES
[1] A. Bodin, J. Ma, X. J. Xin and P. Krishnaswami, “A
Meshless Integral Method Based on Regularized Bound-
ary Integral Equation,” Computer Methods in Applied
Mechanics and Engineering, Vol. 195, No. 44-47, 2006,
pp. 6258-6286. 0doi:10.1016/j.cma.2005.12.005
[2] A. E. Green and P. M. Naghdi, “A General Theory of an
Elasto-Plastic Continuum,” Archive for Rational Me-
chanics and Analysis, Vol. 18, No. 4, 1965, pp. 251-281.
1doi:10.1007/BF00251666
[3] J. H. Chiou, J. D. Lee and A. G. Erdman, “Comparison
between Two Theories of Plasticity,” Computers & Struc-
tures, Vol. 24, No. 1, 1986, pp. 23-37.
2doi:10.1016/0045-7949(86)90332-9
[4] E. H. Lee, “Elastic-Plastic Deformation at Finite Strains,”
Journal of Applied Mechanics, Vol. 36, No. 1, 1969, pp.
1-6. 3doi:10.1115/1.3564580
[5] J. H. Chiou, J. D. Lee and A. G. Erdman, “Development
of a Three-Dimensional Finite Element Program for
Large Strain Elastic-Plastic Solids,” Computers & Struc-
tures, Vol. 36, No. 4, 1990, pp. 631-645.
4doi:10.1016/0045-7949(90)90078-G
[6] J. D. Lee, “A Large-Strain Elastic-Plastic Finite Element
Analysis of Rolling Process,” Computer Methods in Ap-
plied Mechanics and Engineering, Vol. 161, No. 3-4,
1998, pp. 315-347. 5doi:10.1016/S0045-7825(97)00324-1
[7] P. Hu, “Finite-Element Numerical Analysis of Sheet
Metal under Unaxial Tension with a New Yield Crite-
rion,” Journal of Materials Processing Technology, Vol.
31, No. 1-2, 1992, pp. 245-253.
6doi:10.1016/0924-0136(92)90025-N
[8] T. Belytschko, P. Krysl and Y. Krongauz, “A Three-Di-
mensional Explicit Element-Free Galerkin Method,” In-
ternational Journal for Numerical methods in Fluids, Vol.
24, No. 12, 1997, pp. 1253-1270.
7doi:10.1002/(SICI)1097-0363(199706)24:12<1253::AID-
FLD558>3.0.CO;2-Z
[9] R. Rossi and M. K. Alves, “On the Analysis of an EFG
Method under Large Deformations and Volumetric Lock-
ing,” Computational Mechanics, Vol. 39, No. 4, 2007, pp.
381-399. 8doi:10.1007/s00466-006-0035-z
[10] Y. P. Chen, A. Eskandarian, M. Oskard and J. D. Lee,
“Meshless Analysis of High-Speed Impact,” Theoretical
and Applied Fracture Mechanics, Vol. 44, No. 3, 2005,
pp. 201-207. 9doi:10.1016/j.tafmec.2005.09.007
[11] A. Eskandarian, Y. P. Chen, M. Oskard and J. D. Lee,
“Meshless Analysis of Fracture. Plasticity and Impact,”
Proceedings of ASME 2003 International Mechanical
Engineering Congress and Exposition, Washington DC,
15-21 November 2003, pp. 89-97.
[12] Y. Xiong, H. Cui and S. Long, “Meshless local Petrov-
Galerkin Method for Large Deformation Analysis,” Chi-
nese Journal of Computational Mechanics, Vol. 26, No. 3,
2009, pp. 353-357.
[13] D. Hu, S. Long, X. Han and G. Li, “A Meshless Local
Petrov-Galerkin Method for Large Deformation Contact
Analysis of Elastomers,” Engineering Analysis with Bou-
ndary Elements, Vol. 31, No. 7, 2007, pp. 657-666.
1doi:10.1016/j.enganabound.2006.11.005
[14] Z. Han, A. Rajendran and S. Atluri, “Meshless Local
Petrov-Galerkin (MLPG) Approaches for Solving Nonlin-
ear Problems with Large Deformations and Rotations,”
Computer Modeling in Engineering and Sciences, Vol. 10,
No. 1, 2005, pp 1-12.
[15] D. Li, Z. Lu and W. Kang, “A Coupled Finite Element
and Meshless Local Petrov-Galerkin Method for Large
Deformation Problems,” Advanced Materials Research,
Vol. 97-101, 2010, pp. 3777-3780.
1doi:10.4028/www.scientific.net/AMR.97-101.3777
[16] Y. Gu, Q. Wang and K. Lam, “A Meshless Local Kriging
Method for Large Deformation Analyses,” Computer
Methods in Applied Mechanics and Engineering, Vol.
196, No. 9-12, 2007, pp. 1673-1684.
1doi:10.1016/j.cma.2006.09.017
[17] Y. Gu, C. Yan and P. Yarlagadda, “An Advanced Mesh-
less Technique for Large Deformation Analysis of Metal
Forming,” Australian Journal of Mechanical Engineering,
Vol. 7, No. 1, 2009, pp. 25-32.
[18] Y. Gu, “Meshless TL and UL Approaches for Large De-
formation Analysis,” Advanced Materials Research, Vol.
139-141, 2010, pp. 893-896.
1doi:10.4028/www.scientific.net/AMR.139-141.893
[19] H. Gun, S. Caliskan and A. Gun, “A Meshless Formula-
tion of Euler-Bernoulli Beam Theory for Prediction of
Large Deformation,” Textile Research Journal, Vol. 81,
No. 10, 2011, pp. 1075-1080.
1doi:10.1177/0040517511398946
[20] Q. Li and K. Lee, “An Adaptive Meshless Method for
Analyzing Large Mechanical Deformation and Contacts,”
Journal of Applied Mechanics, Transactions ASME, Vol.
75, No. 4, 2008, Article ID: 041014.
1doi:10.1115/1.2912938
[21] Q. Li and K. Lee, “An Adaptive Meshless Method for
Modeling Large Mechanical Deformation and Contacts,”
IEEE International Conference on Robotics and Automa-
tion, Roma, 10-14 April 2007, pp. 1207-1212.
[22] H. Zhu, W. Liu, Y. Cai and Y. Miao, “A Meshless Local
Natural Neighbor Interpolation Method for Two-Dimen-
sion Incompressible Large Deformation Analysis,” Engi-
neering Analysis with Boundary Elements, Vol. 31, No.
10, 2007, pp. 856-862.
1doi:10.1016/j.enganabound.2007.02.003
[23] S. Wang, “A Large-Deformation Galerkin SPH Method
for Fracture,” Journal of Engineering Mathematics, Vol.
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
359
71, No. 3, 2011, pp. 305-318.
1doi:10.1007/s10665-011-9455-7
[24] J. Chen, C. Pan, C. Wu and W. Liu, “Reproducing Kernel
Particle Methods for Large Deformation Analysis of Non-
Linear Structures,” Computer Methods in Applied Me-
chanics and Engineering, Vol. 139, No. 1-4, 1996, pp.
195-227. 1doi:10.1016/S0045-7825(96)01083-3
[25] S. Jun, W. Liu and T. Belytschko, “Explicit Reproducing
Kernel Particle Methods for Large Deformation Prob-
lems,” International Journal for Numerical Methods in
Engineering, Vol. 41, No. 1, 1998, pp. 137-166.
1doi:10.1002/(SICI)1097-0207(19980115)41:1<137::AID-
NME280>3.0.CO;2-A
[26] K. Liew, T. Ng and Y. Wu, “Meshfree Method for Large
Deformation Analysis—A Reproducing Kernel Particle
Approach,” Engineering Structures, Vol. 24, No. 5, 2002,
pp. 543-551. 2doi:10.1016/S0141-0296(01)00120-1
[27] D. Li, J. Xu and W. Kang, “Applying Element-Free Gal-
erkin Method to Simulate Die Forging Problems,” Advan-
ced Materials Research, Vol. 139-141, 2010, pp. 1174-1177.
2doi:10.4028/www.scientific.net/AMR.139-141.1174
[28] W. Quak, A. van den Boogaard and J. Huétink, “Mesh-
less Methods and Forming Processes,” International
Journal of Material Forming, Vol. 2, No. S1, 2009, pp.
585-588.
[29] J. Ma, X. J. Xin and P. Krishnaswami, “An Elastoplastic
Meshless Integral Method Based on Regularized Boun-
dary Integral Equation,” Computer Methods in Applied
Mechanics and Engineering, Vol. 197, No. 51-52, 2008,
pp. 4774-4788. 2doi:10.1016/j.cma.2008.06.019
[30] J. Ma, “Application of Meshless Integral Method to Metal
Forming,” Proceedings of the ASME Design Engineering
Technical Conference, Montreal, 15-18 August 2010, pp.
141-151.
[31] S. N. Atluri, J. Sladeck, V. Sladeck and T. Zhu, “The
Local Boundary Integral Equation (LBIE) and Its Mesh-
less Implementation for Linear Elasticity,” Computa-
tional Mechanics, Vol. 25, No. 2-3, 2000, pp. 180-198.
2doi:10.1007/s004660050468
[32] A. H. Stroud and D. Secrest, “Gaussian Quadrature For-
mulas,” Prentice-Hall, Upper Saddle River, 1966.
[33] Y. Y. Lu, T. Belytschko and L. Gu, “A New Implementa-
tion of the Element Free Galerkin Method,” Computer
Methods in Applied Mechanics and Engineering, Vol.
113, No. 3-4, 1994, pp. 397-414.
2doi:10.1016/0045-7825(94)90056-6
[34] T. Belytschko, Y. Y. Lu and L. Gu, “Element Free Gal-
erkin Method,” International Journal for Numerical Me-
thods in Engineering, Vol. 37, No. 2, 1994, pp. 229-256.
2doi:10.1002/nme.1620370205
[35] L. Gavete, J. J. Benito, S. Falcon and A. Ruiz, “Imple-
mentation of Essential Boundary Conditions in a Mesh-
less Method,” Communications in Numerical Methods.
Engineering, Vol. 16, No. 6, 2000, pp. 409-421.
2doi:10.1002/1099-0887(200006)16:6<409::AID-CNM34
9>3.0.CO;2-Z
[36] T. Zhu and S. N. Atluri, “A Modified Collocation Method
and a Penalty Foumulation for Enforcing the Essential
Boundary Conditions in the Element Free Galerkin Me-
thod,” Computational Mechanics, Vol. 21, No. 3, 1998,
pp. 211-222. 2doi:10.1007/s004660050296
[37] D. N. Arnold, F. Brezzi, B. Cockburn and L. D. Marini,
“Unified Analysis of Discontinuous Galerkin Methods for
Elliptic Problems,” SIMA Journal on Numerical Analysis,
Vol. 39, No. 5, 2002, pp. 1749-1779.
2doi:10.1137/S0036142901384162
[38] D. Hegen, “Element-Free Galerkin Methods in Combina-
tion with Finite Element Approaches,” Computer Meth-
ods in Applied Mechanics and Engineering, Vol. 19, No.
1, 1996, pp. 120-135.
[39] J. Gosz and W. K. Liu, “Admissible Approximations for
Essential Boundary Conditions in the Reproducing Ker-
nel Particle Method,” Computational Mechanics, Vol. 19,
No. 2, 1996, pp. 120-135. 2doi:10.1007/BF02824850
[40] F. C. Gunther and W. K. Liu, “Implementation of Bound-
ary Conditions for Meshless Methods,” Computer Meth-
ods in Applied Mechanics and Engineering, Vol. 163, No.
1-4, 1998, pp. 205-230.
3doi:10.1016/S0045-7825(98)00014-0
[41] C. A. M. Duarte and J. T. Oden, “An h-p Adaptive Me-
thod Using Clouds,” Computer Methods in Applied Me-
chanics and Engineering, Vol. 139, No. 1-4, 1996, pp.
237-262. 3doi:10.1016/S0045-7825(96)01085-7
[42] Y. Y. Lu, T. Belytschko and M. Tabbara, “Element-Free
Galerkin Method for Wave Propagation and Dynamic
Fracture,” Computer Methods in Applied Mechanics and
Engineering, Vol. 126, No. 1-2, 1995, pp. 131-153.
3doi:10.1016/0045-7825(95)00804-A
[43] X. Zhang, X. Liu, M. W. Lu and Y. Chen, “Imposition of
Essential Boundary Conditions by Displacement Con-
straint Equations in Meshless Methods,” Communications
in Numerical Methods in Engineering, Vol. 17, No. 3,
2001, pp. 165-178. 3doi:10.1002/cnm.395
[44] T. Zhu and S. N. Atluri, “A Modified Collocation Method
and a Penalty Foumulation for Enforcing the Essential
Boundary Conditions in the Element Free Galerkin Me-
thod,” Computational Mechanics, Vol. 21, No. 3, 1998,
pp. 211-222. 3doi:10.1007/s004660050296
[45] G. J. Wagner and W. K. Liu, “Application of Essential
Boundary Conditions in Mesh-Free Methods: A Corre-
cted Collocation Method,” International Journal for Nu-
merical Methods in Engineering, Vol. 47, No. 8, 2000, pp.
1367-1379.
3doi:10.1002/(SICI)1097-0207(20000320)47:8<1367::AID
-NME822>3.0.CO;2-Y
[46] C. C. Wu and M. E. Plesha, “Essential Boundary Condi-
tion Enforcement in Meshless Methods: Boundary Flux
Collocation Method,” International Journal for Numeri-
cal Methods in Engineering, Vol. 53, No. 3, 2002, pp.
499-514. 3doi:10.1002/nme.267
[47] T. Belytschko, Y. Y. Lu and L. Gu, “Element Free Gal-
erkin Method,” International Journal for Numerical Me-
thods in Engineering, Vol. 37, No. 2, 1994, pp. 229-256.
3doi:10.1002/nme.1620370205
[48] P. Krysl and T. Belytschko, “Analysis of Thin Plates by
the Element-Free Galerkin Method,” Computational Me-
J. F. MA, X. J. XIN
Copyright © 2012 SciRes. WJM
360
chanics, Vol. 17, No. 1, 1998, pp. 26-35.
[49] W. K. Liu and Y. Chen, “Wavelet and Multiple Scale
Reproducing Kernel Methods,” International Journal for
Numerical Methods in Fluids, Vol. 21, No. 10, 1995, pp.
901-931.
[50] N. R. Aluru, “A Reproducing Kernel Particle Method for
Meshless Analysis of Microelectromechanical Systems,”
Computational Mechanics, Vol. 23, No. 4, 1999, pp. 324-
338. 3doi:10.1007/s004660050413
[51] J. S. Chen, C. Pan and C. T. Wu, “A Lagrangian Repro-
ducing Kernel Particle Method for Metal Forming Analy-
sis,” Computational Mechanics, Vol. 22, No. 3, 1998, pp.
289-338. 3doi:10.1007/s004660050361
[52] J. S. Chen, C. Pan, C. T. Wu and W. K. Liu, “Reproduc-
ing Kernel Particle Methods for Large Deformation Ana-
lysis of Non-Linear Structures,” Computer Methods in
Applied Mechanics and Engineering, Vol. 139, No. 1-4,
1996, pp. 195-227. doi:10.1016/S0045-7825(96)01083-3
[53] T. Zhu, J.-D. Zhang and S. N. Atluri, “A Local Boundary
Integral Equation (LBIE) Method in Computational Me-
chanics, and a Meshless Discretization Approach,” Com-
putational Mechanics, Vol. 21, No. 3, 1998, pp. 223-235.
4doi:10.1007/s004660050297
[54] J. Sladek, V. Sladeck and R. Van Keer, “Meshless Local
Boundary Integral Equation for 2D Elastodynamic Prob-
lems,” International Journal for Numerical Methods in
Engineering, Vol. 57, No. 2, 2003, pp. 235-249.
4doi:10.1002/nme.675
[55] S. Long and Q. Zhang, “Analysis of Thin Plates by the
Local Boundary Integral Equation (LBIE) Method,” En-
gineering Analysis with Boundary Elements, Vol. 26, No.
8, 2002, pp. 707-718.
4doi:10.1016/S0955-7997(02)00025-5
[56] T. Belytschko, W. K. Liu and B. Moran, “Nonlinear Fi-
nite Elements for Continua and Structures,” John Wiley
& Sons, Ltd., Chichester, 2000.