Boundary-only element solutions of 2D and 3D nonlinear and

Engineering Analysis with Boundary Elements 31(2007)974–982

Boundary-only element solutions of 2D and 3D nonlinear

and nonhomogeneous elastic problems

X.W.Gao a,?,Ch.Zhang b ,L.Guo a

Department of Engineering Mechanics,Southeast University,Nanjing,210096China b

Department of Civil Engineering,University of Siegen,D-57068Siegen,Germany

Received 2November 2006;accepted 10May 2007

Available online 18September 2007

Abstract

This paper presents a robust boundary element method (BEM)that can be used to solve elastic problems with nonlinearly varying material parameters,such as the functionally graded material (FGM)and damage mechanics problems.The main feature of this method is that no internal cells are required to evaluate domain integrals appearing in the conventional integral equations derived for these problems,and very few internal points are needed to improve the computational accuracy.In addition,one of the basic ?eld quantities used in the boundary integral equations is normalized by the material parameter.As a result,no gradients of the ?eld quantities are involved in the integral equations.Another advantage of using the normalized quantities is that no material parameters are included in the boundary integrals,so that a uni?ed equation form can be established for multi-region problems which have different material parameters.This is very ef?cient for solving composite structural problems.r 2007Elsevier Ltd.All rights reserved.

Keywords:Boundary element method;Nonhomogeneous problem;Functionally graded materials;Isotropic damage mechanics;Radial integration method

1.Introduction

Boundary element method (BEM),an attractive and promising numerical analysis tool due to its semi-analytical nature and the reduction of the problem dimension,has been successfully applied to homogeneous,isotropic and linear problems.However,its extension to nonhomoge-neous and nonlinear problems is not straight forward,since no universal fundamental solutions are available for these problems.Consequently,one has to apply the fundamental solutions for homogeneous and linear problems to estab-lish integral equations for the nonhomogeneous and nonlinear problems.As a result of this treatment,domain integrals associated with the nonhomogeneity and non-linearity are inevitably introduced in the integral equations.Direct evaluation of the domain integrals requires the discretization of the domain into internal cells [1,2].This severely eliminates the distinct advantage (namely the boundary-only discretization)of BEM.

To circumvent the domain discretization,the commonly used method is to transform the domain integrals into equivalent boundary integrals.In this method,the dual reciprocity method (DRM)by Nardini and Brebbia [3]is extensively used.However,DRM requires particular solutions,which restricts its application to complicated problems [4].Recently,a new transformation method,the radial integration method (RIM),was developed by Gao [5,6],which can transform any domain integral into a boundary integral and a radial integral.The radial integral is independent of geometry,and therefore only the global boundary of the problem needs to be discretized into boundary elements to evaluate the boundary integral.Afterwards,RIM was used to solve thermoelasticity problems [7],elastic inclusion problems [8],elastoplastic frictional contact problems [9]and creep damage me-chanics problems [10].In view of its robustness and simplicity in evaluating domain integrals without using internal cells,Hematiyan [11]gives RIM a very good assessment and Albuquerque et al.[4]compared RIM to DRM numerically through dynamic problems of

https://www.360docs.net/doc/d82453014.html,/locate/enganabound

?Corresponding author.Tel.:+862583792412;fax:+862583615736.

anisotropic plates.More recently,RIM is applied to solve varying material property problems for heat conduc-tion problems[12]and two-dimensional(2D)elasticity problems[13]where the heat conductivity and shear modulus are regarded as spatial functions.In these works, normalized temperature(product of temperature and conductivity)and normalized displacements(product of displacements and shear modulus)are introduced to formulate boundary integral equations.The advantage of using the normalized quantities is that no gradients of the ?eld quantities are involved in the integral equations,nor are there material parameters in the boundary integral terms.Therefore,not only the derived formulation can be easily coded,but also the accuracy of the computational results can be improved considerably.

In this paper,RIM is used to solve2D and three-dimensional(3D)nonhomogeneous elasticity problems in which the shear modulus may be a function of stresses. Since this is a nonlinear problem,the iterative solution procedure using Newton–Raphson method is employed to solve the system of equations.In iterative computation, stresses are required.Therefore,a new integral equation with a free term for evaluating stresses is derived in the paper.The arrangement of this paper is as follows.Firstly, the RIM for transforming a domain integral to the boundary is reviewed.Secondly,boundary-domain integral equations for nonlinear and nonhomogeneous problems are derived using the concept of normalized quantities. Then,the arising domain integrals are converted into boundary integrals by using RIM,resulting in a boundary-only element algorithm.Finally,two numerical examples are given to demonstrate the accuracy and the ef?ciency of the presented method.

