Analysis of functionally graded plates by meshless method: A purely analytical formulation

Engineering Analysis with Boundary Elements 36 (2012) 639–650

Contents lists available at SciVerse ScienceDirect

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

Analysis of functionally graded plates by meshless method: A purely analytical formulation P.H. Wen a, M.H. Aliabadi b,n a b

School of Engineering and Material Sciences, Queen Mary, University of London, London E1 4NS, UK Department of Aeronautics, Imperial College, London SW7 2BY, UK

a r t i c l e i n f o

abstract

Article history: Received 22 August 2011 Accepted 29 November 2011 Available online 11 January 2012

Functionally graded plates under static and dynamic loads are investigated by the local integral equation method (LIEM) in this paper. Plate bending problem is described by the Reissner moderate thick plate theory. The governing equations for the functionally graded material with respect to the neutral plane are presented in the Laplace transform domain and therefore the in-plane and bending problems are uncoupled. Both isotropic and orthotropic material properties are considered. The local integral equation method is developed with the locally supported radial basis function (RBF) interpolation. As the closed forms of the local boundary integrals are obtained, there are no domain or boundary integrals to be calculated numerically in this approach. The solutions of the nodal values for the entire plate are obtained by solving a set of linear algebraic equation system with certain boundary conditions. Details of numerical procedures are presented and the accuracy and convergence characteristics of the method are examined. Several examples are presented for the functionally graded plates under static and dynamic loads and the accuracy for proposed method has been observed compared with 3D analytical solutions. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Functionally graded materials Kirchhoff and Reissner plates Laplace transform method Meshless method Radial basis function Static and dynamic loads

1. Introduction The boundary element method (BEM) is an effective numerical tool for the solution of boundary or initial-boundary value problems [1]. The first application of the boundary integral equation method to Reissner’s plate theory was given by Van der Ween [2]. Dynamic analysis of elastic Reissner plates has been performed by the direct BEM in the frequency domain [3,4]. Unfortunately, the number of applications to analyse functionally graded plates by the use of BEM is very limited due to the lack of availability of fundamental solutions. Meshfree methods have achieved remarkable progress in recent years. Typical methods include the method of fundamental solution [5,6], the diffuse element method [7], the element-free Galerkin method [8] and the meshless local Petrov–Galerkin method [9]. These methods have attracted much attention during the past decade, especially owing to their high adaptivity and low effort to prepare input data for numerical analyses [10]. The meshless collocation method, meshless local Petrov–Galerkin method and local boundary integral equation method (LBIE) with moving least square and radial basis function interpolations has been developed by Ferreira et al. [11], Gilhooley et al. [12], Sladek

n

Corresponding author. E-mail address: [email protected] (M.H. Aliabadi).

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

et al. [13] and Li et al. [14], respectively, for anisotropic non-homogeneous and composite plates. Recently an improved meshless collocation method has been proposed for two elastostatic and elastodynamic problems by Wen and Aliabadi [15] and for nonlinear analysis by Wen and Hon [16]. A comprehensive review of meshless methods can be found by Atluri [17] and Liu [18]. The interpolation methods including the radial basis function method and moving least squares method play an important role in the meshless investigations. Hardy [19] developed a general scattered data approximation method, called multi-quadric (MQ), to approximate two-dimensional geographical surfaces. Numerical results show that the MQ method offers a very accurate interpolation. In Franke’s review [20] paper, MQ was rated one of the best methods among 29 scattered data interpolation schemes, based on accuracy, stability, efficiency, memory requirement and ease of implementation. Kansa [21] derived a modified MQ scheme for solving PDE problems. In the last two decades, the developments of radial basis functions as a truly meshless method has drawn the attention of many investigators (see Golberg and Chen [22]). Functionally graded materials (FGMs) are multi-phase materials with the phase volume fractions varying gradually in space, in a pre-determined profile. This results in continuously graded mechanical properties that vary gradually with location within the material. Generally, the spatial gradients in material behaviour render FGMs as superior to conventional composites. FGMs

640

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

differ from composites wherein the volume fraction of the inclusion is uniform throughout the composite. The closest analogous of FGMs are laminated composites, but the latter possess distinct interfaces. More details about FGMs can be referred to references [23,24]. The most familiar FGMs are compositionally graded from a refractory ceramic to a metal. It can incorporate incompatible functions such as the heat, wear and corrosion resistances of ceramics and the high toughness, high strength and bonding capability of metals without severe internal thermal stresses. In this paper, the meshless local integral equation method is presented for the FGMs plate with Reissner’s plate theory. With the use of radial basis functions, the analytical solutions for all domain integrals in the weak form are derived. The weak formulations for the governing equations with a unit test function are obtained exactly for the local domain integrals. As the closed form of the local domain integrals are obtained, the computational time is reduced significantly. The Laplace transform technique is applied for dynamic problems. Numerical results for isotropic and orthotropic plates with various boundary conditions and static and Heaviside loadings are presented to illustrate the applicability of the proposed method.

2. Governing equations for FGM plates

where material parameter matrix 2 E1 n21 =e 0 E1 =e 6 E n =e E =e 0 2 6 2 12 6 0 0 G12 Dðx,y,z0 Þ ¼ 6 6 6 0 0 0 4 0 0 0

0

0

3

0 7 7 7 0 7 7 0 7 5 G23

0 0 G13 0

ð4Þ

in which e¼1  n12n21, E1 and E2 are the Young’s moduli referring to the axes x and y, respectively. G12, G13 and G23 are shear moduli and nij are Poisson’s ratios. There are two main assumptions for the functionally graded materials. The first assumption is that the material properties are graded along the plate thickness as Pðx,y,z0 Þ ¼ Pb ðx,yÞ þ ½Pt ðx,yÞPb ðx,yÞðz0 =hÞn

ð5Þ

and the section assumption is Pðx,y,z0 Þ ¼ Pb ðx,yÞexpðz0 d=hÞ,

d ¼ lnðP t =P b Þ

ð6Þ

where Pt and Pb denote the properties of the top and bottom faces of the plate, respectively, and n is a parameter that indictates the material variation profile. In addition, Possion’s ratios and mass density of the plate are assumed to be constants. Under pure bending, we have Z h sx ðx,y,z0 Þdz0 ¼ 0 ð7Þ 0

Consider a FGMs rectangular plate of uniform thickness h, length a and width b as shown in Fig. 1. The Cartesian coordinate system is introduced such that the bottom and top surfaces of plate is placed in the plane z0 ¼0 and z0 ¼h. Then, the spatial displacement field has the following form [1,2]:

Therefore, the position of the neutral plane is obtained for the first assumption of material variation profile:

u ¼ u0 þðz0 z0 Þw1 ðx,y,tÞ v ¼ v0 þðz0 z0 Þw2 ðx,y,tÞ w ¼ w3 ðx,y,tÞ

