A meshless local natural neighbour interpolation method to modeling of functionally graded viscoelastic materials

Engineering Analysis with Boundary Elements 52 (2015) 92–98

Contents lists available at ScienceDirect

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

A meshless local natural neighbour interpolation method to modeling of functionally graded viscoelastic materials S.S. Chen n, C.J. Xu, G.S. Tong School of Civil Engineering and Architecture, East China Jiaotong University, NanChang 330013, China

art ic l e i nf o

a b s t r a c t

Article history: Received 26 July 2014 Received in revised form 9 November 2014 Accepted 13 November 2014

A meshless local natural neighbour interpolation (MLNNI) method applied to solve two-dimensional quasi-static and transient dynamic problems in continuously heterogeneous and linear viscoelastic media is presented and discussed in this article. The analysis is performed using the correspondence principle and the Laplace transform technique. In the present method, nodal points are spread on the analyzed domain and each node is surrounded by a polygonal sub-domain. The trial functions are constructed by the natural neighbour interpolation and the three-node triangular FEM shape functions are taken as test functions. The natural neighbour interpolants are strictly linear between adjacent nodes on the boundary of the convex hull. As a result, the essential boundary conditions can be imposed by directly substituting the corresponding terms in the system of equations. To get the viscoelastic solution in the time domain, the inverse Laplace transform algorithm of Stehfest is employed. Some numerical examples are given at the end to demonstrate the availability and accuracy of the present method. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Meshless method MLPG Viscoelasticity Correspondence principle

1. Introduction The interest in functionally graded materials (FGMs) has been steadily growing in recent years due to their continuously and smoothly varying material properties. If properly designed, FGMs can offer various advantages such as reduction of thermal stresses, minimization of stress concentration or intensity factors and attenuation of stress waves [1]. Therefore, they have gained potential applications in a wide range of engineering structures and components such as electronic devices, corrosion-resistant and wearresistant coatings, thermal barrier coatings and biomaterials [2,3]. Therefore, the behavior of FGMs needs to be better understood to fully exploit their characteristics. So far, extensive research has been carried out to efficiently and accurately simulate FGMs. Some representative examples are the works by Zhang and Paulino [4], Arciniega and Reddy [5], Song and Paulino [6], Sladek et al. [7], Wang and Qin [8] and so on. However, most of the previous research has been limited to elastic material behavior. As a matter of fact, FGMs exhibit creep and relaxation behavior [9] and therefore accurate simulation of them necessitates the use of viscoelastic constitutive models. In the framework of linear continuum theory, such a behavior can be described by linear viscoelasticity. Generally speaking, the correspondence principle [10] is one of the most useful tools in viscoelasticity because the Laplace transform of a viscoelastic

n

Corresponding author. Tel.: þ 86 791 87046084. E-mail address: [email protected] (S.S. Chen).

http://dx.doi.org/10.1016/j.enganabound.2014.11.016 0955-7997/& 2014 Elsevier Ltd. All rights reserved.

solution can be directly obtained from the corresponding elastic solution. After solving the transformed problem, the solutions are transformed back to the time domain by numerical methods. Unfortunately, for general linear viscoelastic FGMs, the correspondence principle may not be applicable. To avoid this problem, Paulino and Jin [11] proved that the correspondence principle may be used to obtain viscoelastic solutions only for a class of FGMs when the relaxation moduli are separable in space and time. Based on such particular correspondence principle for FGMs, Paulino and Jin [12,13] have subsequently studied crack problems of FGM strips subjected to antiplane shear loading and later to in-plane loading [14]. Besides, Sladek et al. [15] presented a meshless local Petrov-Galerkin method for static and dynamic analysis of continuously nonhomogeneous and linear viscoelastic solids. In the article, the transformed meshless formulations for solving nonhomogeneous viscoelastic problems are very similar to those for solving nonhomogeneous elastic problems. Although the finite element method (FEM) [16,17] and the boundary element method (BEM) [18] have been successfully established and applied to a variety of engineering problems, there is still a growing interest in developing new advanced numerical methods [19–24]. In recent years, meshless methods have been drawing much attention from researchers due to their higher adaptivity and lower cost in preparing input data for numerical analysis. As a fundamental base for the derivation of many meshless formulations, the meshless local Petrov-Galerkin (MLPG) method [25] differs significantly with other meshless methods in that the governing equation is satisfied node-by-node in a local way of making weighted residual zero over a sub-domain. The MLPG method provides the flexibility in choosing

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