2.Review of radial integration method(RIM)to evaluate domain integrals without internal discretization

In terms of RIM,a domain integral with the integrand f(x,x p)can be transformed into an equivalent boundary integral[5,6]as follows:

Z O fex;x pTd OexT?

Fex q;x pT

r aex q;x pT

q r

q n

d Gex qT,(1)

where a?1for2D and a?2for3D problems,r(x q,x p) denotes the distance between the source point x p and the boundary point x q,and F(x q,x p)is determined by the following radial integral:

Fex q;x pT?

Z rex q;x pT

fex;x pTr a d r.(2)

In Eqs.(1)and(2),O is the domain of the problem under consideration,G is the boundary of O,n denotes the outward normal vector to the boundary G,x is the?eld point in O and the integral variable r in Eq.(2)is the distance from the source point x p to the?eld point x.Fig.1shows the relationship between the source point x p,?eld point x,boundary point x q and the distance r.

The quantities included in Eqs.(1)and(2)can be determined using

q r

q n

?^rán?r;i n i,(3) ^r?

q r

q x

x qàx p

rex q;x pT

r;i?

q r

q x i

x q

àx p

rex q;x pT

,e4T

^rá^r?r;i r;i?1.(5) In order to evaluate the radial integral(2),the co-ordinates x in the function f(x,x p)need to be expressed in terms of the integral variable r using

x?x pt^r r or x i?x p itr;i r.(6) An important note should be made that the quantity r,i appearing in the above equations is a constant vector for the radial integral(2),although it is a function of x q for the boundary integral(1).

Usually,Eq.(2)can be analytically integrated by using Eq.(6).For instance,consider the3D example of f?r,i x j/r.From Eqs.(6)and(2)it follows that F?

R r

r;iex p

tr;j rTr d r?r;i r2ex p

=2tr;j r=3T,and then in terms of Eq.(1),we obtain

r;i x j

d O?

r;i

x p

r;j r

q r

q n

d G.(7)

For some complicated kernel functions,it may be not easy to integrate Eq.(2)analytically.In these cases, numerical integration techniques,such as Gaussian quad-rature,may be used to evaluate this integral for each boundary integral point x q.To use the Gaussian quad-rature,the following variable transformation is required:

rex q;x pT

e1txTeà1p x p1T.(8)

Fig.1.Relationship between points and distances.

X.W.Gao et al./Engineering Analysis with Boundary Elements31(2007)974–982975

Using this transformation and Eq.(6),integral(2)can be expressed as

Fex q;x pT?

rex q;x pT

at1Zt1

à1

fexexT;x pTe1txTa d x

rex q;x pT

at1X N r

l?1

fexex lT;x pTe1tx lTa w l,e9T

where N r is the number of Gaussian points,x l is the Gaussian coordinate,w l is the associated weight and x(x l)is determined by substituting Eq.(8)into Eq.(6).

From Eq.(7),it can be seen that the RIM not only can transform a domain integral to the boundary,but also can remove integral singularities[6].

3.Boundary-domain integral equations for nonhomogeneous elasticity problems

Let us consider an isotropic linear elastic solid with the shear modulus m being dependent on Cartesian coordinates x and Poisson’s ratio n being constant.Under these assumptions,we can write the elasticity tensor as

C ijklexT?mexTC0

ijkl

,(10) where

ijkl ?

1à2n

d ij d kltd ik d jltd il d jk.(11)

In which,d ij denotes the Kronecker delta.The relation-ship between the stress tensor s ij and the displacement gradients u k,l can be expressed as

s ij?m C0

ijkl

u k;l,(12) where u k,l?q u k/q x l.Throughout the paper,a comma after a quantity represents spatial derivatives and repeated subscripts denote summation.The traction vector t i on the boundary of the considered domain is related to the stress components by

t i?s ij n j.(13) Based on the equilibrium equations s ij,j?0without considering body forces,the following domain integral with the weight function U ij(x,x p)can be written as

U ijex;x pTs jk;kexTd O?0.(14)

Substituting Eq.(12)into Eq.(14)and using Gauss’s divergence theorem,it follows that[13]

