A simple solution method to 3D integral nonlocal elasticity: Isotropic-BEM coupled with strong form local radial point interpolation

Engineering Analysis with Boundary Elements 36 (2012) 606–612

Contents lists available at SciVerse ScienceDirect

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

A simple solution method to 3D integral nonlocal elasticity: Isotropic-BEM coupled with strong form local radial point interpolation M. Schwartz a,b,n, N.T. Niane b, R. Kouitat Njiwa a a b

Institut Jean Lamour, Dpt SI2M, UMR 7198 CNRS, Ecole des Mines, Parc de Saurupt, CS 14234, Nancy Cedex 54042, France PSA Peugeot-Citro¨ en, Centre Technique La Garenne Colombes, 18 rue des Fauvelles, La Garenne Colombes 92250, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 15 July 2011 Accepted 8 October 2011 Available online 16 December 2011

Nonlocal theories are of growing interest as they can address problems that lead to unphysical results in the framework of classical models. In this work, a solution procedure for three-dimensional integral nonlocal elastic solid is presented. The approach is based on the partition of the displacement field into complementary and particular parts. The complementary displacement is the solution of a Navier type equation and is obtained by the boundary element method, while the particular displacement is obtained using a local radial point interpolation method. The method is illustrated by comparing the responses to some simple loadings of a solid of finite extent with the original nonlocal model of Eringen and the enhanced model of Polizzotto. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Nonlocal elasticity Isotropic-BEM Meshfree strong form Radial basis function

1. Introduction Classical theory of linear elasticity has proven very efficient for many engineering problems. However, it is increasingly clear that there are situations that cannot be satisfactorily addressed by this theory. As well-known examples one can mention the mechanical modeling of carbon nanotubes where size effect become prominent and have to be accounted for (e.g. Sudak [1], Wang and Varadan [2], Filiz and Aydogdu [3] and Hu et al. [4]), the singularity of stress predicted on a crack front. The latter is physically unacceptable and reflects the inability of the classical linear elasticity (local description of the behavior) to deal with problems in domains with sharp geometric singularities. As noted by Eringen [5], this inability is due to the limitations of the classical theory, which fails to handle physical problems when the influence of events arising at a microstructure or microscopic level is significant. According to Eringen [5], the limitation of the classical elasticity is due to its lack of an internal characteristic length scale. The consideration of nonlocal elasticity appeared in the late sixties in the papers by Kroner [6], Kunin [7] and Krumhanls [8]. These authors noted that within the framework of classical elasticity, internal forces in the body are contact forces, which have zero range, whereas cohesive forces in real materials do have a finite range. In their works, an enriched version of elastic

n ¨ Centre Technique La Garenne Corresponding author at: PSA Peugeot-Citroen, Colombes, Case Courrier: LG075, 18 rue des Fauvelles, La Garenne Colombes 92250, France. Tel: þ 156476945 (from France), þ 33156476945 (from abroad); fax: þ 156474059 (from France) or þ 33156474059 (from abroad). E-mail address: [email protected] (M. Schwartz).

0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.10.004

continuum mechanics, which takes into account the long range interactions and the microstructure effects, allows a more refine description of the elastic behavior of the material. Thereafter, improved formulations of nonlocal constitutive relations were proposed by Edelen and Laws [9], Edelen et al. [10] and Eringen and Edelen [11]. In the case of linear isotropic nonlocal elastic solids, much work was done by Eringen and Kim [12] and Eringen et al. [13]. In the last mentioned studies, the nonlocal constitutive relation involves an integral term of the strain field over the entire volume of the solid. As a consequence, the stress field at a point depends on the strain in the whole solid. Following this approach, substantial works have been done to address crack problems (e.g. Eringen [14,15], Gao and Dai [16], Minghao et al. [17]), dislocations problems (e.g. Eringen [18,19], Lazar et al. [20]), problems with concentrated loads (Povstenko and Kubik [21]), contact problems (Artan et al. [22]). With respect to continuum boundary value problem, Polizzotto [23] pointed out a shortcoming of the integral formulation due to Eringen. Indeed, regarding a solid of finite extent, the later does not provide a uniform stress under a uniform strain field. Polizzotto et al. [24] proposed an enhanced constitutive law which solves the shortcoming. The model has already been applied to 2D crack problem by Jackiewicz and Holka [25]. As pointed out by Pisano et al. [26], regarding nonlocal elasticity of integral type, an issue that should be investigated thoroughly is related with the numerical solution of structural problems. These authors investigate the application of the finite element method to a 2D problem on the basis of the nonlocal finite element method proposed by Polizzotto [23]. The present work proposed a solution method of 3D nonlocal elastic problems. The Eringen

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