the trial and test functions, as well as the sizes and shapes of local sub-domains and therefore brings forward a series of new methods [26–31] with amazing flexibility and efficiency. Among these, the meshless local natural neighbour interpolation (MLNNI) method [30,31] is formulated to combine the advantage of easy imposition of essential boundary conditions of natural element method (NEM) [32] with some prominent features of the MLPG method. Remarkable successes of the MLNNI method have been reported in a wide range of computational problems [33–37]. In this article, the MLNNI method is applied to solve linear viscoelastic problems of FGMs by incorporating the correspondence principle. In our meshless method, a set of properly scattered nodes within the problem domain is employed to represent the continuous problem and no elements or meshes are required. Each node is surrounded by a polygonal sub-domain, which can be easily constructed with Delaunay tessellations. The viscoelastic solids can be effectively treated in the Laplace domain by means of the MLNNI method. Then, the inverse Laplace transform algorithm of Stehfest [38] is applied to obtain the final time-dependent solutions. It is worth mentioning that several quasi-static boundary value problems are required to be solved for various values of the Laplace transform parameter. Finally, several numerical examples are presented to demonstrate the validity and accuracy of the proposed method.

2. Governing equations Some of the governing equations for viscoelastic problems of FGMs are outlined in this section. Let us consider a transient dynamic problem in a continuously heterogeneous and linear viscoelastic domain Ω bounded by a surface Γ with an outward unit normal vector with the components ni . The elastodynamic equilibrium equation can be written as σ ij;j ðx; tÞ  ρðxÞu€ i ðx; tÞ þ bi ðx; tÞ ¼ 0

ð1Þ

where σ ij ðx; tÞ is the stress tensor, bi ðx; tÞis the body force vector, ρðxÞ is the mass density and ui ðx; tÞ is the displacement vector. It should be mentioned here that a comma denotes partial derivatives with respect to spatial variables and the conventional Einstein's summation rule over repeated indices is applied throughout the analysis. All material parameters could be dependent on spatial variables. The quasi-static problems can be considered formally as a special case of the general transient problems by omitting the acceleration term u€ i ðx; tÞin the equilibrium eq. (1). Consequently, both cases are analyzed simultaneously in this article. The kinematic relations for the strain field εij is defined by εij ¼

1 ðu þ uj;i Þ 2 i;j

ð2Þ

and the linear viscoelastic constitutive law Z t deij dτ μðx; t τÞ sij ðx; tÞ ¼ 2 dτ 0 Z

t

ð3aÞ

dεkk dτ dτ

ð3bÞ

1 1 sij ¼ σ ij  σ kk δij ; eij ¼ εij  εkk δij 3 3

ð4Þ

σ kk ðx; tÞ ¼ 3

0

Kðx; t  τÞ

93

The corresponding boundary and initial conditions are given as follows ui ðx; tÞ ¼ u~ i ðx; tÞ; on the essential boundary Γ u

ð5aÞ

t i ðx; tÞ ¼ σ ij ðx; tÞnj ðxÞ ¼ t~ i ðx; tÞ; on the natural boundary Γ t

ð5bÞ

ui ðx; tÞjt ¼ 0 ¼ ui ðx; 0Þ and u_ i ðx; tÞjt ¼ 0 ¼ u_ i ðx; 0Þ in Ω

ð5cÞ

where Γ u is the part of the global boundary with prescribed displacements, while on Γ t the traction vector is prescribed. The correspondence principle is probably the most useful tool in viscoelasticity because the Laplace transform of the viscoelastic solution can be directly obtained from the existing elastic solution. Unfortunately, in general, the correspondence principle does not hold for FGMs. To avoid this problem, for a class of FGMs with the following separable relaxation functions μðx; tÞ ¼ μ0 μðxÞf ~ ðtÞ

ð6aÞ

Kðx; tÞ ¼ K 0 K~ ðxÞgðtÞ

ð6bÞ

~ K~ ðxÞ,f ðtÞ and gðtÞ where μ0 and K 0 are material constants, and μðxÞ, are nondimensional functions. Paulino and Jin [11] have shown that the correspondence principle still holds. For simplicity, the only restrictive requirement employed in the present work is the separation of the spatial and temporal variables in the relaxation functions. Substituting Eqs. (6a)-(6b) into Eqs. (3a)-(3b) results in Z t deij sij ðx; tÞ ¼ 2μ0 μðxÞ dτ ð7aÞ ~ f ðt  τÞ dτ 0 σ kk ðx; tÞ ¼ 3K 0 K~ ðxÞ

Z

t 0

gðt  τÞ

dεkk dτ dτ

ð7bÞ