~u iex pT?

G U ijex;x pTt jexTd Gà

T ijex;x pT~u jexTd G

V ijex;x pT~u jexTd O,e15T

where U ij and T ij are the Kelvin displacement and traction fundamental solutions[1,2]for an isotropic, homogeneous and linear elastic medium with m?1,and

V ij?

à1

4pae1ànTr a

f~m;k r;k?e1à2nTd ij

tb r;i r;j te1à2nTe~m;i r;jà~m;j r;iTge16Tin which b?2for2D and b?3for3D problems and a?bà1,~u jexTand~m are normalized displacements and shear modulus,i.e.,

~u iexT?mexTu iexT;~mexT?log mexT.(17) Note here that Eq.(15)is an integral representa-tion for arbitrary internal points.As described in many text books on BEM(e.g.,[1,2]),boundary integral equations for boundary points can be obtained by letting x p-G.

4.Integral equation for internal stresses

According to Eq.(17),it can be shown that

q u k

q x p

1q~u k

q x p

à~u k

q~m

q x p

.(18)

Using this equation and according to Eq.(12),it follows that

s ij?C0

ijkl

q~u k

q x p

àC0

ijkl

~m;l~u k.(19)

To evaluate stresses,the term q~u k=q x p

needs to

be evaluated.Differentiating Eq.(15)with respect to x p

l yields

q~u kex pT

q x p

q U ksex;x pT

q x p

t sexTd G

q T ksex;x pT

q x p

~u sexTd G

q x p

V ksex;x pT~u sexTd O.e20T

All of these differentials,except the last,can be evaluated without dif?culty.However,the last domain integral is strongly singular and a‘jump’term arises during this process.To deal with this integral,a small sphere O e with radius e,and centred at x p,is separated from O,as

Ωε

Γε

Fig.2.A small sphere O e cut out from O.

X.W.Gao et al./Engineering Analysis with Boundary Elements31(2007)974–982 976

depicted in Fig.2.The last term in Eq.(20)is separated into two parts,thus

q q x p l

Z O V ks ex ;x p

T~u s ex Td O ?lim !0Z O àO q V ks ex ;x p Tq x p

~u s ex Td O t~u s ex p Tlim !0Z O q V ks ex ;x p T

q x p l

d O .

e21T

Noting that q eáT=q x p l ?àq eáT=q x l and using Gauss’

theorem,it follows that Z O q V ks ex ;x p T

q x p

l d O ?àZ

G V ks ex ;x p Tn l ex Td G .(22)Noting that on the boundary G e ,n l ?r ,l and r ?e ,the

last term can be analytically integrated to get Z

O q V ks ex ;x p Tq x p

d O ?1

8e1àn Tfe3à4n Ted ks ~m ;l td sl ~m

;k Tàe1à4n Td kl ~m ;s g e23a T

for 2D problems,and

O q V ks ex ;x p Tq x p

d O ?1

15e1àn Tfe4à5n Ted ks ~m ;l td sl ~m

;k Tàe1à5n Td kl ~m ;s g e23b T

for 3D problems.

Combining these results with the differentiations indi-cated in Eq.(20),we can express Eq.(19)in the form

s ij ex p T?Z G U ijk ex ;x p

Tt k ex Td G àZ G T ijk ex ;x p T~u

k ex Td G tZ

V ijk ex ;x p T~u k ex Td O tF ijk ex p T~u k ex p T,e24T

where the kernels U ijk and T ijk are the same as in the

homogeneous problems [1,2]with m ?1,and V ijk ?

2pa e1àn Tr f b ~m

;m r ;m ?e1à2n Td ij r ;k tn ed ik r ;j td jk r ;i Tàg r ;i r ;j r ;k

tbn e~m

;i r ;j t~m ;j r ;i Tr ;k àe1à4n T~m ;k d ij te1à2n Teb ~m

;k r ;i r ;j t~m ;j d ik t~m ;i d jk Tg ,e25T

F ijk ?

à14e1àn Tf d ij ~m

;k td ik ~m ;j td jk ~m ;i g ;for 2D ;à1

15e1àn Tfe2t10n Td ij ~m

;k te7à5n Ted ik ~m ;j td jk ~m ;i Tg ;for 3D :

8<:(26)