607

model of nonlocal elasticity and the enhanced version by Polizzotto are considered. The proposed strategy couples the conventional boundary element method (BEM) for elastostatic problems with the local radial point collocation method. Details of the approach are given in Section 2 below. In Section 3, the effectiveness of the method is demonstrated on examples of tension and shear loading of a unit cube and torsion of a cylindrical specimen.

2. Problem statement and solution method 2.1. Problem statement Throughout, indicial notation and associated summation on repeated indices is adopted. Vectors and matrices in sentences are boldface while scalars are italics. Consider a homogeneous elastic solid whose material occupies a geometrical domain O. According to Eringen [5], a simplified relation between the stress and the strain fields is given by Z sij ðxÞ ¼ C ijkl gðx,x0 Þekl ðx0 ÞdVðx0 Þ ð1Þ O

needs to be improved. This has been done by Polizzotto [24] who proposed a stress–strain relation of the type:

sij ðxÞ ¼ C ijkl Rkl ðeÞ

ð4Þ

where

In relation (1), C stands for the fourth order tensor of elastic constants of the local theory. The function g(x,x0 ) is positive and is known as attenuation function. e is the usual small strain tensor, which is defined by the displacement vector u as

ekl ðxÞ ¼ ðuk,l ðxÞ þ ul,k ðxÞÞ=2

Fig. 1. Schematic representation of the influence domain intersecting the solid boundary.

Rkl ðeÞ ¼ ½1gðxÞekl ðxÞ þ and, gðxÞ ¼

Of

gðx,x0 Þekl ðx0 ÞdVðx0 Þ

gðx,x0 ÞdVðx0 Þ,

0 r gðxÞ r1

Of

ð2Þ

Eq. (1) states that the stress at a point x in the solid depends on the value of the strain at all points in the medium. One can easily conceived that if the field point x0 is far enough from the evaluation point x, the stress s(x) is practically unaffected by the value of e(x0 ). This fact is taken into account through the definition of the attenuation function g(x,x0 ) which has the following properties (Lazar et al. [20]):

Z

Z

This stress–strain relation will also be investigated in this work. 2.2. Solution method In the absence of body forces, the equilibrium of the solid is described by the following local differential equation:

sij,j ¼ 0 - The nonlocal kernel has the dimension of (length)  3. Then it must depend on an internal characteristic length - It reaches its maximum at x ¼x0 and decay to zero at large distances R - It has to satisfy a normalization condition: O1 gðx,x0 ÞdVðx0 Þ ¼ 1, where ON is the infinite domain embedding the finite domain O. It is a continuous function of position. In the limit as the internal length goes to zero, it must become a Dirac delta function allowing to recovering the classical elasticity (local). In this work, the following attenuation function is adopted: 3

gðx,x0 Þ ¼ ð8p l Þ1 expð9x0 x9=lÞ where l is the internal length characteristic of the material under study. Having selected the value of l, a spherical region Of with radius R is defined such that: Z gðx,x0 ÞdVðx0 Þ  1 ð3Þ Of

Hereafter, R is taken to be 9 times the value of l. In this case, the normalization condition (3) is verified with an error less than 1% and for 9x0  x9 ZR, the attenuation function is practically zero (see Polizzotto et al. [24]). Note that this value of R is valid only for the attenuation function chosen in this study and must be recalculated if another function is used. Consider a solid of finite extent, there are points located near the boundary (e.g. points B and C in Fig. 1) for which the spherical influence region Of is truncated by the solid boundary. For such points, relation (3) is no more valid and the stress-strain relation

ð5Þ

Eq. (5) must be supplemented by proper boundary conditions given in terms of known displacement and traction. In the case of an isotropic solid considered in this work, relation (4) is put in the form:

sij ðxÞ ¼ C ijkl ekl ðxÞ þ C ijkl Q kl ðeÞ

ð6Þ

with Qkl(e) ¼Rkl(e)  ekl The field Eq. (5) are written as @ @ ðC e Þ þ ðC Q ðeÞÞ ¼ 0 @xj ijkl kl @xj ijkl kl