For the second assumption of variation profile, one holds   1 1 þ h ð9Þ z0 ¼ 1b ln b

ð1Þ

where z0 indicates the position of neutral plane where the in-plane stresses are zero, u0 and v0 displacements of in-plane, w1, w2 are rotations around x and y axes, respectively, and w3 is the out-of-plane deflection. The linear strains are given by

ex ¼ u0,1 þ ðz0 z0 Þw1,1 ey ¼ v0,2 þ ðz0 z0 Þw2,2 gxy ¼ u0,2 þ v0,1 þ ðz0 z0 Þðw1,2 þ w2,1 Þ gxz ¼ w1 þw3,1 gyz ¼ w2 þw3,2

ð2Þ

For a plane stress problem in orthotropic materials, the Hook’s law can be written in matrix form as 2 3 2 3

sx

ex

6 sy 7 6 ey 7 6 7 6 7 6 7 6 7 6 txy 7 ¼ Dðx,y,zÞ6 gxy 7 6 7 6 7 6t 7 6g 7 4 xz 5 4 xz 5

tyz

ð3Þ

gyz

z' b

y'

z0 ¼

ðn þ1Þðnb þ 2Þ h: 2ðn þ 2Þðnb þ 1Þ

where b ¼ Pb(x,y)/Pt(x,y). It should be noticed that the neutral plane does not exist on the middle plane (z0 ¼h/2) of the plate except b ¼1. New coordinate system is established on the neutral plane as shown in Fig. 1, i.e. z¼z0  z0, and the bending moments Mab and shear forces Qa are defined as [1] 2 3 2 3 " # M 11 Z hz0 " t # Z hz0 sx Q1 13 6M 7 6s 7 y k ¼ z dz, ¼ dz 4 22 5 4 5 Q2 t23 z0 z0 txy M 12 2 3 2 3 N11 Z hz0 sx 6 7 6s 7 ð10Þ and 4 N22 5 ¼ 4 y 5dz z0 txy N12 where k ¼5/6 in the Reissner plate theory. Substituting (2) and (3) into moment and force resultants (10) yields M11 ¼ D11 w1,1 þD12 w2,2 M22 ¼ D21 w1,1 þD22 w2,2 M12 ¼ Sðw1,2 þ w2,1 Þ

q ( x, y , t ) Q

z

M

z

h

ð8Þ

y M

x

a

M

x'

M

Q Fig. 1. A rectangular plate and coordinate system and sign convention of bending moments and shear forces for the plate on neutral plane.

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

Os

ð11Þ

in which D, S and C are material parameters. As far as the first material assumption is concerned, one has E1t f, e

D11 ¼

D12 ¼ D21 ¼

C 1 ¼ kG13t g, A11 ¼

E2t g , e

E1t n21 , e

C 2 ¼ kG23t g, A11 ¼

D22 ¼ E1t g , e

E2t f, e A11 ¼

B ¼ G12 g

S ¼ G12t f , E1t n21 g ¼ A21 , e

By ignoring coupling effect by inertia forces between in-plane and bending cases, one has the governing equations in the following form: @N 11 @N 12 þ þq1 ¼ J0 u€ 0 @x1 @x2 @N 12 @N 22 þ þq2 ¼ J0 v€ 0 @x1 @x2 @M 11 @M 12 €1 þ Q 1 ¼ I0 w @x1 @x2 @M 12 @M 22 €2 þ Q 2 ¼ I0 w @x1 @x2 @Q 1 @Q 2 €3 þ þq0 ¼ J 0 w @x1 @x2

ð14Þ

rh 3

ðhz0 Þ3 þ z30

i

0

Overall variables occurring in system equations, from Eq. (14), we have ~ 11 @M ~ 12 @M ~1 þ Q~ 1 ¼ I0 p2 w @x1 @x2 ~ 12 @M ~ 22 @M ~2 þ Q~ 2 ¼ I0 p2 w @x1 @x2 @Q~ 1 @Q~ 2 ~3 þ þ q~ 0 ¼ J0 p2 w @x1 @x2

and consider the equilibrium equation: Z Z Z ~ 3 q~ 0 ÞdO x1 ðQ~ 1 n1 þ Q~ 2 n2 ÞdG x1 ðJ 0 s2 w Q~ 1 dO ¼

ð20Þ

In the same way, one has Z Z ~ 3 q~ 0 ÞdO x2 ðQ~ 1 n1 þ Q~ 2 n2 ÞdG x2 ðJ 0 s2 w Q~ 2 dO ¼

ð21Þ

Os

Z

@Os

Os

@Os

Os

Therefore domain integrals in Eqs. (17a), (17b) and (17c) can be transformed into local boundary integrals and rewritten in a weak form in the Laplace transform space as Z ~ 11 n1 þ M ~ 12 n2 x1 ðQ~ 1 n1 þ Q~ 2 n2 ÞdG ½M @Os Z ~ 1 x1 ðJ 0 s2 w ~ 3 q~ 0 ÞdO ¼ ½I0 s2 w ð22aÞ Os

~ 12 n1 þ M ~ 22 n2 x2 ðQ~ 1 n1 þ Q~ 2 n2 ÞdG ½M Z ~ 2 x2 ðJ 0 s2 w ~ 3 q~ 0 ÞdO ¼ ½I0 s2 w

Z @Os

ðQ~ 1 n1 þ Q~ 2 n2 ÞdG ¼

Z Os

~ 3 q~ 0 ÞdO ðJ 0 s2 w

ð22cÞ

In the case of static circumstance, the local integral equations become Z Z ½M 11 n1 þ M 12 n2 x1 ðQ 1 n1 þQ 2 n2 ÞdG ¼ x1 q0 dO ð23aÞ @Os

Z Z @Os

ð16Þ

ð22bÞ

Os

@Os

Os

½M 12 n1 þ M 22 n2 x2 ðQ 1 n1 þQ 2 n2 ÞdG ¼

Z

Z ðQ 1 n1 þ Q 2 n2 ÞdG ¼  q0 dO

In the local boundary integral equation approaches, the weak form of differential equation over a local integral domain Os can be written as ! Z ~ 11 @M ~ 12 @M 2 ~ ~ þ Q 1 I0 p w 1 un dO ¼ 0 ð17aÞ @x1 @x2 Os

ð17bÞ

Os

x2 q0 dO

ð23bÞ

ð23cÞ

Os

The boundary conditions can be introduced directly and are given as ~ i¼w ~ 0i ði ¼ 1,2,3Þ on GD w ~ ij nj ¼ t~ 0 and ði ¼ 1,2Þ Q~ j nj ¼ t~ 0 on GT M i 3

3. Local integral equation method

! ~ 12 @M ~ 22 @M ~ 2 un dO ¼ 0 þ Q~ 2 I0 p2 w @x1 @x2