The domain integrals in Eq.(24)are interpreted in the Cauchy principal value sense.This means that cutting out an in?nitesimal region around the source point x p does not change the integration result.Based on this property,the conventional singularity-separation scheme [1]can be used to regularize the strongly singular domain integral involved in Eq.(24).

The stress integral representation formula (24)is only applicable to evaluate stresses at internal points.For

boundary points when the source point x p approaches the

?eld point x ,high order singularities arise.Although these high-order singular integrals can be directly evaluated using the methods described in [14,15],the existing Fortran subroutines described in [1]using the ‘traction recovery’method is adopted in this paper to computed the boundary stresses.

5.Transformation of domain integrals to the boundary using RIM

In order to avoid the discretization of the domain into internal cells for evaluating the domain integrals appearing in Eqs.(15)and (24),the RIM described in the preceding section is used to transform them into boundary integrals.To do this,the normalized displace-ments in Eqs.(15)and (24)are approximated by a series of prescribed basis functions as commonly used in DRM [16].Thus,

~u i ex T?X

a A i f A eR =S A Tta k i x k ta 0

i ,

(27)X

a A i ?0,

(28a)

a A i x A

j ?0,

(28b)

where R ?x àx A

is the distance from the application point x A to the ?led point x ,S A is the size of the support region for the radial basis function (RBF)at point x A ,and

a A i and a k

i are coef?cients to be determined.Investigations show that the 4th-order spline RBF [17]can give very stable results,which is given by

f A eR =R A T?

1à6R S A 2t8R S A 3à3R S A

40p R p S A ;

0R X S A

8<:(29)

In this study,the support size S A is determined by

specifying the number of application points around x A ,while

the coef?cients a A i and a k

i are determined by collocating x in Eq.(27)at all application points.As a result,together with Eq.(28),a set of algebraic equations can be resulted in as follows in the matrix form f ~u

g ??f f a g ,(30)

where f ~u

g is a vector of the normalized displacements at all application points,and {a }is a vector consisting of the

coef?cients a A i for all application points and a k

i .If there are no any two application points having the same coordinates,the matrix [f ]is invertible and thereby f a g ??f à1f ~u

g .(31)

Substitution of Eq.(27)into the displacement domain

integral of Eq.(15)and applying RIM Eqs.(1)and (2),it

X.W.Gao et al./Engineering Analysis with Boundary Elements 31(2007)974–982

977

follows that

Z O V ij~u j d O?

a A

r a

q r

q n

F A

d Gta k

r;k

r a

q r

q n

d G

tea k

x p

ta0

q r

q n

d G,e32T

F A ij ?

Z r

r a V ij f A d r,(33a)

F1 ij ?

Z r

r b V ij d r,(33b)

F0 ij ?

Z r

r a V ij d r.(33c)

Similarly,the domain integrals included in the stress integral equations(see Eq.(24))can be written as

Z O V ijk~u k d O?

a A k

q r

q n

F A

ijk

d Gta l k

r;l

q r

q n

ijk

d G

tea l

x p

ta0

r a

q r

q n

ijk

d G,e34T

F A ijk ?

Z r

r a V ijk f A d r,(35a)

F1 ijk ?

Z r

r b V ijk d r,(35b)

F0 ijk ?

Z r

r a V ijk d r.(35c)

To evaluate the radial integrals shown in Eqs.(33)and (35),the following relationships are used:

???????????????????????????

r2t2srtˉR2

;s?r;iˉR i;ˉR i?x p iàx A

.(36)

From Eqs.(16)and(33)it can be seen that all radial integrals in Eq.(33)are regular and can be integrated numerically by using the standard Gaussian quadrature formulation(9)without any dif?culty.However,from Eqs.(25)and(35)we can see that the radial integrals for

F A ijk and F0ijk in Eqs.(35)are strongly singular with the

order(1/r)when the source point approaches the?eld point.Therefore,a special technique,referred as singularity separation technique,is required to evaluate these radial integrals.The operation procedure is exactly the same as in the treatment of initial stress domain integrals in plasticity (see Ref.[5]for details).6.Application to functionally graded materials(FGMs) FGMs represent a new class of composite materials designed to achieve high performance levels superior to that of homogeneous materials by combining the desirable properties of each constituent[18,19].Such materials possess continuously graded properties with gradual change in a microstructure and have wide-range innovative applications in engineering sciences.The macro-material properties of FGMs, e.g.the shear modulus m,are functions of spatial https://www.360docs.net/doc/d82453014.html,ually,exponential spatial variations are used[13],i.e.,