ð7Þ

It is tempting to try to solve Eq. (7) by the recently introduced numerical methods known as meshless or meshfree methods. Detailed of some of these methods can be found in the papers by Belytschko et al. [27], Atluri et al. [28], Liu et al. [29] and Liu et al. [30]. However, most of them use either a local weak form or a local integral equation of the field equations. Then either a volume or a surface domain of integration is required. Regarding nonlocal integral elasticity, the final strategy is similar to that of nonlocal FEM. An attractive meshless approach is the point collocation method applied to strong form field equations. The main shortcoming of this method as pointed out by Liu and Gu [31] is the deterioration of the solution accuracy in the presence of the Neumann type boundary condition (natural boundary conditions). In order to overcome this shortcoming, they proposed the weak–strong form approach which uses a local weak form of field equations at such boundary collocation points. In this work, we present an alternative to the weak–strong form point collocation method.

608

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

Now assume that the displacement field is the sum of a complementary field uC and a particular displacement uP. This partition is introduced into the first term of relation (7) and the following is considered:

For an internal collocation node, interpolation (13) is adopted for the particular integral and the displacement field. After some algebraic manipulations, Eq. (9) is written in the following matrix form:

@ ðC ec Þ ¼ 0 @xj ijkl kl

ð8Þ

½Q up fup=L g þ½Q u fu=L g ¼ f0g

@ @ ðC ep Þ þ ðC Q ðeÞÞ ¼ 0 @xj ijkl kl @xj ijkl kl

ð9Þ

2.2.1. Complementary solution Eq. (8) is similar to the standard Navier equations. Its boundary integral formulation is well established and reads: Z Z U mi ðx,yÞt ci ðxÞdGðxÞ ¼ T mi ðx,yÞðuci ðxÞuci ðyÞÞdGðxÞ ð10Þ G

G

Umi is the Kelvin fundamental solution and Tmi is the corresponding traction. Their expressions are provided in many textbooks (e.g. Brebbia [32], Balas et al. [33], Bonnet [34]). ! is the boundary of O. Eq. (10) is valid for boundary points as well as internal points. The application of usual approximation leads to a system of equations of the form: ½Hfuc g ¼ ½Gft c g

ð11Þ

ð14Þ

Collecting Eq. (14) for all internal collocation points, one obtains ( p) ( ) uB uB II IB II Q  Q  ½Q IB þ ½Q ¼ f0g ð15Þ p up up u u uI uI In relation (15), I stands for internal points while B indicates boundary nodes. Since we are looking for a particular solution, it can be selected such that fuPB g ¼ 0. It then follows that: ( ) uB p II IB II ¼ f0g ð16Þ ½Q up fuI g þ ½Q u Q u  uI

2.2.3. Boundary traction Let n denotes the outer unit normal at the boundary point x. The traction at point x is expressed, according to the proposed partition as: t i ðxÞ ¼ sij ðxÞnj ðxÞ ¼ C ijkl eckl ðxÞnj ðxÞ þC ijkl epkl ðxÞnj ðxÞ þ C ijkl Q kl ðeÞnj ðxÞ or more compactly

2.2.2. Particular displacement field The particular solution satisfies Eq. (9). Let the domain O and its boundary G be represented by properly scattered nodes. A kinematical field nk at any point x in the domain is interpolated by the local radial point interpolation as explained in the work by Liu and Gu [30]. That is:

nhk ðxÞ ¼

N X

Ri ðrÞaki þ

i¼1

M X

Pj ðxÞbkj

ð12  aÞ

P j ðxi Þaki ¼ 0

ðj ¼ 1 to MÞ

p ft p g ¼ ½K BI tp fuI g

ð17  aÞ (

uB uI

) ¼ ½K dt fug

ð17  bÞ

ð12  bÞ

i¼1

Ri(r) is the selected radial basis function, N the number of nodes in the neighborhood of point x and M is the number of monomials Pj(x) in the basis of the selected augmented polynomial degree. Coefficients aki and bkj can be determined by enforcing Eq. (12-a) to be satisfied at the N nodes in the neighborhood of point x. This leads to the solution of a system of equations of the form:    ( ) ak nk=L R P ¼ bk 0 PT 0 where {nk/L} is the vector of nodal values of nk(x): fnk=L g ¼ ðn11