ð19Þ

@Os

and r is the mass density. The dots indicate differentiations with respect to time t. Therefore, the resultants of in-plane force are uncoupled with the moments and shear forces of bending (u0 ¼v0 ¼0 in this case). Applying the Laplace transform, Z 1 L½f ðx,y,z,tÞ ¼ f~ ðx,y,z,pÞ ¼ f ðx,y,z,tÞept dt ð15Þ

Os

By use of divergence theorem, one has  Z  Z @ @ ðx1 Q~ 1 Þ þ ðx1 Q~ 2 Þ dO ¼ x1 ðQ~ 1 n1 þ Q~ 2 n2 ÞdG @x2 Os @x1 @Os

Os

where

Z

where un is an arbitrary test function. The local weak forms in Eqs. (17a), (17b) and (17c) are a starting point to derive local boundary integral equations if appropriate test functions are selected. For a simple text function, all domain integrals in the local integral domain can be transformed to boundary integrals, which can be evaluated analytically. For example, a step functions can be used as the test functions uni in each integral domain: ( 1 at x A ðOs [ @Os Þ n u ðxÞ ¼ : ð18Þ 0 at x= 2Os

Z I0 ¼

ð17cÞ

ð12Þ

where t in the subscript indicates the property of material on the top surface of the plate and   nb þ 3 nb þ2 nb þ1 2 h z0 h , g ¼ h: ð13Þ f¼ 3ðn þ 3Þ 2ðn þ 2Þ nþ1

J 0 ¼ rh,

! @Q~ 1 @Q~ 2 2 ~ ~ þ þ q 0 J0 p w 3 un dO ¼ 0 @x1 @x2

Z

Q 1 ¼ C 1 ðw1 þ w3,1 Þ Q 2 ¼ C 2 ðw2 þ w3,2 Þ N11 ¼ A11 u0,1 þA12 v0,2 N22 ¼ A21 u0,1 þA22 v0,2 N12 ¼ Bðu0,2 þv0,1 Þ

641

ð24Þ

0 u~ 0i andt~ i

are the prescribed displacements and tractions, in which respectively, on the displacement boundary GD and on the traction boundary GT, and ni is the unit normal outward to the boundary G.

4. Radial basis function approximation Consider a local domain qOs shown in Fig. 2, which is the neighbourhood of a point x( ¼{x1,x2}) and is considered as the

642

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

node x

support domain of x

Local integral domain Ω

Node in support domain ξ

Fig. 2. Arbitrary distributed node, support domain of x, local integral domain for weak formulation.

domain of definition of the RBF approximation for the trial function at x and also called as support domain to an arbitrary point x. Generally the support domain is chosen as a circle centred at x. To interpolate the distribution of function w in the local domain qOs over a number of randomly distributed nodes n[¼ {n1, n2, y, nK}, nk ¼ (xk1, xk2), k¼1,2, y, K], the approximation of function u at the point x can be expressed by wðxÞ ¼

K X

Rk ðx, nk Þak ¼ RðxÞaðxÞ

ð25Þ

k¼1

where R(x) ¼{R1(x,n), R2(x,n), y, RK(x,n)} is the set of radial basis functions centred around the point x[¼(x1,x2)], fak gKk ¼ 1 are the unknown coefficients to be determined. The radial basis function selected multi-quadrics [20,21] as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð26Þ Rk ðx, nk Þ ¼ l þ ðx1 xk1 Þ2 þ ðx2 xk2 Þ2 where l is a free parameter. In order to guarantee unique solution of the interpolation problem, the displacement field can be interpolated by wðxÞ ¼

K X

Rk ðx, nÞak þ

T X

P t ðxÞbt ¼ RðxÞa þ PðxÞb

ð27Þ

t¼1

k¼1

and 2

P1 ðn1 Þ 6 P 6 1 ðn2 Þ 6 6 : 6 P0 ðnÞ ¼ 6 : 6 6 6 : 4 P1 ðnK Þ

P2 ðn1 Þ

...

P2 ðn2 Þ

...

:

...

:

...

:

...

P 2 ðnK Þ

...

P T ðn1 Þ

3

7 P T ðn2 Þ 7 7 7 : 7 7 : 7 7 7 : 5 PT ðnK Þ

ð32Þ

KT

Solving these equations in Eq.(30) gives a ¼ aw

b ¼ bw,

1 b ¼ ðPT0 R1 P0 T R1 0 P0 Þ 0 ,

T 1 1 a ¼ R1 P0 T R1 ð33Þ 0 ½IP0 ðP0 R0 P0 Þ 0 

where I denotes the diagonal unit matrix. Substituting the coefficients a and b into Eq. (27), we can obtain the approximation of the field function in terms of the nodal values:

UðxÞ ¼ RðxÞa þ PðxÞb,

wðxÞ ¼ UðxÞw,

ð34Þ

in which U(x) [¼(j1,j2, y, jK)] is defined as shape function, matrix R(x) and P(x) are scale (1  K) and (1  6) matrix, respectively. It is worth noting that the shape function depends uniquely on the distribution of scattered nodes within the support domain and it has the Kronecker Delta property.

along with the constraints: K X

5. Analytical forms of local boundary integrals P t ðnk Þak ¼ 0,

1 r t rT

ð28Þ

k¼1

where fP t gTt ¼ 1 is a basis for PT  1, the set of d-variate polynomials of degree rT 1. In this paper, following polynomials are considered (T ¼6): P ¼ f1,x1 ,x2 ,x21 ,x1 x2 ,x22 g

ð29Þ

Firstly, consider a unit test function, i.e. ji(x)¼ 1 and the local domain is enclosed by several straight lines as shown in Fig. 3, therefore, the local boundary integral equation becomes for the collocation point in the domain xj[¼(xj1,xj2)]: Z h i ~ 12 n2 x1 ðQ~ 1 n1 þ Q~ 2 n2 Þ dG ~ 11 n1 þ M M @Os

A set of linear equations can be written, in the matrix form, as ¼ R0 a þP0 b ¼ w,

P0 a ¼ 0

ð30Þ

K

where w ¼ fwi gi ¼ 1 indicates the vector of nodal value and matrix: 2 3 R1 ðn1 Þ R2 ðn1 Þ . . . RK ðn1 Þ 6 R ðn Þ R ðn Þ . . . R ðn Þ 7 6 1 2 K 2 7 2 2 6 7 6 7 : : ... : 6 7 ð31Þ R0 ðnÞ ¼ 6 7 : : . . . : 6 7 6 7 6 7 : : ... : 4 5 R1 ðnK Þ R2 ðnK Þ . . . RK ðnK Þ KK

 Z Z L  X ~ 11 x1 Q~ 1 ÞdG þnðlÞ ðM ~ 12 x1 Q~ 2 ÞdG nðlÞ ðM 1 2 Gl

l¼1

ð35Þ

Gl

where L is number of straight line. Consider the radial basis function interpolation in Eq. (34) and relationship between stress and strain in Eq. (11), the first equation Eq. (22a) becomes ( " K L K X X X ðlÞ ðlÞ ðD11 F 1ilj nðlÞ 1 þ SF 2ilj n2 C 1 H 1ilj n1 Þaik k¼1

þ

l¼1 T X

i¼1

#

ðkÞ 2 ðD11 G1tlj n1ðlÞ þSG2tlj n2ðlÞ C 1 I1tlj nðlÞ 1 Þbtk I0 s dkj Os gw1

t¼1

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

643

Support domain L b (x , x )

Local domain 1 x

j

s

2 n x

a (x , x )

Fig. 3. Local integral domain centred at xj with straight boundary lines and its support area.

þ

" K X L K X X

þ

T X t¼1



K X k¼1

ðlÞ ðlÞ ðD12 F 2ilj nðlÞ 1 þ SF 1ilj n2 C 2 H2ilj n2 Þaik

i¼1

k¼1l¼1

#

ðlÞ ðlÞ ðD12 G2tlj nðlÞ 1 þ SG1tlj n2 C 2 I 2tlj n2 Þbtk

(

" L K X X l¼1

J 1ilj aik þ

T X

wðkÞ 2

#

)

K 1tlj btk þ xj1 J 0 s2 dkj Os w3ðkÞ ¼ xj1 q~ 0 Os

t¼1

i¼1

ð36aÞ where dkj is the Kronecker delta function. For the second equilibrium equation in Eq. (22b), one has " K X L K X X ðlÞ ðlÞ ðD12 F 1ilj nðlÞ 2 þ SF 2ilj n1 C 1 H 2ilj n1 Þaik k¼1l¼1

þ

i¼1

T X

ðlÞ ðlÞ ðD12 G1tlj nðlÞ 2 þ SG2tlj n1 C 1 I2tlj n1 Þbtk

t¼1

(

K X

þ

#

k¼1

þ

T X t¼1 K X



" L K X X

l¼1

k¼1

wðkÞ 1

J2ilj aik þ

~ ðkÞ ~ ðkÞ ðD11 fk,1 w 1 þ D12 fk,2 w 2 Þ,

2

I0 s dkj Os

~ 12 ¼ M

wðkÞ 2

K 2tlj btk þ xj2 J 0 s dkj Os

Q~ 1 ¼

wðkÞ 3

K X

~ ðkÞ ~ ðkÞ Sðfk,2 w 1 þ fk,1 w 2 Þ,

K X

~ ðkÞ ~ ðkÞ C 1 ðfk w 1 þ fk,1 w 3 Þ,

k¼1

t¼1

¼ xj2 q~ 0 Os

~ ðkÞ ~ ðkÞ ðD12 fk,1 w 1 þ D22 fk,2 w 2 Þ

k¼1

) 2

K X

k¼1

)

#

T X

K X

k¼1

a

i¼1

l¼1

where ~ 11 ¼ M

#

"

k ¼ 1,2, . . ., M T ð38Þ

ðlÞ ðlÞ ðD22 F 2ilj nðlÞ 2 þ SF 1ilj n1 C 2 H 2ilj n2 Þ ik

i¼1

L K X X

~ ij nj ¼ t~ 0 ðn Þ ði ¼ 1,2Þ and Q~ j nj ¼ t~ 0 ðn Þ on GT , M k i 3 k

~ 22 ¼ M

ðlÞ ðlÞ ðD22 G2tlj nðlÞ 2 þ SG1tlj n1 C 2 I2tlj n2 Þbtk

(

and MO is number of nodes in the domain O. All closed forms of integrals are given in Appendix A. A part from the collocation points in the domain, we need to consider the traction/ displacement boundary conditions. Suppose there are nodes both in the domain and on the boundary,M¼MO þMT þMD, where MT and MD are numbers of nodes on the traction/displacement boundaries. The traction boundary conditions can be introduced directly and are given as

ð36bÞ

Q~ 2 ¼

K X

~ ðkÞ ~ ðkÞ C 2 ðfk w 2 þ fk,2 w 3 Þ

ð39Þ

k¼1

and for the third equilibrium equation Eq. (22c) " # K X L K T X X X ðlÞ ðkÞ C 1 Lilj nðlÞ a þ C M n b 1 tlj 1 tk w1 1 ik k¼1l¼1

þ

k¼1l¼1

(

K X

þ

k¼1

þ

t¼1

i¼1

" K X L K X X

T X

C 2 Lilj nðlÞ 1 ik þ

a

l¼1

# C 2 M tlj nðlÞ 1 btk

~ i ðnk Þ ¼ w ~ 0i , w

wðkÞ 2

t¼1

i¼1

" L K X X

T X

and the shape functions ji are defined in Eq. (34). For the displacement boundary (for example clamped boundary), we also can introduce the displacement equation as

ðC 1 F 1ilj n1ðlÞ þ C 2 F 2ilj nðlÞ 2 Þ ik

a

i¼1

#

ðlÞ ðC 1 G1tlj nðlÞ 1 þ C 2 G2tlj n2 Þbtk

) 2

J 0 s Os dkj wðkÞ 3 ¼ q0 Os

t¼1

ð36cÞ for xj , j ¼ 1,2, . . ., M O , where the integrals Z sl Z sl Z sl @Ri @P t ds, Grtlj ¼ ds, Hrilj ¼ xr Ri ds, F rilj ¼ @xr @xr 0  Z 0sl Z 0sl  @R @Ri ðlÞ Irilj ¼ xr P t ds, J rilj ¼ xr C 1 i nðlÞ n2 ds, 1 þ C2 @x1 @x2 0 0  Z sl  @P t ðlÞ @P t ðlÞ K rtlj ¼ xr C 1 n þ C2 n ds, @x1 1 @x2 2 0 Z sl Z sl Lilj ¼ Ri ds, Mtlj ¼ Pt ds, 0

0

ð37Þ

i ¼ 1,2,3;

k ¼ 1,2, . . ., M D :

ð40Þ

Therefore, there are total 3  (MO þMT þMD) linear algebraic equations, which are used to determine the same number unknowns of displacements either in the domain or on the traction boundary. For the dynamic problems, the total number of samples (Lp þ1) in the transformation space pk, k¼0, 1, 2, y, Lp, are selected in the Laplace transform domain. Transformed variables are evaluated for these specified Laplace parameters by the numerical procedure above. Then, each variable in time domain can be determined by the Laplace inversion technique. Here, the method proposed by Durbin [25] is adopted. The formula of inversion used is written as " # Lp X 2eyt 1 f ðtÞ ¼  f~ ðyÞ þ Reff~ ðy þ 2kpi=TÞe2kpti=T g ð41Þ 2 T k¼0 where f~ ðpk Þ denotes the transformed variable in the Laplace domain; the parameter of the Laplace transform is chosen as

644

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

pffiffiffiffiffiffiffi pk ¼ y þ2kp i/T ð i ¼ 1Þ. There are two free parameters in pk: y and T. The selection of parameters T depends on the observing period in the time domain.

6. Numerical examples In this section, numerical implementations are carried out for the simply supported and clamped plates under static load and/or impact load with the Heaviside time dependence. In order to test the accuracy of the proposed method, the numerical results are compared with the 3D analytical solutions and BEM for both static and dynamic cases.

a21 ¼ a12 ,

2

a22 ¼ D22 b þ Sa2 þ op2 ,

a31 ¼ a13 , a32 ¼ a23 , i rh o ¼ ðhz0 Þ3 þz30 , 3

2

a ¼ mp=a, b ¼ np=b

For the static problems, the Laplace transform parameter s is taken to zero and the symbol ‘‘~’’ for each variable should be removed. To examine the influence by the shape of the local integral domain, firstly we consider a square plate with of thickness h/a¼0.1 and material profile n¼1 in details. The neutral plane position z0 ¼0.01411 m (z0/h¼0.5556) and bending stiffness are D11 ¼ 1:4241  104 Nm, D12 ¼ 2:1362  103 Nm, D22 ¼ 7:1207  103 Nm,

We firstly consider a simply supported orthotropic square plate under uniformly distributed load. A square orthotropic plate with a side-length a ¼0.254 m and the thickness h is subjected to a uniformly distributed static load q0 ¼2.07  106 N/m2 in the vertical direction. The following material parameters are used in numerical analysis: Young’s moduli on the top surface E2t ¼6.895  109 N/m2, E1t ¼2E2t, the ratio of Young’s moduli b ¼E1b/ ¼E1t ¼E2b/E2t ¼ 0.5 (constant) and Poisson’s ratios n21 ¼0.15 and n12 ¼0.3. The used shear moduli are G44t ¼ G55t ¼ G66t ¼E2t/2(1þ n12). The supported domain is selected as a circle of radius dk centred at field point xk[ ¼(xk1,xk2)], which is determined such that the minimum number of nodes in the support domain Lk ¼12. The number of polynomial basis {Pt} is chosen to 5 in radial basis function interpolation. For the rectangular plate, a set of uniformly distributed collocation points (Nx  Ny ¼21  21) are selected in the domain. Free parameter l in (26) is chosen as a/Nx. It has been shown that the selection of this parameter has very little influence on the accuracy of numerical solution. For example, the stable numerical solution can be obtained in the region of parameter l[0.01a/Nx,5a/Nx] in these computation procedures. In general case, the selection of the minimum number Lk is more sensitive than the free parameter l. The numerical results become to unstable when Lk 425. More detailed discussions were presented in the books [17] and [18] by Atluri and Liu, respectively. For the simply supported rectangular plate with composite material, the analy~ thin of moderate thick plate theory (Reisnner’s tical solution w theory) in the transformed domain can be obtained as [1]     1 1 X X sin mapx sin npb y 16q~ 0 ~ thin ¼ ð42Þ w mnBmn p2 m ¼ 1,3,5, ... n ¼ 1,3,5, ...

C 1 ¼ C 2 ¼ 4:2099  107 N=m:

b1 ¼ a12 a23 a13 a22 , 2

b2 ¼ a13 a12 a11 a23 ,

a11 ¼ D11 a2 þ Sb þ op2 ,

a12 ¼ ðD12 þSÞab,

b0 ¼ a11 a22 a212 a13 ¼ C 1 a

S ¼ 2:6155  103 Nm,

To observe the accuracy and convergence of this method, we consider four different shapes for the local integral domain, i.e. triangle, square, octagon and circle (L¼3, 4, 8 and 128) as shown in Fig. 4. The deflection at centre of plate is analysed and the uniformly distributed nodes 441(21  21) is selected. The numerical result is listed in Table 1 and analytical solution is 3 w3D m. In addition, the parameter D=Dmin 3 ða=2,a=2Þ ¼ 4:120  10 varies from 0.1 to 1, D is character parameter of integral domain as shown in Fig. 4. To examine the accuracy of thin plate theory and Reissner moderate thick plate theory, the analytical solutions of deflection on the middle plane by 3D elasticity analysis [27] are presented in the tables for comparison. The relative errors percentage for deflection of plate defined by

D¼

9w3 w3D 3 9  100 w3D 3

ð43Þ

are shown in Table 1. Table 1 Deflection at centre of plate w3(a/2, a/2) with different local integral domain shape. D/Dmin

Triangle

Square

Octagon

Circle

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

4.1105 4.1104 4.1101 4.1100 4.1098 4.1094 4.1091 4.1087 4.1084 4.1082

4.1105 4.1103 4.1101 4.1098 4.1094 4.1090 4.1086 4.1081 4.1076 4.1071

4.1104 4.1103 4.1100 4.1096 4.1092 4.1086 4.1081 4.1077 4.1073 4.1070

4.1104 4.1103 4.1099 4.1095 4.1090 4.1085 4.1080 4.1076 4.1072 4.1068

Error (%)

0.26

0.27

0.27

0.27

D

L=3

L=4

2

a33 ¼ C 1 a þ C 2 b þ rhp2

6.1. A simply supported rectangular plate with orthotropic materials under static load

where

a23 ¼ C 2 b

L=8

Fig. 4. Different shapes of local integral domain and character parameter D.

L=128

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

645

Deflection of the plate on the central line y¼0.5a is plotted in Fig. 5 with different numbers of collocation points, i.e. 5  5, 11  11 and 21  21, respectively. It is difficult to see the difference between 3D analytical solutions [27] and LIEM results when Nx and Ny are larger than 11. Finally, the variation of the direct stresses sx and sy along thickness at central point of the plate are presented in Fig. 6 when the collocation numbers is taken to Nx  Ny ¼21  21. It is clear that the neutral axes (sx ¼0 and sy ¼0) on two sections shift up with the same distance when the parameter of material profile n is not zero (Et 4 Eb). The agreement with 3D analytical solution is shown to be excellent.

The numerical results given by LIEM, analytical solutions in (42) for moderate thick plate theory and 3D analytical solutions are given in Tables 2–4 for different parameters of material profile n and ratios of plate height and width b/a. It is clear that the satisfactory accurate solution can be obtained using moderate thick plate theory as long as h/ao0.3 for different ratio of a/b. In addition, the relative errors between analytical solution (42) and solutions given by LIEM are less than 1% when h/a 40.05. However, the relative errors increase significantly to use LIEM when h/ao0.025. In this case, accuracy can be improved by increasing the number of collocation point.

Table 2 Relative errors compared with 3D analytical solution when b/a¼ 1. h/a

0.010 0.025 0.050 0.075 0.100 0.125 0.150 0.200 0.250 0.300 0.400 0.500

n¼1

n¼ 3

3D [27]

Analytical

%

LIEM

%

3D [27]

Analytical

%

LIEM

%

3.888Eþ 00 2.497E  01 3.158E  02 9.526E  03 4.120E  03 2.175E  03 1.306E  03 6.004E 04 3.397E  04 2.191E  04 1.158E  04 7.396E  05

3.884Eþ 00 2.494E  01 3.154E  02 9.525E  03 4.125E  03 2.182E  03 1.312E  03 6.065E  04 3.452E  04 2.243E  04 1.208E  04 7.896E  05

0.1 0.1 0.1 0.0 0.1 0.3 0.5 1.0 1.6 2.4 4.3 6.8

– – 3.159E  02 9.502E  03 4.109E  03 2.173E  03 1.306E  03 6.039E  04 3.439E  04 2.235E  04 1.204E  04 7.878E  05

– – 0.1 0.3 0.2 0.1 0.1 0.6 1.2 2.0 4.0 6.5

4.324E þ00 2.776E  01 3.494E  02 1.057E-02 4.586E  03 2.432E  03 1.466E  03 6.811E  04 3.894E  04 2.537E  04 1.365E  04 8.824E  05

4.329E þ00 2.780E  01 3.519E  02 1.064E  02 4.618E  03 2.448E  03 1.476E  03 6.862E  04 3.930E  04 2.568E  04 1.397E  04 9.205E  05

0.1 0.2 0.7 0.7 0.7 0.7 0.7 0.8 0.9 1.2 2.4 4.3

– – 3.523E  02 1.062E  02 4.600E  03 2.438E  03 1.470E  03 6.833E  04 3.915E  04 2.559E  04 1.393E  04 9.185E  05

– – 0.8 0.5 0.3 0.3 0.3 0.3 0.5 0.9 2.1 4.1

Table 3 Relative errors compared with 3D analytical solution when b/a¼ 0.5. h/a

0.010 0.025 0.050 0.075 0.100 0.125 0.150 0.200 0.250 0.300 0.400 0.500

n ¼1

n¼3

3D [27]

Analytical

%

LIEM

%

3D [27]

Analytical

%

LIEM

%

7.402E  01 4.767E  02 6.087E  03 1.866E  03 8.241E  04 4.461E  04 2.751E  04 1.340E  04 8.013E  05 5.424E  05 3.069E  05 2.014E  05

7.392E  01 4.761E  02 6.087E  03 1.871E  03 8.288E  04 4.504E  04 2.791E  04 1.376E  04 8.348E  05 5.753E  05 3.416E  05 2.400E  05

0.1 0.1 0.0 0.2 0.6 1.0 1.5 2.7 4.2 6.1 11.3 19.2

– – 6.063E  03 1.854E  03 8.206E  04 4.459E  04 2.764E  04 1.364E  04 8.285E  05 5.716E  05 3.401E  05 2.392E  05

– – 0.4 0.6 0.4 0.0 0.5 1.8 3.4 5.4 10.8 18.8

8.225E  01 5.301E  02 6.750E  03 2.079E  03 9.238E  04 5.035E  04 3.128E  04 1.547E  04 9.375E  05 6.418E  05 3.688E  05 2.439E  05

8.239E  01 5.310E  02 6.800E  03 2.095E  03 9.315E  04 5.082E  04 3.163E  04 1.572E104 9.613E  05 6.669E  05 4.001E  05 2.830E  05

0.2 0.2 0.7 0.8 0.8 0.9 1.1 1.6 2.5 3.9 8.5 16.0

– – 6.768E  03 2.077E  03 9.221E  04 5.031E  04 3.132E  04 1.559E  04 9.544E  05 6.629E  05 3.984E  05 2.821E  05

– – 0.3 0.1 0.2 0.1 0.1 0.8 1.8 3.3 8.0 15.6

Table 4 Relative errors compared with 3D analytical solution when b/a¼ 2. h/a

0.010 0.025 0.050 0.075 0.100 0.125 0.150 0.200 0.250 0.300 0.400 0.500

n¼1

n¼ 3

3D[27]

Analytical

%

LIEM

%

3D[27]

Analytical

%

LIEM

%

7.290Eþ 00 4.671E  01 5.895E  02 1.773E  02 7.640E  03 4.015E  03 2.397E  03 1.089E  03 6.089E  04 3.881E  04 2.014E  04 1.271E  04

7.268Eþ 00 4.664E  01 5.887E  02 1.772E  02 7.640E  03 4.019E  03 2.402E  03 1.095E  03 6.146E  04 3.936E  04 2.067E  04 1.324E  04

0.3 0.1 0.1 0.1 0.0 0.1 0.2 0.5 0.9 1.4 2.6 4.2

– – 5.923E  02 1.766E  02 7.587E  03 3.986E  03 2.381E  03 1.085E  03 6.091E  04 3.903E  04 2.052E  04 1.317E  04

– – 0.5 0.4 0.7 0.7 0.7 0.4 0.0 0.6 1.9 3.6

8.102Eþ 00 5.191E  01 6.522E  02 1.966E  02 8.493E  03 4.479E  03 2.684E  03 1.231E  03 6.944E  04 4.468E  04 2.357E  04 1.507E  04

8.101Eþ 00 5.200E  01 6.567E  02 1.979E  02 8.547E  03 4.505E  03 2.699E  03 1.237E  03 6.977E  04 4.493E  04 2.383E  04 1.539E  04

0.0 0.2 0.7 0.7 0.6 0.6 0.5 0.5 0.5 0.6 1.1 2.1

– – 6.600E  02 1.972E  02 8.485E  03 4.467E  03 2.674E  03 1.225E  03 6.915E  04 4.456E  04 2.366E  04 1.531E  04

– – 1.2 0.3 0.1 0.3 0.4 0.5 0.4 0.3 0.4 1.6

646

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

0.0045 0.0040 0.0035

w3 (m)

0.0030 0.0025 5X5 11X11 21X21 Analytical

0.0020 0.0015 0.0010 0.0005 0.0000 0.00

0.05

0.10

0.15

0.20

0.25

x (m) Fig. 5. Variation of deflection vs the x coordinate for orthotropic functionally graded plate (n¼ 1, h/a ¼0.1).

0.030

0.025

σy

σx

z’ (m)

0.020

0.015

0.010

LIEM 3D analytical [27]

0.005

0.000 -8.0E+07 -6.0E+07 -4.0E+07 -2.0E+07 0.0E+00

2.0E+07

4.0E+07

6.0E+07

8.0E+07

1.0E+08

1.2E+08

2

σ (N/m ) 0

Fig. 6. Variation of stresses along z coordinate for orthotropic functionally graded plate (n¼ 1,h/a¼ 0.1) with 21  21 collocation points.

6.2. A square plate with functionally graded material with different types of supported edges Considering dynamic problem, a simply supported square plate with the side-length a¼b¼0.254 m, thickness h/a¼0.1 and graded orthotropic material is studied firstly. The plate is subjected to a uniformly distributed load q(x,y)¼q0H(t), where H(t) is Heaviside step function and uniform load q0 ¼2.07  106 N/m2. The orthotropic material parameters are the same as that in static case above. The exact solution by 3D analysis [27] orthotropic homogenous material (n¼0) under static uniformly distributed load 3 q0 gives w3D mm on the middle plane at the 3 ¼ 2:9715  10 central of plate and stresses s3D,b ¼ 8:4824  107 N=mm2 and x 7 2 s3D,t ¼ 8:4971  10 N=mm on the bottom and top surfaces, x respectively. Two parameters used in the inversion of the Laplace transformation are selected as: y ¼5/t0 and T¼60t0, where t0 ¼a/c0 as defined inffi [26] (c0 is the velocity of longitudinal elastic wave: pffiffiffiffiffiffiffiffiffiffiffiffiffi c0 ¼ c33t =r), c33t ¼E2t(1 n12n21)/D, D ¼ 1n222 2n12 n21 ð1 þ n22 Þ.

The total number of collocation points is taken to L¼21  21 and the total number of samples in the Laplace domain K¼ 25. The normalised deflection w3 ðtÞ=w3D and stresses sbx ðtÞ=s3D,b , x 3 stx ðtÞ=s3D,t at the centre of the bottom and top surfaces of the x plate are plotted, respectively, in Figs. 7 and 8. Two cases are considered: (a) n¼0 for constant orthotropic homogenous material and (b) n¼ 2 for quadratic functionally graded orthotropic material. In case (a), the stresses on the centre bottom and top surfaces are the same by Reissner’s theory. However in case (b), the absolute value of the stress on the bottom surface is much larger than that on the top surface due to the much stiffer material in the bottom of the plate. 6.3. A clamped circular plate under static load In order to show the validity of proposed method, a circular thick plate with clamped edge boundary is examined. To compare with the boundary element method [4], the isotropic homogenous

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

647

3.5 3.0

n=2

2.5

w3(t) /w33D

2.0 1.5

n=0

1.0 0.5 0.0

0

5

10

15

20

-0.5

25

30

35

40

45

50

c0t/a

Fig. 7. Variation of deflection vs normalised time ct/a for orthotropic uniformly and functionally graded plate (n¼ 2, h/a¼ 0.1). Solid lines indicate the results by 3D analysis [27] and dash lines by LIEM.

3.0 n = 2, σ / σ

2.5

n = 0, σ /

2.0

 (t) /

1.5

1.0 n = 2, /

0.5

0.0

0

5

10

15

20

-0.5

25

30

35

40

45

50

c t /a

Fig. 8. Variation of deflection vs the x coordinate for orthotropic uniformly and functionally graded plate (n¼ 2, h/a¼ 0.1). Solid lines indicate the results by 3D analysis [27] and dash lines by LIEM.

thick plate is considered. The geometric and material properties are as follows: R¼1 m, the ratio of Young’s modulus and uniform pressure E/q0 ¼210,000, E1 ¼E2, n¼0, Poisson’s ratios n21 ¼ n12 ¼0.3. A circle of radius dn for support domain is selected such that the minimum number of nodes in the support domain LZ ¼12. The distribution of collocation points in the domain is shown in Fig. 9 with total number of M¼390. Fig. 10 shows the results of deflection w3(x,0) and rotation w1(x,0) when thickness h ¼0.2 m and h ¼0.3 m, respectively. The numerical results given by BEM are presented for comparison. It can be seen that the agreement with BEM is excellent. Finally a clamped circular orthotropic functionally graded plate with the same material parameters as in Section 6.1 above is considered. The radius of the circular plate R¼ 0.127 m (a/2), the thickness h¼ 0.2R (h/a¼0.1) and the uniformly distributed static load q0 ¼ 2.07  106 N/m2. The collocation point distribution is the same as in above analysis. The variations of deflection and rotation along the centre line are shown in Figs. 11 and 12, respectively.

Fig. 9. Distribution of the collocation point for a circular plate.

648

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

2.0E-04 w (m), h/R=0.2

1.5E-04 w (m), h/R=0.2

1.0E-04

w (m), h/R=0.3 5.0E-05

-1.0

-0.8

-0.6

-0.4

0.0E+00 -0.2 0.0

0.2

0.4

-5.0E-05 w (m), h/R =0.3

0.6

0.8

1.0

x (m)

-1.0E-04

LIEM BEM[4]

-1.5E-04 -2.0E-04

Fig. 10. Variations of deflection and rotation along the centre line of the clamped supported circular plate.

1.4E-03 n=0

1.2E-03

n=1 1.0E-03

n=2

6.0E-04

w3 (m)

8.0E-04

4.0E-04 2.0E-04

-1.0

-0.8

-0.6

-0.4

0.0E+00 -0.2 0.0

0.2

0.4

0.6

0.8

1.0

x (m) Fig. 11. Variations of deflection along the centre line of the clamped supported circular plate for different parameter of material graded profile n.

0.015

0.010

w1

0.005

x (m) -1.0

-0.8

-0.6

-0.4

-0.2

0.000

-0.005

0.0

0.2

0.4

0.6

0.8

1.0

n=0 n=1

-0.010

n=2

-0.015 Fig. 12. Variations of rotation along the centre line of the clamped supported circular plate for different parameter of material graded profile n.

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

The maximum value of deflection at the centre of plate is w03 ð0,0Þ ¼ 0:8045  103 m when n¼0.

649

Then the closed form for Jrilj can be obtained: J1ilj ¼ C 1 ðxa1 ðxa1 xi1 ÞT 0 þð2xa1 xi1 ÞT 1 cos bl þ T 2 cos2 bl Þsin bl C 2 ðxa1 ðxa2 xi2 ÞT 0 þ ðxa1 sin bl þðxa2 xi2 Þ cosbl ÞT 1 þ T 2 sin bl cos bl Þcos bl

7. Conclusions

J2ilj ¼ C 1 ðxa2 ðxa1 xi1 ÞT 0 þðxa2 cos bl þðxa1 xi1 Þ sin bl ÞT 1 The meshless local integral equation method (LIEM) is presented to solve bending problems for the moderate thick plates with functionally graded material in this paper. Firstly, the governing equations for functionally graded material plate were derived with Reissner plate theory. The coupled effect between in-plane and out-plane deformations disappears in the Laplace space. For LIEM, all domain integrals in the weak forms of governing equations are transformed to the boundary integrals. In this paper, the analytical formulations for the meshless local integral equation method were presented. Hence there are not any domain and boundary integrals to be determined numerically in the computation process and consequently the computation time will be reduced. Many observations were carried out for several examples, such as the influence of the parameter of integral domain (D=Dmin ) and the number of distributed nodes (M) in the domain. Comparisons were made with the analytical solution for 3D investigation and excellent accuracy was achieved. By observing several numerical examples, it can be conclude with the following: (1) meshless method LIEM can be fulfilled without any integrals to be evaluated numerically; (2) computation time is reduced; (3) numerical solutions are stable and convergent for suitable selection of free parameters. We can also conclude that LIEM is one of the simplest methods in the group of mesh free methods and it should be extended to solve more complicated problems including elastoplasticity, nonlinear and contact problems.

Appendix A Consider n1ðlÞ ¼ sin bl , n2ðlÞ ¼ cos bl , we have solutions in closed form for Frilj: F 1ilj ¼ ðxla1 xi1 ÞT 0 þ cos bl T 1 F 2ilj ¼ ðxla2 xi2 ÞT 0 þ sin bl T 1

ðA1Þ

where T 0 ¼ ln

T1 ¼ a¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b þ s þ a2 þ 2bs þ s2 bþa

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bþ s þ a2 þ 2bs þ s2 a2 þ 2bs þ s2 ab ln bþ a

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l2 þ ðxa1 xi1 Þ2 þ ðxa2 xi2 Þ2 ,

ðA2Þ

b ¼ ðxa1 xi1 Þcos bl þ ðxa2 xi2 Þsin bl ,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s ¼ ðxb1 xa1 Þ2 þðxb2 xa2 Þ2 Similarly we have solutions in closed form for Hrilj: H1ilj ¼ xa1 S0 þcos bl S1 ,

H2ilj ¼ xa2 S0 þ sin bl S1

ðA3Þ

where

S0 ¼

S1 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðb þ sÞ a2 þ2bs þ s2 ab a2 b b þ sþ a2 þ 2bs þs2 þ ln 2 bþa 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 bðb þsÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a3 ab ða2 þ2bs þ s2 Þ3  a2 þ 2bs þ s2  þ 2 3 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 bðb a2 Þ b þs þ a2 þ2bs þs2 ln þ 2 b þa 1 3

ðA4Þ

þ T 2 sin bl cos bl Þsin bl C 2 ðxa2 ðxa2 xi2 ÞT 0 þ ð2xa2 xi2 ÞT 1 sin bl þ T 2 sin2 bl Þcos bl

ðA5Þ

where T2 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðs3bÞ a2 þ 2bs þ s2 þ 3ab b þ sþ a2 þ 2bs þs2 2 þð3b a2 Þln 2 bþa ðA6Þ

Finally for Lilj, we have pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðb þ sÞ a2 þ 2bs þ s2 ab ða2 b Þ b þ s þ a2 þ 2bs þ s2 ln þ Lilj ¼ 2 2 bþa ðA7Þ It is very easy to find analytical solutions for Grtlj, Irtlj, Krtlj and Mrtlj. References [1] Aliabadi MH. The Boundary Element Method, Application in Solid and Structures. Chichester: Wiley; 2002. [2] Van der Ween F. Application of the boundary integral equation method to Reissner’s plate model. Int J Numer Methods Eng 1982;18:1–10. [3] Antes H. Static and dynamic analysis of Reissner–Mindlin plates. In: Beskos DE, editor. Boundary Element Analysis of Plates and Shells. Berlin: Springer; 1991. [4] Wen PH, Aliabadi MH. Boundary element frequency domain formulation for dynamic analysis of Mindlin plates. Int J Numer Methods Eng 2006;67: 1617–40. [5] Wen PH. Point intensity method of solving circular plate resting on elastical subgrade. Eng Mech 1987;4(2):18–26. [6] Wen PH. The numerical method for complex restrained problem of rectangular plate on elastic base. J Cen South Inst Min Metall 1988;19(3): 346–8. [7] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10: 307–18. [8] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [9] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [10] Wen PH, Aliabadi MH, Lin YW. Meshless method for crack analysis in functionally graded materials with enriched radial base functions. CMESComput Modeling Eng Sci 2008;30:133–47. [11] AJM Ferreira, Batra RC, Roque CMC, Qian LF, Martins PALS. Static analysis of functionally graded plates using third-order shear deformation theory and a meshless method. Compos Struct 2005;69:449–57. [12] Gilhooley DF, Batra RC, Xiao JR, McCarthy MA, Gillespie Jr JW. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Compos Struct 2007;80:539–52. [13] 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 Methods Eng 2006;67:1830–50. [14] Li LY, Wen PH, Aliabadi MH. Meshfree modeling and homogenization of 3D orthogonal woven composites. Compos Sci Technol 2011;71:1777–88. [15] Wen PH, Aliabadi MH. An improved meshless collocation method for elastostatic and elastodynamic problems. J Commun Numer Methods Eng 2008;24(8):635–51. [16] Wen PH, Hon YC. Geometrically nonlinear analysis of Reissner–Mindlin plate by meshless computation. CMES: Comput Modeling Eng Sci 2007;21(3): 177–91. [17] Atluri SN. The Meshless Method (MLPG) for Domain and BIE Discretizations. Forsyth, GA, USA: Tech Science Press; 2004. [18] Liu GR. Meshfree Methods Moving Beyond the Finite Element Method. NY: CRC Press; 2010. [19] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;176:1905–15. [20] Franke R. Scattered data interpolation: test of some methods. Math Comput 1982;38:181–200. [21] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics – I. Comput Math Appl 1990;19(8/9):127–45.

650

P.H. Wen, M.H. Aliabadi / Engineering Analysis with Boundary Elements 36 (2012) 639–650

[22] Golberg MA, Chen CS. The theory of radial basis functions applied to the BEM for inhomogeneous partial differential equations. Bound Elem Commun 1994;5:57–61. [23] Suresh S, Mortensen A. Fundamentals of Functionally Graded Materials. London: IOM Communications; 1998. [24] Miyamoto Y, Kaysser WA, Rabin BH, et al. Functionally Graded Materials: Design, Processing and Applications. MA: Kluwer; 1992.

[25] Durbin F. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method. Comput J 1975;17:371–6. [26] Timoshenko SP, Woinowsky-Krieger S. Theory of Plates and Shells. Singapore: McGraw-Hill; 1959. [27] Wen PH, Sladek J, Sladek V. Three-dimensional analysis of functionally graded plates. Int J Numer Methods Eng 2011;87:923–42.