mexT?m0e a xtb ytg z(37) or

mexT?m0e cr,(38) where a,b,g and c are constants,x,y,z,are the spatial coordinates and r is a distance measured from a speci?ed position.Apparently,for the relationship shown in Eq.(37),the term~m;i is constant and all radial integrals in Eqs.(33)and(35)can be analytically integrated, resulting in a considerable saving of the computational time.The relationship(37)have been used in[13]for a2D plate analysis and relationship(38)will be used in this paper for a3D hollow cylinder under internal pressure in the section of numerical examples.

7.Application to isotropic damage mechanics

In damage mechanics,the material properties change as the damage evolves.Let us consider a homogeneous, isotropic and linear elastic solid with the shear modulus m being dependent on the damage variable D which is regarded as the ratio of the damage volume to the total volume of the material.The following rule is adopted:

meDT?m0e1àDT,(39) where m0is the shear modulus corresponding to the undamaged material state.

In view of the dependence of the stress components on the orientation of the coordinate system,the damage variable D is assumed to be a function of stress invariants [1].The equivalent stress s to the Von Mises yield criterion is taken as the variable of D.The following damage evolution equation is used:

D?1àeàes=s0Tm,(40) where m and s0are material constants,and

?????????????

s0ij s0ij

;s0

?s ijà

s kk d ij.(41) Thus,the shear modulus can be expressed as

mesT?m0eàes=s0Tm.(42) The derivatives of the normalized shear modulus~m with respect to global coordinates in Eqs.(16)and(25)can be

X.W.Gao et al./Engineering Analysis with Boundary Elements31(2007)974–982 978

computed using

~m;i?àm

mà1q

q x i

.(43)

To compute spatial derivatives involved in the right-hand side term of Eq.(43),RBF is used to approximate the equivalent stress,i.e.,

s s0?

b A f AeR=S ATtb k x ktb0,(44)

b A?0,(45a)

X A b A x A

?0,(45b)

where the coef?cients b A and b k can be determined by collocating the?eld point through all application points similar to the determination of the coef?cients in Eq.(27). Thus

q q x i

b A

q f A

q R

q x i

tb i.(46)

Substituting Eq.(29)into the above equation yields

q q x i

à12

b A1àR

S A

ex iàx A

tb i0p R p S A;

b i R X S A:

(47)

Fig.3.A hollow cylinder under internal

pressure.

Fig.4.BEM mesh for the hollow cylinder.

0.4

0.8

1.2

1.6

b/ a

Fig. 5.Variation of radial displacement on inner surface with non-

homogeneity ratio.

0.4

0.8

1.2

1.6

012345

b/ a

Fig. 6.Variation of radial displacement on outer surface with non-

homogeneity

ratio.

(

)

0.4

0.8

1.2

1.6

b/ a

Fig.7.Variation of hoop stress on inner surface with nonhomogeneity

ratio.

X.W.Gao et al./Engineering Analysis with Boundary Elements31(2007)974–982979

8.System of equations and solution scheme

By discretizing the boundary G into a series of boundary elements,introducing some internal nodes and collocating the source point at every boundary and internal node,a system of equations expressed in terms of boundary

unknowns and internal normalized displacements can be formed [13].For FGMs problems,usual linear equation set solvers can be directly applied to solve the system of equations.However,for the damage mechanics problems,since the shear modulus depends on the stresses,the resulting system of equations is a nonlinear equation set.An iterative solution scheme is required to solve the equation set.In this paper,the Newton–Raphson iterative scheme (e.g.[1])is adopted and the norm of the equation set is used to judge if the iteration convergence is achieved.9.Numerical example for isotropic damage mechanics Two numerical examples are considered:the ?rst one deals with a 3D hollow cylinder under internal pressure consisting of FGMs;and the second one is concerned with a 2D plate under shear deformation for isotropic damage analysis.In both examples,four points per direction are used to determine the support size of RBF.To validate the correctness,computational results are compared with existing ones and/or with the ?nite element method (FEM)results.

9.1.Functionally graded material problem