Using the local radial point interpolation, the last two terms of the right hand side of the above relation, are written for all boundary nodes in the matrix form:

BI fdtg ¼ ½K BB dt K dt 

j¼1

with the constraints: N X

t i ðxÞ ¼ t ci ðxÞ þ t pi ðxÞ þ dt i ðxÞ

2.2.4. Final equation Consider the partition of the displacement field, Eq. (11) is rewritten as ½Hfug½Gftg ¼ ½Hfup g½Gft p þ dtg Using the expression of the particular integral as obtained from Eqs. (16) and (17), one obtains: ½Hfug½Gftg ¼ ð½H½UQ ½G½K tp ½UQ ½G½K dt Þfug II with ½UQ  ¼ ½Q IIup 1 ½Q IB u Qu The final form of the equation is similar to that of the conventional BEM, that is

~ ½Hfug½Gftg ¼ f0g

n12 n13 ::: nN1 nN2 nN3 Þ

ð18Þ

It now remains to take into account boundary conditions in the usual way.

It follows that, fbk g ¼ ð½PT ½R1 ½PÞ1 ½PT ½R1 fnk=L g ¼ ½F b fnk=L g and,

3. Numerical examples

fak g ¼ ½R1 ð½I½P½F b Þfnk=L g ¼ ½F a fnk=L g Approximation (12-a) is now written as

nhk ðxÞ ¼ ½R1 R2 UUU RN ½F a fnk=L g þ ½P1 P2 UUU PM ½F b fnk=L g or more compactly as

nhk ðxÞ ¼ ½FðxÞfnk=L g

ð13Þ

In this work, the generalized radial basis functions Ri ðrÞ ¼ ðr 2 þc2 Þq are adopted. r is the Euclidian distance between the field point x and the center xi. c and q are known as shape parameters. Let us point out that the approach uses two distinct domains of influence. The first one is related to the constitutive relation and

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

is selected smaller than the spherical region of influence for the local point interpolation. 3.1. Case of unit cube under uniform tension In the following, the stress–strain relation given by Eq. (1) is called Eringen’s model while the improved relation (4) is called Polizzotto’s model. In order to pinpoint the difference between the two models, we consider first the case of a unit cube under uniform tension load.

Fig. 2. Schematic of the loaded unit cube: (a) Uniform tensile loading; (b) uniform shear loading.

609

The homogeneous isotropic solid has a Young modulus E¼220 GPa and a Poisson ratio n ¼0.3. The characteristic internal length l is selected to be 10 mm. The latter value is chosen arbitrarily. Rigorously, it must be selected as a representative length of the inner microstructure of the material. A scheme of the considered specimen is given in Fig. 2(a). The unit cube is subjected to uniform displacement u in the z direction on the upper face (z¼0.5). The lower face (z ¼  0.5) is displacement constrained in the z direction and traction free in the tangential direction. The other faces of the cube are free of traction. This type of Dirichlet loading is selected since it allows a simple comparison between the two selected models. The analytical solution for szz is given by: szz ¼ EUu (since the length of one edge of the cube is 1). Note that for the presented results, the boundary of the cube is subdivided into 24 nine-node quadrilateral elements, i.e. four elements per face. The choice of the shape parameters is of great importance in the local radial point interpolation. Liu and Gu [30] propose to use the value q¼1.03 for static problems. The parameter c is problem dependent. The latter is written as c¼ ad0, where a is a dimensionless coefficient and d0 is the minimal distance between collocation nodes in the domain. As a first step in our analysis, the value q¼1.03 is adopted and a is set to 1.0.

Table 1 Uniform tensile loading of a unit cube: numerical results of strain and stress values at an internal point and a boundary point. Influence of the number of internal collocation points. Point

Fields

Analytical

Numerical (27 points)

Numerical (45 points)

Numerical (125 points)

Eringen’s model

Polizzotto’s model

Eringen’s model

Polizzotto’s model

Eringen’s model

Polizzotto’s model

(0, 0, 0)

exx (  10  3) eyy (  10  3) ezz (  10  3) szz (MPa)

 0.1363636  0.1363636 0.4545454 100

 0.1363634  0.1363634 0.4545446 99.05985

 0.1363635  0.1363635 0.4545450 99.99989

 0.1363630  0.1363630 0.4545448 99.06001

 0.1363635  0.1363635 0.4545450 99.99991

 0.1363639  0.1363639 0.4545447 99.05973

 0.1363635  0.1363635 0.4545450 99.99990