Applying the Laplace transform to Eqs. (1), (2), (7a), (7b) and boundary conditions (4a)-(4b), one obtains σ ij;j ðx; pÞ  ρðxÞp2 ui ðx; pÞ þ F i ðx; pÞ ¼ 0

ð8aÞ

1 εij ðx; pÞ ¼ ðui;j ðx; pÞ þ uj;i ðx; pÞÞ 2

ð8bÞ

sij ðx; pÞ ¼ 2μeij ðx; pÞ ¼ 2μ0 μðxÞpf ~ ðpÞeij ðx; pÞ

ð8cÞ

σ kk ðx; pÞ ¼ 3Kεkk ðx; pÞ ¼ 3K 0 K~ ðxÞpgðpÞεkk ðx; pÞ

ð8dÞ

ui ðx; pÞ ¼ u~ i ðx; pÞ

ð8eÞ

t i ðx; pÞ ¼ t~ i ðx; pÞ

on Γ u on Γ t

ð8fÞ

where a bar over a variable means the Laplace transform, p is the Laplace transform parameter and F i ðx; pÞ ¼ bi ðx; pÞ þ ρðxÞpui ðx; 0Þ þ ρðxÞu_ i ðx; 0Þ

ð9Þ

is the redefined body force in the Laplace domain. Notice that the set of Eqs. (8a)-(8f) has identical form as the Laplace transformed elastodynamic equations of FGMs with the time-independent shear modulus μ ¼ μ0 μðxÞ ~ and the bulk modulus K ¼ K 0 K~ ðxÞ, provided that the transformed viscoelastic quantities are associated with the corresponding transformed elastodynamic quantities and μ0 pf ðpÞ and K 0 pgðpÞ are associated with μ0 and K 0 , respectively [15].

with

where sij and eij are deviatoric components of the stress and strain tensors, respectively, δij is the Kronecker delta, and μðx; tÞ and Kðx; tÞ are the relaxation functions in shear and dilatation, respectively. Note that the relaxation functions μðx; tÞ and Kðx; tÞ depend on spatial coordinates, whereas in homogeneous viscoelasticity they are only functions of time, that is, μ ¼ μðtÞ and K ¼ KðtÞ[10].

3. Implementation of the MLNNI method 3.1. Brief description of the natural neighbour interpolation This section gives a brief summary of the natural neighbour interpolation (NNI), of which excellent illustrations can be referred to Sukumar et al. [32]. The NNI is one of many possible schemes for an interpolation of discrete data with reasonable accuracy, which is

94

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

based on the Voronoi diagram and its dual Delaunay tessellationof  the domain. Consider a set of nodes N ¼ n1 ; n2 ; ⋯; nM in two-dimensional Euclidean space R2 . The first-order Voronoi diagram of the set N is a partition of the plane into regions T I , where each region T I is associated with a node nI , such that any point x in T I is closer to nI (nearest neighbour) than to any other nodes, in mathematical terms T I ¼ fx A R2 : dðx; xI Þ o dðx; xJ Þ

8 J a Ig

ð10Þ

where dðx; xI Þ is the distance between x and xI . For the definition of Sibson interpolation, it is necessary to previously introduce the concept of second-order Voronoi cell. This second-order diagram is a subdivision of cells T IJ , where each region is defined as the locus of the points that have the node nI as the closet node and the node nJ as the second closet node. In mathematical terms, the second-order Voronoi cell T IJ is defined as 2

T IJ ¼ fx A R : dðx; xI Þ o dðx; xJ Þ o dðx; xK Þ

8 J a I a Kg

ð11Þ

It should be noted that T IJ is non-empty only if nI and nJ are natural neighbours. The natural neighbour shape functions of x with respect to one of its natural neighbours I are defined in two dimensions as the ratio of the area of T xI and T x ϕI ðxÞ ¼ AI ðxÞ=AðxÞ

ð12Þ

where n

AðxÞ ¼ ∑ AJ ðxÞ

ð13Þ

J¼1

and n is the number of natural neighbours of point x. From this definition, and in the context of two-dimensional approximations, the unknown function uðxÞ is approximated in the form n

uh ðxÞ ¼ ∑ ϕI ðxÞuI

ð14Þ

I¼1

where uI ðI ¼ 1; 2; ⋯; nÞ are the vectors of nodal values at the n natural neighbours. It should be stressed here that the natural neighbour (Sibson) shape functions possess some remarkable properties such as positivity, interpolation, partition of unity and linear completeness. 3.2. Weak formulation and discretization

bounded by Γ Is can be written as Z vI ðσ ij;j ðx; pÞ  ρðxÞp2 ui ðx; pÞ þ F i ðx; pÞÞdΩ ¼ 0