In the ?rst example,a nonhomogeneous hollow cylinder under a pressure of p ?1on the internal surface as shown in Fig.3is analyzed.The cylinder has the dimensions of a ?4,b ?5,h ?0.8.The ?nite-length hollow cylinder is considered as a part of the in?nite-length tube.In this case,the plane strain condition holds true.An exponential spatial variation of shear modulus is considered as m ex T?m 0e c er àa T;where c ?e1=eb àa TTln em b =m a T.The Poisson’s ratio is assumed to be constant n ?0.25and m 0?4000which corresponds to Young’s modulus E 0?10000.

Exploiting symmetry,only a quarter of the cylinder needs to be discretized into boundary elements.The discretization scheme consists of 136quadratic boundary

H o o p s t r e s s (r =b )

0.40.81.21.62

Fig.8.Variation of hoop stress on outer surface with nonhomogeneity ratio.

X 1

Fig.9.A rectangular plate subjected to a tensile loading.

Fig.10.BEM and FEM meshes:(a)BEM mesh with 406internal nodes;(b)BEM mesh with 28internal nodes;(c)BEM mesh without internal nodes;and (d)FEM mesh with 861nodes.

X.W.Gao et al./Engineering Analysis with Boundary Elements 31(2007)974–982

980

elements,de?ned by 410nodes,as depicted in Fig.4.This is a 3D problem.To simulate the plane strain condition,the top and bottom surfaces are applied with sliding

condition along radial direction and no deformation occurs along z -direction.The outer surface of the cylinder with the normal along the radial direction is traction free.

The computational results are shown in Figs.5–8,in which Figs.5and 6are plots of the radial displacements on inner and outer surfaces of the cylinder,respectively,versus the radio of m b to m a and Figs.7and 8are the plots of the hoop (circumferential)stresses.The displacements and stresses in these ?gures have been divided by corresponding analytical solutions for homogenous materials [20].For comparison purpose,the results using MLPG method [21]and using FEM software ABAQUS are also plotted in ?https://www.360docs.net/doc/d82453014.html,parison shows that good agreements between the three methods are obtained,which veri?es the correctness and accuracy of the present method.9.2.2D isotropic elastic damage problem

The second example is concerned with the isotropic elastic damage problems.A rectangular plate with the dimensions of L ?W is analyzed (Fig.9).The investigated plate is subjected to a uniform tensile loading s 1as depicted in Fig.9.The plate is discretized into 90linear boundary elements with 90boundary nodes.The shear modulus changes according to Eq.(42).The Poisson’s ratio is selected as n ?0.25,undamaged shear modulus is m 0?0.4which corresponds to Young’s modulus E 0?1,and the damage parameters are s 0?15and m ?1.

To examine the dependence of the computational results on the number of internal nodes,three BEM models with 406,28and 0internal nodes,respectively,are computed (see Fig.10(a)–(c)).For the purpose of veri?cation,this problem was also computed using the FEM software ABAQUS.Fig.10(d)shows the used FEM mesh.

Fig.11shows the computed displacement distribu tion over the top side using different internal nodes,while Fig.12is the damage variable D along the bottom side.Fig.13gives the contour plot of the damage variable D

0.1

0.20.30.40.50.60.7

x 1 coordinate

D i s p l a c e m e n t u 1

Fig.11.Displacement distribution along x 1-direction (x 2?0.4).

00.050.10.150.2

x 1 coordinate

D a m a g e f a c t o r D

Fig.12.Damage distribution along x 1-direction (x 2?0)

Fig.13.Contour plot of damage variable D .

X.W.Gao et al./Engineering Analysis with Boundary Elements 31(2007)974–982

981

over the entire region.The iteration convergence history is shown in Fig.14.From Figs.11and 12it can be seen that the current boundary-only method is quite insensitive to the number of internal nodes.10.Conclusions and discussions

In this paper,a boundary element analysis approach for solving 2D and 3D nonlinear and nonhomogeneous problems is presented.A new integral equation is derived for evaluating internal stresses.Normalized displacements are introduced in the boundary-domain integral formula-tion,which avoids displacement gradients in the domain integrals.To transform domain integrals into boundary integrals along the global boundary of the analyzed domain,the RIM is applied,which results in a boundary-only element scheme.The normalized displacements in the domain integrals are approximated by the 4th-order spline-type RBF.Detailed application to FGMs and isotropic elastic damage problems is described.Two numerical examples are given to show that the present boundary-only element method is easy to implement and requires only very few internal nodes.