(0.5, 0.25, 0.25)

exx (  10  3) eyy (  10  3) ezz (  10  3) szz (MPa)

 0.1363636  0.1363636 0.4545454 100

 0.1363630  0.1363638 0.4545453 49.53003

 0.1363634  0.1363634 0.4545451 99.99994

 0.1363624  0.1363633 0.4545448 49.52999

 0.1363632  0.1363634 0.4545444 99.99977

 0.1363649  0.1363642 0.4545444 49.52972

 0.1363645  0.1363631 0.4545441 99.99952

Fig. 3. Evolution of szz (Eringen model): (a) On the top of the cube; (b) in the plane Oxy.

610

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

The displacement field in nonlocal elasticity is identical to that of conventional elasticity. Then for the considered loading, the analytical solution is available and is used for comparison purpose. The calculation has been done with 27, 45 and 125 regularly distributed internal collocation nodes. In all cases, the displacement values at all points are practically analytical solution. The non-zero components of the strain and stress tensors are shown in Table 1 for an internal node and a boundary point. For this loading state, the calculated strains are practically analytical solution. The effect of solid boundary (truncation of the nonlocality domain of influence), corrected by the improved model is clearly demonstrated. In Fig. 3(b), szz calculated by Eringen’s model, is represented for points in the middle plane perpendi-

cular to the loading direction. As can be observed, the calculated values are expected ones except at those points near the solid boundary. This effect is more clear when we consider the upper surface of the cube (Fig. 3(a)). Near and at the solid boundaries, the calculated stresses are lower than the expected value. When the nonlocality influence domain is not truncated by the solid boundary, the model leads to results that are highly accurate. There is no solid boundary effect when using the improved model of Polizzotto. szz profiles along a line on the top of the cube and another through the solid are provided as a further illustration (Fig. 4(a) and (b)). Generalized multi-quadrics radial basis functions depend on two shape parameters q and c, which can affect the accuracy

Fig. 4. (a) Evolution of szz between (  0.5;0;0.5) and (0.5;0;0.5); (b) evolution of szz between ( 0.5;0;0) and (0.5;0;0).

Table 2 Influence of shape parameter c on the strain field at point (0.5;0.25;0.25) for uniform traction loading (Polizzotto’s model). c (q ¼1.03)

Fields (  10  3)

Numerical (27 internal nodes)

1  d0

exx eyy ezz

 0.1363630  0.1363638 0.4545453

0.01  d0

exx eyy ezz

 0.1363634  0.1363636 0.4545451

0.0001  d0

exx eyy ezz

 0.1363634  0.1363636 0.4545451

Table 3 Uniform shear loading of a unit cube: numerical results of strain and stress values at an internal point and a boundary point. Influence of the number of internal collocation points. Point

Fields

Analytical

Numerical (27 points)

Numerical (125 points)

Eringen’s model

Polizzotto’s model

Eringen’s model

Polizzotto’s model

(0, 0, 0)

eyz (  10  3) syz (MPa)

0.0590909 100

0.0590927 99.06135

0.0590908 99.99981

0.0590905 99.05974

0.0590906 99.99981

(0.5, 0.25, 0.25)

eyz (  10  3) syz (MPa)

0.0590909 100

0.0590913 49.52990

0.0590908 99.99983

0.0590906 49.52988

0.0590910 99.99992

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

of the solution. It is well known that the parameter c affects more the solution that q. Thus, it is desirable that reasonable accuracy be maintained within a large range of variation of c. Keeping the value of q at 1.03, c has been varied from 0.0001  d0 to 1  d0. In this case of uniform loading, the results are practically insensitive to the value of this parameter. As an example, in Table 2, the strain values calculated at a boundary point are presented for different values of c and the worst case of number of internal collocation points. Numerical simulations with different values of the shape parameter q have led to the same conclusions.

Table 4 Influence of shape parameter c on the strain field at point (0.5;0.25;0.25) for uniform shear loading (Polizzotto’s model). c(q¼ 1.03)

Fields (  10  3)

Numerical (27 internal nodes)

10  d0 1  d0 0.01  d0 0.0001  d0

eyz eyz eyz eyz