where vI is the test function. As the essential boundary conditions can be imposed directly as in the FEM, there is no need to introduce Lagrange multipliers or penalty parameters in Eq. (15). Using Gaussian divergence theorem and considering t i ðx; pÞ ¼ σ ij ðx; pÞnj ðxÞ

ð16Þ

the local weak form (15) can be written as R R R Γ Isi vI t i ðx; pÞdΓ þ Γ Isu vI t i ðx; pÞdΓ þ Γ Ist vI t i ðx; pÞdΓ  Z ðvI;j σ ij ðx; pÞ þ vI ρðxÞp2 ui ðx; pÞ  vI F i ðx; pÞÞdΩ ¼ 0

ð17Þ

ΩIs

where Γ Isi [ Γ Isu [ Γ Ist ¼ Γ Is . Here,Γ Isi is the local boundary that is totally inside the global domain,Γ Isu ¼ Γ Is \ Γ u is the part of the local boundary which lies on the global boundary with prescribed displacements, and similarly, Γ Ist ¼ Γ Is \ Γ t is the part of the local boundary that lies on the global boundary with prescribed tractions. In the MLNNI method, the domain is subdivided into Delaunay tessellations. Therefore, the local polygonal sub-domain ΩIs can be conveniently constructed by collecting all the surrounding Delaunay triangles T iI with node I being their common vertices, as shown in Fig. 1. To simplify Eq. (17), we can deliberately select the test function vI such that they vanish over Γ Isi . This can be easily accomplished in the present work by choosing the test functions vI to be the three-node triangular FEM shape functions N I in each Delaunay triangle T iI . By doing so, Eq. (17) can be simplified as R R R ~ Γ Isu N I t i ðx; pÞdΓ þ Γ Ist N I t i ðx; pÞdΓ þ ΩIs N I F i ðx; pÞdΩ Z ðN I;j σ ij ðx; pÞ þ N I ρðxÞp2 ui ðx; pÞÞdΩ ¼ 0 ð18Þ  ΩIs

Obviously, there are two cases as shown in Fig. 2 for the integrals over Γ Isu in Eq. (18). For Fig. 2(a), the test functions NI are equal to zero over Γ Isu , and thus the integrals over Γ Isu are equal to zero; For Fig. 2(b), as already pointed out, essential boundary conditions must be enforced on node I and therefore the integrals over Γ Isu require no computing. Consequently, Eq. (18) can be further simplified as Z Z Z ðN I;j σ ij ðx; pÞ þ N I ρðxÞp2 ui ðx; pÞÞdΩ ¼ N I t~i ðx; pÞdΓ þ N I F i ðx; pÞdΩ ΩIs

Instead of writing the global Galerkin weak form for the governing equations presented in section 2, the MLNNI method is derived from a local weak form over a set of overlapping subdomains ΩIs ðI ¼ 1; 2; ⋯; MÞ, where I indicates a node, and M is the total number of nodes. Using the local weighted residual technique, the local weak form of Eq. (8a) over a local sub-domain ΩIs

ð15Þ

ΩIs

Γ Ist

ΩIs

ð19Þ Assume that the problem domain Ω is represented by properly scattered nodes. Then the Laplace transformed displacements in Eq. (19) can be approximated using Eq. (14). Substitution of Eq. (14) into Eq. (19) for all nodes leads to the following discretized system equations KðpÞuðpÞ ¼ fðpÞ

ð20Þ

where uðpÞ is the nodal vectors of Laplace transformed displacements, and KðpÞ and fðpÞ are the matrices of stiffness and load, respectively. They are in the following forms Z Z KIJ ðpÞ ¼ VTI Dðx; pÞBJ dΩ þ N I ρðxÞp2 ΦJ dΩ ð21aÞ ΩIs

Fig. 1. Local polygonal sub-domain and boundary for node I.

ΩIs

Fig. 2. Essential boundary condition Γ Isu over sub-domain ΩIs .

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

Z f I ðpÞ ¼

Γ Ist

~ pÞdΓ þ N I tðx;

Z ΩIs

N I Fðx; pÞdΩ

ð21bÞ

y

105

85

where the matrices of VI ,ΦJ ,BJ and Dðx; pÞ are given as follows 2 3 N I;x 0 6 0 N I;y 7 VI ¼ 4 ð22aÞ 5 NI;y N I;x " ΦJ ¼

ϕJ

0

0

ϕJ

2

ϕJ;x 6 BJ ¼ 4 0 ϕJ;y

P

21

L

ð22bÞ

x

Fig. 3. A viscoelastic strip subjected to a uniform tension.

3

0 ϕJ;y 7 5

ð22cÞ

4.1. A viscoelastic strip subjected to a uniform loading

ϕJ;x

E 1  v2

1

v

6v 4 0

1 0

0

3

0 7 5 for plane stress condition

ð22dÞ

1v 2

with v¼

w

1

#

2 Dðx; pÞ ¼

95

3K  2μ

ð23aÞ

6K þ 2μ

E ¼ 2μð1 þvÞ

ð23bÞ

here, E and v are the Laplace transform of Young's modulus and Poisson's ratio, respectively. Although nonsymmetric stiffness matrices may require more computational time and memory storage, the flexibility, ease of implementation, and accuracy embodied in the present formulation will still make it attractive.

Once the nodal unknowns have been determined numerically by solving Eq. (20), the inversion of Laplace transform is another important issue. In the present analysis, the inverse Laplace transform algorithm of Stehfest [38] is employed, which is easy to implement and leads to accurate results for many problems. Accordingly, the approximation of the time domain solution can be obtained by   ln 2 N ln 2 f ðtÞ ¼ i ð24Þ ∑ V if t i¼1 t where V i is given by the following equation N=2

k ð2kÞ! k ¼ ½ði þ 1Þ=2 ðN=2  kÞ!k!ðk  1Þ!ði kÞ!ð2k  iÞ! minði;N=2Þ

∑

E ¼ E0 f ðtÞ ¼ E0 ½E1 =E0 þð1 E1 =E0 Þexpð t=t 0 Þ

ð25Þ

The accuracy and the computational effort of the method are dependent on the truncation number Nof the series in Eq. (24), which should be optimized by trial and error. Increasing N increases the accuracy of the result up to a point, and then the accuracy declines because of increasing round-off errors. An optimal choice of 10 r N r20 has been recommended by Wang [39] and N ¼ 16 is adopted in the present work.

ð26Þ

where E1 =E0 ¼ 0:5,E0 ¼ 1,t 0 ¼ 1and Poisson's ratio v ¼ 0:0. According to the correspondence principle, the Young's modulus E0 is replaced by E0 pf ðpÞ ¼ E1 þ ðE0  E1 Þ

3.3. Numerical inversion of Laplace transform

V i ¼ ð  1ÞN=2 þ i

A viscoelastic strip, as depicted in Fig. 3, is considered when subjected to a uniform traction loading P ¼ 1:0HðtÞ at the end x ¼ L, wherein HðtÞ is the Heaviside unit step function. The opposite end of the strip is fixed in the x-direction, while its other boundaries are free of tractions. The length of the strip is L ¼ 3 and the width w ¼ 1. As its analytical solution is very simple and easily obtained, this example is recommended to verify the accuracy of the numerical model proposed in this work. A standard linear viscoelastic solid with homogeneous material and constant relaxation time is first considered in our numerical analysis

pt 0 pt 0 þ 1

ð27Þ

Under the quasi-static assumption, the analytical solution for the longitudinal displacement at the loaded end of the viscoelastic strip can be obtained as follows       Lðpt 0 þ 1Þ ðE1  E0 Þ E1 t 1 u ¼ L1 exp  ¼L þ ð28Þ 2 E1 E0 E0 t 0 E1 E1 p þ E0 p t 0 in which L  1 represents the inverse Laplace transform. In the numerical computation by the MLNNI method, regularly distributed nodes 21  5 are employed here, which is also shown in Fig. 3. To study the influence of number of Gaussian point on accuracy, six Gaussian points, as well as three Gaussian points, are used for domain integrals in each Delaunay triangular region. The numerical results for the longitudinal displacement at node 21 are illustrated and compared with analytical results from Eq. (28) in Fig. 4. As it may be readily deduced from Fig. 4, the present method produces a result indistinguishable from the analytical solution. It is also observed that the number of Gaussian point has little influence on the numerical accuracy. For the consideration of computational inefficiency, we recommend using three Gaussian points in each Delaunay triangular region for numerical integration. In the next stage, nonhomogeneous material properties and quasi-static assumption are considered. The material parameters are described by

4. Numerical examples

E ¼ E0 expðγxÞ½E1 =E0 þ ð1  E1 =E0 Þexpð  t=t 0 Þ

The present MLNNI method has been coded in Fortran. Case studies are carried out in order to examine the performance of the proposed solution procedure for linear viscoelastic problems of functionally graded materials. Except specially mentioned, three Gaussian points are employed in each Delaunay triangular region for domain integration.

where E1 =E0 ¼ 0:5,E0 ¼ 1,t 0 ¼ 1,v ¼ 0:0 and γ is known as the material gradient parameter. By increasing the exponent γ, the Young's modulus increase, and in turn, the longitudinal displacements of the functionally graded viscoelastic strip decrease. In such case, the analytical solution for the longitudinal displacement at the loaded end can be obtained by using the correspondence principle

ð29Þ

96

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

6.5

5

6.0

4 Displacement

Displacement

5.5 5.0 4.5

Analytical Present-3GP Present-6GP

4.0 3.5 3.0

3 2 Sladek et al. [15] Present

1

2.5

0

2.0 0

2

4

6

8

10

0

2

Time Fig. 4. Longitudinal displacement in a homogeneous viscoelastic strip (quasi-static analysis).

4

6

8 Time

10

12

14

16

Fig. 6. Transient response of a homogeneous viscoelastic strip (dynamic analysis).

Petrov-Galerkin method. In Fig. 6, numerical results for the time variation of the longitudinal displacement at node 21 are compared with those of Sladek et al. [15]. One can find from Fig. 6 that our results are, in general, close to the available numerical results. 4.2. A viscoelastic hollow cylinder subjected to internal pressure The second example regards a thick hollow cylinder of inner radius a ¼ 8 and outer radius b ¼ 16 subjected to an internal pressure P ¼ P 0 HðtÞwith P 0 ¼ 1, as shown in Fig. 7. Here HðtÞis the Heaviside unit step function. The cylinder is assumed to be sufficiently long such that a plane strain condition assumption holds. For comparison purposes, the same homogeneous material properties as used by Lee and Westmann [41] are assumed in the first step. The viscoelastic material is characterized by an elastic bulk modulus and a shear relaxation modulus as follows Fig. 5. Longitudinal displacement in a functionally graded viscoelastic strip (quasistatic analysis).

KðtÞ ¼ K 0 gðtÞ ¼ K 0

ð32aÞ

as follows

  G1 expð  t=t 0 Þ GðtÞ ¼ G0 f ðtÞ ¼ G0 1 þ G0

ð32bÞ

    1  expð  γLÞ ðE1  E0 Þ E1 t 1 u¼ exp  þ γ E1 E0 E0 t 0 E1

ð30Þ

For three different material gradient parameters γ ¼ 0:2, 0.4 and 0.6, numerical results for the time variation of the longitudinal displacement at node 21 are illustrated in Fig. 5. For comparison, the corresponding analytical results from Eq. (30) are also plotted in the same graph. From this figure, one can observe that the computational results of the present MLNNI method are in good agreement with the analytical results. Now the influence of the mass density or the inertial term on the time variation of the longitudinal displacement is investigated, where the inertial term in the governing Eq. (1) has to be considered. Here, the homogeneous material properties are assumed and the same geometry of the strip is considered as in the quasi-static case. The generalized three parameters viscoelastic model [40] employed here can be written as   α α d E d q1 α sij þ sij ¼ eij þ q2 α eij ð31aÞ 1þv dt dt   α α d E d εkk þ q2 α εkk q1 α σ kk þ σ kk ¼ 1  2v dt dt

ð31bÞ

where E ¼ 1,v ¼ 0,ρ ¼ 1,q1 ¼ 0:5,q2 ¼ 1 and α ¼ 1:0. The same problem was also analyzed by Sladek et al. [15] using the meshless local

where G0 ¼ 120,G1 ¼ 360,K 0 ¼ 1280 and t 0 ¼ 2:5. Due to the symmetry, only the upper right quadrant of the cylinder has to be discretized. In order to study the influence of the irregular node distribution on the accuracy of the present method, both the regular and irregular node distribution shown in Fig. 8 are employed. Numerical results for the radial displacements are plotted in Fig. 9. It is evident that the results from the regular and the irregular node distributions almost coincide with each other, and the results from both node distributions agree well with those of Lee and Westmann [41]. Obviously, the present method works well for linear viscoelastic problems with the help of the correspondence principle and the Laplace transform technique. In addition, we can conclude that the irregular node distribution does not affect the numerical accuracy in the present meshless method much. Now we proceed to consider a functionally graded viscoelastic cylinder under the quasi-static assumption. The material parameters are described by KðtÞ ¼ K 0 expðγðr  aÞÞ

ð33aÞ

GðtÞ ¼ ½G0 þ G1 expð  t=t 0 Þexpðγðr  aÞÞ

ð33bÞ

where G0 ¼ 120,G1 ¼ 360,K 0 ¼ 1280, t 0 ¼ 2:5 and the material gradient parameter γ is selected as γ ¼ 0:1 and 0.2. Figs. 10 and 11 depict the time variation of the radial displacements on inner and outer surfaces of the cylinder respectively. Also shown in these figures are

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

97

b a

P Fig. 9. Time variation of the radial displacement in a homogeneous viscoelastic hollow cylinder (quasi-static analysis).

Fig. 7. A viscoelastic hollow cylinder subjected to internal pressure.

Fig. 10. Time variation of the radial displacement on the inner surface of a functionally graded viscoelastic hollow cylinder (quasi-static analysis).

Fig. 11. Time variation of the radial displacement on the outer surface of a functionally graded viscoelastic hollow cylinder (quasi-static analysis).

5. Conclusions Fig. 8. Nodes for a viscoelastic hollow cylinder.

the numerical results of the homogeneous cylinder. One can observe from Figs. 10 and 11 that the increasing of the material gradient parameter of FGMs decreases the radial displacements.

In this article, we proposed analyzing the linear viscoelastic problems of FGMs by the MLNNI method together with the numerical inversion of the Laplace transform. Through the use of the correspondence principle, the MLNNI developed for the linear elastic solid can be applied directly to the linear viscoelastic solids in the Laplace domain. After obtaining the physical responses in the

98

S.S. Chen et al. / Engineering Analysis with Boundary Elements 52 (2015) 92–98

Laplace domain, their associated value in time domain are calculated by the inverse Laplace transform algorithm of Stehfest. In the present meshless method, the analyzed domain is divided into small overlapping polygonal sub-domains. The three-node triangular FEM shape functions and the NNI are chosen as the test and trial functions, respectively. Compared to the conventional MLPG method, the NNI has Kronecker delta function property, and the construction of local sub-domain is convenient both for internal and boundary nodes. More importantly, the present meshless method seems to be more flexible than the standard FEM due to the fact that an adaptation of the nodal density is easier than a mesh adaptation. All the numerical examples given above successfully demonstrate that the proposed solution technique is quite efficient. Acknowledgements This article is supported by the Doctoral Startup Foundation of East China Jiaotong University (Grant No. 09132020) and the Foundation of Jiangxi Educational Committee (Grant No. KJLD14041). References [1] Ching HK, Yen SC. Transient thermoelastic deformations of 2-D functionally graded beams under nonuniformly convective heat supply. Compos Struct 2006;73:381–93. [2] Oonishi H, Noda T, Ito S, Kohda A, Yamamoto H, Tsuji E. Effect of hydroxyapatite coating on bone growth into porous titanium alloy implants under loaded conditions. J Biomater Appl 1994;5:23–7. [3] Sladek J, Sladek V, Zhang Ch. An advanced numerical method for computing elastodynamic fracture parameters in functionally graded materials. Comp Mater Sci 2005;32:532–43. [4] Zhang Z, Paulino GH. Wave propagation and dynamic analysis of smoothly graded heterogeneous continua using graded finite elements. Int J Solids Struct 2007;44:3601–26. [5] Arciniega RA, Reddy JN. Large deformation analysis of functionally graded shells. Int J Solids Struct 2007;44:2036–52. [6] Song SH, Paulino GH. Dynamic stress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. Int J Solids Struct 2006;43:4830–66. [7] Sladek J, Sladek V. Zhang Ch. Application of meshless local Petrov-Galerkin (MLPG) method to elastodynamic problems in continuously nonhomogeneous solids. CMES: Comput Model Eng Sci 2003;4(6):637–47. [8] Wang H, Qin QH. Meshless approach for thermo-mechanical analysis of functionally graded materials. Eng Anal Bound Elem 2008;32:704–12. [9] Paulino GH, Jin ZH, Dodds RH. Failure of functionally graded materials. In: Karihaloo B, Knauss WG, editors. Comprehensive structural integrity, vol. 2. Amsterdam: Elsevier; 2003. p. 607–44 ([chapter 13]). [10] Christensen RM. Theory of viscoelasticity. New York: Academia Press; 1971. [11] Paulino GH, Jin ZH. Correspondence principle in viscoelastic functionally graded materials. ASME J Appl Mech 2001;68:129–32. [12] Paulino GH, Jin ZH. Viscoelastic functionally graded materials subjected to antiplane shear fracture. ASME J Appl Mech 2001;68:284–93. [13] Paulino GH, Jin ZH. A crack in a viscoelastic functionally graded material layer embedded between two dissimilar homogeneous viscoelastic layers-antiplane shear analysis. Int J Fracture 2001;111:283–303. [14] Jin ZH, Paulino GH. A viscoelastic functionally graded strip containing a crack subjected to in-plane loading. Eng Fract Mech 2002;69:1769–90. [15] Sladek J, Sladek V, Ch Zhang, Schanz M. Meshless local Petrov-Galerkin method for continuously nonhomogeneous linear viscoelastic solids. Comput Mech 2006;37:279–89.

[16] Long YQ, Cen S, Long ZF. Advanced finite element method in structural engineering. Berlin, Heidelberg, Beijing: Springer-Verlag Gmbh, Tsinghua University Press; 2009. [17] Cen S, Fu XR, Zhou MJ. 8- and 12- node plane hybrid stress-function elements immune to severely distorted mesh containing elements with concave shapes. Comput Methods Appl Mech Eng 2011;200:2321–36. [18] Huang S, Liu YJ. A fast multipole boundary element method for solving the thin plate bending problem. Eng Anal Bound Elem 2013;37(6):967–76. [19] Hu DA, Long SY, Liu KY, Li GY. A modified meshless local Petrov-Galerkin method to elasticity problems in computer modeling and simulation. Eng Anal Bound Elem 2006;30:399–404. [20] Gu YT, Wang QX, Lam KY, Dai KY. A pseudo-elastic local meshless method for analysis of material nonlinear problems in solids. Eng Anal Bound Elem 2007;31:771–82. [21] Peng MJ, Li RX, Cheng YM. Analyzing three-dimensional viscoelasticity problems via the improved element-free Galerkin (IEFG) method. Eng Anal Bound Elem 2014;40:104–13. [22] Wang DD, Zhang HJ. A consistently coupled isogeometric-meshfree method. Comput Methods Appl Mech Eng 2014;268:843–70. [23] Ren HP, Cheng J, Huang AX. The complex variable interpolating moving leastsquares method. Appl Math Comput 2012;219:1724–36. [24] Sladek J, Sladek V, Ch Zhang, Krivacek J, Wen PH. Analysis of orthotropic thick plates by meshless local Petrov-Galerkin (MLPG) method. Int J Numer Meth Eng 2006;67:1830–50. [25] Atluri S N, Zhu T. A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [26] Barry W. A wachspress meshless local Petrov-Galerkin method. Eng Anal Bound Elem 2004;28(5):509–23. [27] Liu GR, Gu YT. A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246(1):29–46. [28] Xiao JR. Local Heaviside weighted MLPG meshless method for twodimensional solids using compactly supported radial basis functions. Comput Methods Appl Mech Eng 2004;193:117–38. [29] Li H, Wang QX, Lam KY. Development of a novel meshless local Kriging (LoKriging) method for structural dynamic analysis. Comput Methods Appl Mech Eng 2004;193:2599–619. [30] Cai YC, Zhu HH. A meshless local natural neighbour interpolation method for stress analysis of solids. Eng Anal Bound Elem 2004;28(6):607–13. [31] Wang K, Zhou SJ, Shan GJ. The natural neighbour Petrov-Galerkin method for elasto-statics. Int J Numer Meth Eng 2005;63:1126–45. [32] Sukumar N, Moran B, Belytschko T. The natural element method in solid mechanics. Int J Numer Meth Eng 1998;43:839–87. [33] Liu YH, Chen SS, Li J, Cen ZZ. A meshless local natural neighbour interpolation method applied to structural dynamic analysis. CMES: Comput Model Eng Sci 2008;31(3):145–56. [34] Chen SS, Li QH, Liu YH, Xia JT, Xue ZQ. Dynamic elastoplastic analysis using the meshless local natural neighbour interpolation method. Int J Comput Methods 2011;8(3):463–81. [35] Zhu HH, Liu WJ, Cai YC, Miao YB. A meshless local natural neighbour interpolation method for two-dimension incompressible large deformation analysis. Eng Anal Bound Elem 2007;31(10):856–62. [36] Li SL, Liu KY, Long SY, Li GY. The natural neighbour Petrov-Galerkin method for thick plates. Eng Anal Bound Elem 2011;35(4):616–22. [37] Chen SS, Li QH, Liu YH, Xue ZQ. A meshless local natural neighbour interpolation method for analysis of two-dimensional piezoelectric structures. Eng Anal Bound Elem 2013;37:273–9. [38] Stehfest H. Algorithm 368: numerical inversion of Laplace transform. Commun Assoc Comput Mach 1970;13:47–9. [39] Wang Y. A note on the solution to the diffusion equation. Int J Numer Anal Met 1996;20:443–52. [40] Schanz M, Antes H. A new visco- and elastodynamic time domain boundary element formulation. Comput Mech 1997;20:452–9. [41] Lee S, Westmann RA. Application of high-order quadrature rules to timedomain boundary element analysis of viscoelasticity. Int J Numer Meth Eng 1995;38:607–29.