Computational results show that the current boundary-only method is quite insensitive to the number and distribution of internal nodes.For simple boundary value problems as shown in the numerical examples presented in this paper,no internal nodes are even needed to obtain moderately accurate numerical results.However,for some problems with stress concentration such as the fracture mechanics problems,some internal nodes may be necessary around the large displacement gradient region.There are two ways to determine the support size of the RBF.One is directly to specify the support size and the other is to specify the number of application points along each https://www.360docs.net/doc/d82453014.html,putational experience shows that the latter can give more stable computational results and four to six points in each direction are appropriate.

Acknowledgments

The support of the Natural Science Foundation of Jiangsu Province,PR China (no.:BK2006091)is gratefully acknowledged.References

[1]Gao XW,Davies TG.Boundary element programming in mechanics.

Cambridge,London,New York:Cambridge University Press;2002.[2]Banerjee PK.The boundary element methods in engineering.

London:McGraw-Hill Book Co;1994.

[3]Nardini D,Brebbia CA.A new approach for free vibration analysis

using boundary elements.In:Brebbia CA,editor.Boundary element methods in engineering.Berlin:Springer;1982.p.312–26.

[4]Albuquerque EL,Sollero P,Paiva WP.The radial integration method

applied to dynamic problems of anisotropic https://www.360docs.net/doc/d82453014.html,mun Numer Meth Eng 2007;23:805–18.

[5]Gao XW.A boundary element method without internal cells for two-dimensional and three-dimensional elastoplastic problems.ASME J Appl Mech 2002;69:154–60.

[6]Gao XW.The radial integration method for evaluation of domain

integrals with boundary-only discretization.Eng Anal Boundary Elements 2002;26:905–16.

[7]Gao XW.Boundary element analysis in thermoelasticity with and

without internal cells.Int J Numer Meth Eng 2003;57:975–90.

[8]Dong CY,Lo SH,Cheung YK.Numerical solution for elastic

inclusion problems by domain integral equation with integration by means of radial basis functions.Eng Anal Boundary Elements 2004;28:623–32.

[9]Gun H.An effective BE algorithm for 3D elastoplastic frictional

contact problems.Eng Anal Boundary Elements 2004;28:859–67.[10]Gun H.3D boundary element analysis of creep continuum damage

mechanics problems.Eng Anal Boundary Elements 2005;29:749–55.[11]Hematiyan MR.A general method for evaluation of 2D and 3D

domain integrals without domain discretization and its application in https://www.360docs.net/doc/d82453014.html,putational Mechanics 2007;39:509–20.

[12]Gao XW.A meshless BEM for isotropic heat conduction problems

with heat generation and spatially varying conductivity.Int J Numer Meth Eng 2006;66:1411–31.

[13]Gao XW,Zhang Ch,Sladek J,Sladek V.An ef?cient meshless BEM

for 2D nonhomogenoeus elastic solids.In:Sladek J,Sladek V,Atluri SN,editors.Advances in meshless methods.Tech Science Press;2006.[14]Sladek V,Sladek J,editors.Singular integrals in boundary element

methods.Southampton:Computational Mechanics Publications;1998.

[15]Gao XW.Numerical evaluation of two-dimensional singular

boundary integrals-theory and Fortran code.J Comput Appl Math 2006;188:44–64.

[16]Golberg MA,Chen CS,Bowman H.Some recent results and

proposals for the use of radial basis functions in the BEM.Eng Anal Boundary Elements 1999;23:285–96.

[17]Atluri SN,Shen S.The meshless local PetrovGalerkin (MLPG)

method.Tech Science Press;2002.

[18]Suresh S,Mortensen A.Fundamentals of functionally graded

materials.London:Institute of Materials;1998.

[19]Miyamoto Y,Kaysser WA,Rabin BH,Kawasaki A,Ford RG.

Functionally graded materials,design,processing and applications.Dordrecht:Kluwer Academic Publishers;1999.

[20]Lubliner J.Plasticity theory.New York:Macmillan Publishing

Company;1990.

[21]Sladek J,Sladek V,Zhang Ch.Meshless local Petrov Galerkin

method for stress analysis in functionally graded materials.Building Res J 2005;53:71–91.

1.E-05

1.E-041.E-031.E-021.E-011.E+001.E+01

Number of iterations

R e s i d u e R

Fig.14.Iteration history.

X.W.Gao et al./Engineering Analysis with Boundary Elements 31(2007)974–982

982