0.05909436 0.05909084 0.05909073 0.05909072

611

3.2. Case of a unit cube under simple shear stress We consider now the case of the same unit cube under shear loading. Always for comparison purpose, the loading is of the Dirichlet type. As shown in Fig. 2(b), the cube is subjected to uniform displacement in the y direction on the upper and lower faces (z¼ 70.5). The same constant displacement is applied in the z direction on lateral faces (y¼ 70.5) with free traction in the tangential direction. The lateral faces (x ¼ 70.5) are constrained in the x direction. The numerical results obtained with Eringen’s and Polizzotto’s models are shown in Table 3 for an internal collocation point and a boundary node. Once again, Eringen’s model gives low stress values at boundary nodes whereas the strain and stress obtained by Polizzotto’s model are practically analytical solutions whatever the position of the point and the number of internal collocation centers. In this case also, high accuracy of results is maintained for a wide range of values of c, from 0.0001  d0 to 10  d0. As an example, the value of eyz at a boundary point is provided in Table 4 for different values of the shape parameter c in the case of 27 internal collocations points. Compared with analytical solution, the maximum error (for c ¼10  d0) is less than 0.006%. 3.3. Case of a cylindrical bar under torsion

Fig. 5. Schematic of the cylinder under torsion loading.

A cylindrical bar with radius 1 mm and height 10 mm is considered (see Fig. 5). In order to simulate torsion, the tangential displacement at the upper and lower face of the specimen is specified. More specifically, let H be the half height of the cylinder and at be the twist angle per unit length. At the upper face ux ¼  atHy and uy ¼ atHx while at the lower ux ¼ atHy and uy ¼  atHx. The lateral surface of the cylinder is free of traction. In this case, the analytical solution for syz is given by syz ¼ mrat, where m is the shear modulus, linked to E and n by: m ¼E/2(1þ n). Along the x-axis, syz can be expressed in terms of x: syz ¼ m9x9at. The simulation uses 80 nine-node quadrilateral boundary elements and 100 internal nodes.

Fig. 6. (a) Evolution of syz between ( 1;0;10) and (1;0;10); (b) evolution of syz between ( 1;0;7.5) and (1;0;7.5).

Table 5 Torsion loading of a cylinder: numerical results of stress values at an internal point and a boundary point. Influence of the number of internal collocation points. Point

Stress (MPa)

Analytical

Numerical (100 internal points)Polizzotto’s model l¼ 1 mm

l¼ 10 mm

l¼ 20 mm

(1, 0, 8)

sxz syz

0 0.338462

o10  5 0.3382272

o 10  4 0.3360592

o10  4 0.3336628

(0.3535, 0.3535, 8)

sxz syz

 0.119664 0.119664

 0.1196661 0.1196689

 0.1197072 0.1196875

 0.1197038 0.1197099

612

M. Schwartz et al. / Engineering Analysis with Boundary Elements 36 (2012) 606–612

Several values of the characteristic internal length have been tested l ¼1, 10 and 20 mm using both Eringen’s and Polizzotto’s models. In Fig. 6(a), the variation of the stress component syz along a diameter of the top surface is presented. Once more, the effect of truncation of the nonlocality domain of influence is evidenced and the higher the internal length the smaller the stress value. In Fig. 6(b), the variation of syz along a diameter of a cross section within the cylinder is presented. The conclusions are similar to those of the loading of the unit cube. As shown in Table 5, Polizzotto’s model leads to results very close to the classical solution with a difference less than 2%. Note that, for this loading also, the quality of the results is not affected by the parameter c.

4. Conclusion Consider classical homogeneous isotropic problem governed by linear differential equation, the boundary element method generally leads to highly accurate results. In the case of non- local integral elasticity a fundamental solution is not yet available and the method seems inapplicable. The present work proposes a promising BEM based solution procedure. It couples the conventional approach for isotropic elastostatic with the local radial point interpolation of a strong form differential equation. The overall approach is simple to implement. Its effectiveness and accuracy has been demonstrated on the example of loading of a unit cube and torsion of a cylindrical bar. The application of the method to fracture mechanics problems and nano-materials will be the subject of a future work. References [1] Sudak LJ. Column buckling of multiwalled carbon nanotubes using nonlocal continuum mechanics. J Appl Phys 2003;94:7281–7. [2] Wang Q, Varadan VK. Application of nonlocal elastic shell theory in wave propagation analysis of carbon nanotubes. Smart Mater Struct 2008;16: 178–90. [3] Filiz S, Aydogdu M. Axial vibration of carbon nanotube heterojunctions using nonlocal elasticity. Comp Mater Sci 2010;49:619–27. [4] Hu Yan-Gao, Liew KM, Wang Q, He XQ, Yakobson BI. Nonlocal shell model for elastic wave propagation in single and double walled carbon nanotubes. J Mech Phys Solids 2008;56:3475–85. [5] Eringen AC. Theory of nonlocal elasticity and some applications. Res Mech 1987;21:313–42. ¨ [6] Kroner E. Elasticity theory of materials with long range cohesive forces. Int J Solids Struct 1967;3:731–42.

[7] Kunin IA. The theory of elastic media with microstructure and the theory of ¨ dislocation. In: Kroner E, editor. Mechanics of Generalized Continua; 1968. [8] Krumhansl JA. Some consideration on the relations between solid state ¨ physics and generalized continuum mechanics. In: Kroner E, editor. Mechanics of Generalized Continua; 1968. [9] Edelen DGB, Laws N. On the thermodynamics of systems with nonlocality. Arch Rationale Mech Anal 1971;43:24–35. [10] Edelen DGB, Green AE, Laws N. Nonlocal continuum mechanics. Arch Rationale Mech Anal 1971;43:36–44. [11] Eringen AC, Edelen DGB. On nonlocal elasticity. Int J Eng Sci 1972;10:233–48. [12] Eringen AC, Kim BS. Stress concentration at the tip of a crack. Mech Res Commun 1974;1:233–7. [13] Eringen AC, Speziale CG, Kim BS. Crack tip problem in nonlocal elasticity. J Mech Phys Solids 1977;25:339–55. [14] Eringen AC. Line Crack subjected to shear. Int J Fract 1978;14:367–79. [15] Eringen AC. Line Crack subjected to anti-plane shear. Eng Fract Mech 1979;12:211–9. [16] Gao J, Dai TM. On the study of a penny-shaped crack in nonlocal elasticity. Acta Mech Solids Sin 1989;2:253–69. [17] Minghao Z, Changjun C, Yuanjie L, Guoning L, Shishan Z. The method of analysis of crack problem in three dimensional non-local elasticity. Appl Math Mech 1999;20. [18] Eringen AC. Edge dislocation in nonlocal elasticity. Int J Eng Sci 1977;15: 177–83. [19] Eringen AC. Screw dislocation in nonlocal elasticity. J Phys D Appl Phys 1977;10:671–8. [20] Lazar M, Maugin GA, Aifantis EC. On a theory of nonlocal elasticity of biHelmholtz type and some applications. Int J Solids Struct 2006;43:1404–21. [21] Povstenko YZ, Kubik I. Concentrated ring loading in a nonlocal elastic medium. Int J Eng Sci 2005;43:457–71. [22] Artan R. The nonlocal solution of the elastic half plane loaded by a couple. Int J Eng Sci 1999;37:1389–405. [23] Polizzotto C. Nonlocal elasticity and related variational principles. Int J Solids Struct 2001;38:7359–80. [24] Polizzotto C, Fushi P, Pisano AA. A strain-difference-based nonlocal elasticity model. Int J Solids Struct 2004;41:2383–401. [25] Jackiewicz J, Holka H. Computational simulation of crack problems by means of the contour element method. Eng Fract Mech 2008;75:461–74. [26] Pisano AA, Sofi A, Fuschi P. Nonlocal integral elasticity: 2D finite element based solutions. Int J Solids Struct 2009;46:3836–49. [27] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [28] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [29] Liu WK, Jun S, Li S, Andee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numer Methods Eng 1995;38: 1655–79. [30] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246:29–46. [31] Liu GR, Gu YT. A meshfree method: meshfree weak–strong (MWS) form method for 2-D solids. Comput Mech 2003;33:2–14. [32] Brebbia CA, Dominguez J. Boundary Elements. An Introductory Course. Computational Mechanics Publications; 1992. [33] Balas J, Sladek J, Sladek V. Stress Analysis by Boundary Element Methods. Elsevier; 1989. [34] Bonnet M. Boundary Integral Equation Methods for Solids and Fluids. NewYork: John Wiley and Sons; 1999.