- Email: [email protected]
Three-dimensional static analysis of thick functionally graded plates by using meshless local Petrov–Galerkin (MLPG) method
ARTICLE IN PRESS Engineering Analysis with Boundary Elements 34 (2010) 564–573
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Three-dimensional static analysis of thick functionally graded plates by using meshless local Petrov–Galerkin (MLPG) method R. Vaghefi, G.H. Baradaran , H. Koohkan Department of Mechanical Engineering, Shahid Bahonar University of Kerman, Kerman, Iran
a r t i c l e in f o
a b s t r a c t
Article history: Received 30 July 2009 Accepted 20 January 2010 Available online 20 February 2010
In this paper, a version of meshless local Petrov–Galerkin (MLPG) method is developed to obtain threedimensional (3D) static solutions for thick functionally graded (FG) plates. The Young’s modulus is considered to be graded through the thickness of plates by an exponential function while the Poisson’s ratio is assumed to be constant. The local symmetric weak formulation is derived using the 3D equilibrium equations of elasticity. Moreover, the field variables are approximated using the 3D moving least squares (MLS) approximation. Brick-shaped domains are considered as the local sub-domains and support domains. In this way, the integrations in the weak form and approximation of the solution variables are done more easily and accurately. The proposed approach to construct the shape and the test functions make it possible to introduce more nodes in the direction of material variation. Consequently, more precise solutions can be obtained easily and efficiently. Several numerical examples containing the stress and deformation analysis of thick FG plates with various boundary conditions under different loading conditions are presented. The obtained results have been compared with the available analytical and numerical solutions in the literature and an excellent consensus is seen. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Three-dimensional meshless local Petrov–Galerkin method Functionally graded material Three-dimensional theory of elasticity Thick plate Moving least squares approximation
1. Introduction Functionally graded materials (FGMs) are advanced composites in which the material properties vary smoothly and continuously by a predetermined function. Generally, those are composed of ceramics and metals so that the material properties such as the Young’s modulus vary from the metallic side to the ceramic one. Initially, they were designed as thermal barrier coatings in space applications. However, recently FGMs have gained lots of applications in nuclear reactors [1], dental and medical implants [2], piezoelectric and thermo electric devices [3–5] and fire retardant doors [6]. Investigation of the mechanical behavior of new materials such as FGMs is an attractive research area in mechanics. Many researchers have studied the mechanical behavior of FGMs having various geometries and loading conditions. Among the geometries, which are considered to analyze, the plates are the most significant due to numerous applications in engineering structures. It should be mentioned that the plate is a three-dimensional (3D) structure in which one dimension is relatively much smaller than the other dimensions. In 2D theories such as the classical plate theory (CPT), in the first order shear deformation plate
Corresponding author. Tel.: + 98 341 2111763; fax: + 98 341 2120964.
E-mail address: [email protected] (G.H. Baradaran). 0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.01.005
theory (FSDT) and the higher order shear deformation plate theory (HSDT), various assumptions are made to obtain a 2D formulation for plates. Obviously, it is easier to find a solution for the plate problems in 2D formulation. However, some errors occur in solutions due to the assumptions considered in 2D theories. With the increase in the thickness of plates, the errors amplify. If no simplifying assumptions are considered, the 3D elasticity equations must be used for elastic analysis of plates. The 3D solution does not involve any limitation of 2D solutions. It is obvious that if 3D solutions are attainable, it will be more accurate than the solutions obtained by the mentioned 2D theories. Analytical and numerical methods have been used for 3D analyses of plates [7–13]. Analytical methods are applicable for some simple cases and generally, the analytical solutions cannot be found for plates with complicated geometries and boundary conditions. Consequently, numerical methods such as the finite element method (FEM) and differential quadrature (DQ) method are the alternatives that have been employed for 3D analysis of plates. Recently, meshless methods, which do not need burdensome effort of mesh generation have gained lots of attention. A variety of these methods have been developed, such as the element-free Galerkin (EFG) method [14], the reproducing kernel particle method (RKPM) [15], hp-clouds [16], and the partition of unity method (PUM) [17], which have been successfully applied for the solutions of some engineering problems. However, most of
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
these methods are not truly meshless, because, they use background cells for the evaluation of integrals in the weak formulation of problems. The meshless local Petrov–Galerkin (MLPG) method proposed by Atluri et al. [18–21] considers local weak formulation instead of the global weak formulation, and does not require a background mesh for evaluation of the integrals. In the MLPG method, the integrals of the weak form are evaluated over local subdomains that partly cover each other. The trial and test functions are chosen from totally different functional spaces. Furthermore, the physical size of the test and trial domains is not necessary to be the same, which makes the MLPG a very flexible method. Based on the concept of the MLPG, six different methods have been introduced, which are labeled as MLPG1–MLPG6 [21]. These six methods differ due to the type of test function considered in the weak formulation. Among them, the MLPG1 and MLPG5 seem to be the most promising formulations; so they are used in the present study. In the MLPG1 a namely fourth-order spline function is considered as test function, while the Heaviside step function is employed as test function in MLPG5. It is noticeable that the MLPG5 does not involve any domain integration or singular integrals. The MLPG methods have been employed in a wide range of applications, for example elasto-statics [22], elastodynamics [23], fluid mechanics [24], convection–diffusion problems [25], thermoelasticity [26], beam problems [27,28], plate problems [29–31], FGM problems [32,33], fracture mechanics [34–36], and strain gradient theory [37]. So far, the applications of the MLPG method are mostly limited to 2D problems. The major reason for the 3D analysis is considered less, related to difficulties in evaluating the integrals over the local sub-domains, especially when a sub-domain intersects the global boundary of problem. However, significant efforts have been carried out to develop the MLPG method to solve 3D problems. Li et al. [38] applied a combination of MLPG2 and MLPG5 for the solution of two classical 3D problems, viz., the Boussinesq problem and the Eshelby’s inclusion problem. They used MLPG5 for nodes inside the domain and MLPG2 for nodes on the boundary of the problem domain. Han and Atluri [39] used the MLPG method to solve 3D elastic fracture problems. Also, Han and Atluri [40] developed three kinds of the MLPG methods using different test functions for analysis of thick beams and spheres and examined the performance of the methods. They used Delaunay triangulation algorithm and introduced a suitable mapping procedure for evaluating the integrals over the spherical sub-domains. The integrals on the volume of the sub-domains were calculated by dividing the domains into cone shapes. The MLPG domain discretization method has been employed to 3D elasto-dynamic problems of impact and fragmentation by Han and Atluri [41]. Vaghefi et al. [42] developed two different 3D MLPG procedures including MLPG1 and MLPG5 for the elasto-static analysis of thick rectangular plates with various boundary conditions. Brick-shaped domains are considered as sub-domains and support domains. The purpose of the present paper is to develop a complete 3D MLPG procedure for the elasto-static analysis of thick FG plates with various boundary conditions. The Young’s modulus of the plate is assumed to vary exponentially through the thickness, and the Poisson’s ratio is assumed to be constant. The local symmetric weak form (LSWF) is used to formulate the problem and 3D MLS approximation is used to approximate the field variables. To impose the essential boundary conditions, the penalty method is adopted. Since the global domain of the problem is cubic, brickshaped domains are considered as sub-domains and support domains for evaluating the integrals of the weak form and approximating the solution variables, respectively. Moreover, the approach used for the construction of the shape and the test functions makes it possible to add nodes in any direction of the plate (through length, width and thickness) depending upon the
565
required accuracy. Therefore, sufficient number of nodes has been added in the thickness direction of the plate to represent the normal and shear stress variations through the thickness with good accuracy. Consequently, more precise solutions can be obtained easily and efficiently. In order to illustrate the accuracy of the proposed approach, several numerical examples are presented and the results are compared with the known solutions in the literature.
2. MLPG formulation for 3D elasticity The 3D elasto-static equilibrium equations in a domain of the volume O, which is bounded by the surface G, are given by
sij;j þ bi ¼ 0; in O
ð1Þ
where sij is the stress tensor and bi is the body force vector. The indices i, j which take the values of 1, 2, 3 refer to the Cartesian coordinates x, y, z, respectively. The boundary conditions are assumed to be ui ¼ u i on Gu ;
ð2aÞ
sij nj ¼ t i on Gt
ð2bÞ
where ui and ti represent the displacement components and the surface traction components, respectively. Gu is the boundary with the prescribed displacement u i while Gt is the boundary with the prescribed traction t i . nj denotes the unit outward vector normal to the boundary G. The weak formulation is constructed over the local subdomains, which are located inside the global domain O. The local sub-domains may overlap with each other and must cover the whole global domain. Various shapes with different sizes can be chosen as sub-domains, but appropriate ones should be considered to obtain precise outcomes. Since the global domain of the rectangular plate is hexahedral, brick-shaped local domains are considered as sub-domain and support domain (see Fig. 1). Selecting the brick-shaped local domain makes the mapping procedure easy and no special treatment is needed when the local sub-domain intersects the global boundary. 2.1. Local symmetric weak form for 3D elastic body The generalized local weak form of the equilibrium equations over a local sub-domain around node I is written as follows: Z Z ðsij;j þ bi ÞnI dOs a ðui u i ÞnI dG ¼ 0 ð3Þ OIs
GIsu
where ui, nI are the trial and the test functions, respectively and can be chosen from different functional spaces. In Eq. (3), the term
local boundary ∂Ωs=Γsi for an internal node
sub-domain Ωs of node I
x
Γu
Γ
Γt
I Γsi
h
z
Γs ∂Ωs=ΓsiUΓs
Ω
y
Fig. 1. Brick-shaped sub-domains in the global domain of a rectangular plate.
ARTICLE IN PRESS 566
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
GIsu illustrates a part of the boundary @OIs over which the essential boundary conditions are prescribed and a denotes a penalty parameter, which is introduced in order to impose the essential boundary conditions (a b1). Using the relation sij,jnI = (sijnI),i sijnI,j and the divergence theorem, Eq. (3) leads to Z Z Z sij nj nI dG ðsij nI;j bi nI Þ dOs a ðui u i ÞnI dG ¼ 0 ð4Þ
where fI(x) is called the MLS approximation shape function and is defined as
It should be noted that the boundary @OIs of the local subdomain is divided into three parts, @OIs ¼ GIsi [ GIsu [ GIst , where GIsi is a part of the local boundary located inside the global domain and has no contact with the global boundary of the problem. Also, GIst is a part of the local boundary that coincides with the global traction boundary and GIsu is a part of the local boundary that coincides with the global displacement boundary. Writing Eq. (4) in terms of GIsi , GIst and GIsu and applying the natural boundary condition ti ¼ sij nj ¼ t i on GIst , the following local symmetric weak form (LSWF) for the linear 3D elasticity problem can be obtained Z Z Z Z sij nI;j dO ti nI dG ti nI dG þ a ui nI dG
uh;i ¼
@OIs
OIs
OIs
GIsi
Z
¼ GIst
GIsu
GIsu
Z
t j nI dG þ a
GIsu
u i nI dG þ
GIsu
Z OIs
bi nI dO
ð5Þ
The test function nI(x) is chosen such that it is positive inside local domain Os and vanishes outside Os.
2.2. 3D approximation using moving least square (MLS) approximation The MLS approximation is considered such that the unknown trial approximant uh(x) of the function u(x) is defined by uh ðxÞ ¼ pT ðxÞaðxÞ
8x A Ox
ð6Þ
T
where p ¼ ½p1 ðxÞ; p2 ðxÞ; . . . ; pm ðxÞ is a complete monomial basis of order m and Ox is the domain of definition of the MLS approximation for the trial function at point x. The quadratic basis vector pT for 3D problems are defined as pT ðxÞ ¼ ½1; x; y; z; x2 ; y2 ; z2 ; xy; yz; zx
quadratic basis
ð7Þ
The coefficient vector a(x) is determined by minimizing a weighted discrete L2-norm, which is defined as JðaðxÞÞ ¼
N X
I
^ T GðxÞ½PaðxÞu ^ gI ðxÞ½pðxI ÞaðxÞu^ 2 ¼ ½PaðxÞu
ð8Þ
fI ðxÞ ¼
m X
pj ðxÞ½A1 ðxÞBðxÞjI
ð11Þ
j¼1
The partial derivatives of the trial function are introduced as follows N X
fI;i ðxÞu^ I
ð12Þ
I¼1 I
where (.),i denotes q(.)/qxi and f;i ðxÞ are derivatives of the MLS shape function and can be obtained as
fI;i ðxÞ ¼
m X
½pj;i ðA1 BÞjI þ pj ðA1 B;i þ A1 ;i BÞjI
ð13Þ
j¼1 1 where A1 Þ;i represents the derivative of the inverse of the ;i ¼ ðA matrix A with respect to xI. A fourth-order spline weight function for a one-dimensional domain is defined as follows: 8 2 3 4 > < 16 dI þ8 dI 3 dI 0 rdI r rI dI rI rI rI ¼ ð14Þ gI ðxÞ ¼ f > rI :0 d I Z rI
In Eq. (14), dI =jx xIj while rI is the size of the support for the I weight function gI defined as rI ¼ ai :d . ai is the dimensionless size I of the support domain and d is the average nodal spacing between two neighboring nodes in the vicinity of node I. Using the brick-shaped support domain, the weight function for the 3D problem can be obtained by a simple extension of the onedimensional function of (14) as follows: ! I dIy d dI ð15Þ gI ðxÞ ¼ f Ix f I f Iz ¼ gx gy gz rx ry rz where the functions gx, gy and gz are obtained by replacing dI/rI with dIx =rxI , dIy =ryI and dIz =rzI in Eq. (14), respectively. The parameters dIx ; rxI ; dIy ; ryI ; dIz and rzI are defined as I
I
I
dIx ¼ jxxI j; dIy ¼ jyyI j; dIz ¼ jzzI j; rxI ¼ ai d x ; ryI ¼ ai :d y ; rzI ¼ ai d z ð16Þ I
I
I
In Eq. (16), d x , d y and d z are the average nodal spacings between two neighboring nodes in the vicinity of the node I in the x-, y- and z-direction, respectively. Then, the dimensions of the support domain used for constructing the MLS shape functions I I I become ð2ai d x Þ ð2ai d y Þ ð2ai d z Þ. 2.3. Test function
I¼1
where gI(x) is the weight function related to node I and gI(x) 40 for all points placed in the support of the weight function. In Eq. I (8), xI is the position of the Ith node, u^ ðI ¼ 1; 2; . . . ; N A Ox Þ is the fictitious nodal value, which may not be equal to uh(xI) and P, G and u^ have the following form 2
2
3 pT ðx1 Þ 6 T 7 6 p ðx2 Þ 7 7 P¼6 6 ^ 7 4 5 pT ðxN Þ
6 6 ; G¼6 6 4
Nm
g1 ðxÞ 0 ^ 0
0 &
& 0
0 ^
3
7 7 7 0 7 5 gN ðxÞ
; u^ ¼ ½u^ 1 ; u^ 2 ; . . . ; u^ N 1N
NN
ð9Þ The coefficient vector a(x) is determined and substituted into I Eq. (6) to obtain the approximation function uh ðu^ Þ as follows: uh ðxÞ ¼ Uu^ ¼
N X I¼1
fI ðxÞu^ I
ð10Þ
The test function is considered to be either a fourth-order spline function or a Heaviside step function and the corresponding MLPG formulations are called MLPG1 and MLPG5, respectively. Brick-shaped sub-domain is also chosen for the support of the test function. In a same manner discussed in the previous section, the dimensions of the sub-domain for node I can be I I I defined as ð2as d x Þ ð2as d y Þ ð2as d z Þ, Where as is a constant, greater than 0 and less than 1, chosen as 0.75 in the present study. Considering the MLPG1 method, the weight function of Eq. (14) has been used to define the test function. The test function for the brick-shaped sub-domain is constructed following the same procedure as mentioned in the previous part for constructing the weight function. Therefore, the test function for MLPG1 is defined as 0 1 0 1 0 1 dIx A @ dIy A @ dIz A I @ nðxÞ ¼ c ¼ f ð17Þ f f ¼ nx ny nz I I I as d x as :d y as :d z
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
The node I is located at the center of the local sub-domain, therefore, the support of the test function equals to the size of the brick-shaped sub-domain. The LSWF stated in Eq. (5) can be written as Z
Z
sij nI;j dO I
Os
Z
GIsu
Z ti nI dG þ a Z
t j nI dG þ a
¼ GIst
GIsu
GIsu
u i nI dG þ
ui nI dG Z OIs
2
0
with (
bi nI dO
ð18Þ
Si ¼
8 9 8 J 9 8 9 J 3 > u^ > u^ > > > > 0 > > > < b1 > = < = = > < x> J J J 07 v^ ; b ¼ b2 ; 5; u^ ¼ u^ y ¼ > > > > > > > :b > ; > > :w Sz ; > : u^ J > 3 ^J;
0 Sy
Sx 6 S¼4 0
567
0
ð23aÞ
z
1
if ui is prescribed on Gu
0
if ui is not prescribed on Gu
i ¼ x; y; z
;
and 2 In the concept of MLPG5, the Heaviside step function (see Eq. (19)) is chosen as the test function. (
nðxÞ ¼
1
at x A Os
0
ð19Þ
at x= 2Os
D11
6 6 D12 6 EðzÞð1nÞ 6 6 D12 D ¼ DðzÞ ¼ 6 ð1 þ nÞð12nÞ 6 0 6 6 0 4 0
D12
D12
0
0
D11
D12
0
0
D12
D11
0
0
0
0
D22
0
0
0
0
D22
0
0
0
0
0
3
7 0 7 7 0 7 7 7; 0 7 7 0 7 5 D22
n
12n with D11 ¼ 1; D12 ¼ : ; D22 ¼ 1n 2ð1nÞ 2.4. Discretization of the weak form and numerical implementation for the 3D plate problem Substituting the MLS appproximation function, Eq. (10) into Eq. (5) and summing up of all nodes, the following discretized system of linear equations is obtained. The discretized system for the MLPG1 method can be written as Z OIs
M X
BI
J
DBJ u^ dO þ a
Z
J¼1
Z ¼
GIst
WI t d G þ a
Z GIst
GIsu
WI S
M X
J
UJ u^ dG
Z GIsu
J¼1
WI u dG þ
Z
WI SND
M X
J
BJ u^ dG
J¼1
WI bdO
OIst
ð20Þ
Simplifying Eq. (20), the following system of linear algebraic J equations in terms of u^ is obtained. M X
J KIJ u^ ¼ f I ;
I ¼ 1; 2; . . . ; M:
ð21Þ
J¼1
where M is the total number of nodes. The so-called ‘‘stiffness’’ matrix K and ‘‘load’’ vector f are defined as Z Z Z BI DBJ dO þ a WI SUJ dG WI SNDBJ dG ð22aÞ KIJ ¼ OIs
fI ¼
Z
GIsu
Z
GIst
WI t dG þ a
GIsu
GIsu
WI u dG þ
Z OIs
WI b dO
where nx, ny and nz represent the component of outward normal vector on the boundary of the local sub-domain. E and n are the Young’s modulus and the Poisson’s ratio, respectively. The Young’s modulus is assumed to vary exponentially through the thickness according to EðzÞ ¼ Eh egðz=h1Þ ;
fJ;x
6 6 0 6 6 6 0 6 J B ¼6 6 fJ 6 ;y 6 6 6 0 4
0
fJ;x fJ;z
fJ;z
2
7 2 I 0 7 7 7 c;x fJ;z 7 7 I 6 7; B ¼ 6 6 0 4 0 7 7 7 0 fJ;y 7 7 5
0
g ¼ lnðE0 =Eh Þ
ð25Þ
where E0 = E(0) and Eh = E(h). Also, u, v and w are displacement components of a point along the x, y and z axes, respectively and cI is the test function for the MLPG1, which was defined in Eq. (17). By using relations (10) and (12), the unknown variables and their derivatives are obtained. Considering the MLPG5 method, the stiffness matrix and the load vector are given by Z KIJ ¼
ð22bÞ fI ¼
Z GIst
GIsi
Z NDBJ dG
Z t dG þ a
GIsu
GIsu
Z SNDBJ dG þ a
u dG þ
Z OIs
GIsu
SUJ dG
ð26aÞ
b dO
ð26bÞ
O 0
0
cI;y
cI;y
0
cI;x
cI;z
0
cI;z
0
cI;y
0
cI;z
3
nx
0
0
ny
0
6 N¼4 0 0
ny
0
nx
nz
0
nz
0
ny
nz
3
a
b
cI;x
z
y Side 4
2
fJ
6 J 07 5;U ¼ 6 4 0 nx 0
0
fJ 0
0
3
2
Side 2
Side 1
7 0 7 7; 5
fJ;x
0
ð24Þ
3
0
fJ;y
0 r z rh
where g is the inhomogeneity parameter. Poisson’s ratio is considered to be constant. Since a FGM is assumed to be isotropic and the Poisson’s ratio to be constant, it is possible to express the inhomogeneity parameter g in terms of the values of the Young’s modulus E= E(z) on the top (z= 0) and the bottom (z =h) surfaces of the plate as
where in 3D space, the matrices and the vectors in Eqs. (22) can be expressed as 2
ð23bÞ
cI
6 7 I 6 0 7 5; W ¼ 4 0 J f 0
0
cI 0
0
c
E0
3
7 0 7 5; I
x Side 3
h
z
Fig. 2. 3D coordinate system and problem dimensions.
Eh
ARTICLE IN PRESS 568
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
3. Defining the problem domain and boundary conditions A rectangular Cartesian coordinate system, shown in Fig. 2, is used to describe the geometry of a rectangular elastic plate with the volume O defined by 0oxoa, 0oyob and 0ozoh. The axes of x and y are positioned in the top surface while, z axis is directed along the thickness. Nodes are located in all three directions to discretize the plate continuum. By the present method, unequal distribution of nodes in different directions is possible. Therefore, nodes can be added in each direction easily to obtain desirable accuracy. The edges x=0,a and y=0,b are boundaries in which the essential boundary conditions or the natural boundary conditions are defined. The top and bottom surfaces of the plate (z=0,h) are loading surfaces. Four types of edge boundary conditions are considered, which are nominated as the simply supported (S), clamped (C), free (F) and guided (G) boundary conditions. These boundary conditions on the edges of the plate are defined as follows:
(ii) Clamped edge (C) u ¼ v ¼ w ¼ 0 on x ¼ 0ðor x ¼ aÞ; or y ¼ 0ðor y ¼ bÞ: (iii) Free edge (F)
sxx ¼ sxy ¼ sxz ¼ 0 on x ¼ 0ðor x ¼ aÞ; or syy ¼ syx ¼ syz ¼ 0 on y ¼ 0ðor y ¼ bÞ: sxy ¼ sxz ¼ 0; u ¼ 0; on x ¼ 0ðor x ¼ aÞ; or syx ¼ syz ¼ 0; v ¼ 0; on y ¼ 0ðor y ¼ bÞ:
1.04 1.03 1.02
on y ¼ 0ðor y ¼ bÞ:
ð27aÞ
1.03
w/w*
1.01 w/w*
1.04
1 0.99
1.02
0.98
1.01
0.97
1
0.96
0.99
0.95 γ =10-1 (MLPG1) γ =10-1 (MLPG5) γ =-10-1 (MLPG1) γ =-10-1 (MLPG5) γ =0 (MLPG1) γ =0 (MLPG5)
0.98 0.97 0.96 5x5x5
7x7x7
g
3D elasticity solution [10]
3D MLPG1
Error (%)
3D MLPG5
Error (%)
10 1 10 2 10 3 10 4 10 5 10 6 10 10 10 10 10 10 0
1.4146 1.3496 1.3433 -1.3426 1.3426 1.3426 1.2740 1.3355 1.3419 1.3425 1.3425 1.3425 1.3430
1.4096 1.3448 1.3384 1.3378 1.3377 1.3377 1.2693 1.3312 1.3370 1.3377 1.3377 1.3377 1.3377
0.35 0.36 0.36 0.36 0.36 0.36 0.37 0.32 0.37 0.36 0.36 0.36 0.39
1.4229 1.3573 1.3509 1.3503 1.3502 1.3502 1.2811 1.3432 1.3495 1.3502 1.3502 1.3502 1.3502
0.58 0.57 0.57 0.57 0.57 0.57 0.56 0.58 0.57 0.57 0.57 0.57 0.54
2 3 4 5 6
13x13x5
13x13x7 13x13x9 13x13x11 13x13x13 13x13x15 Nodal density
0
9x9x9 11x11x11 13x13x13 15x15x15 Nodal density
-1
Eh/E0:
-2
0.1 Analytical[10] 0.1 Present(MLPG1) 0.1 Present(MLPG5) 10 Analytical [10] 10 Present(MLPG1) 10 Present(MLPG5) 1 Analytical [10] 1 Present(MLPG1) 1 Present(MLPG5)
⎯w
Table 1 Non-dimensional out-of-plane displacement wða=2; b=2; h=2Þ in FG square plates with SSSS boundary conditions and SL loading condition. (a/h= 3).
γ= 10-1 (MLPG1) γ= 10-1 (MLPG5) γ= -10-1 (MLPG1) γ= -10-1 (MLPG5) γ= 0 (MLPG1) γ= 0 (MLPG5)
Fig. 4. Convergence of the results for the central deflection of an SSSS square plate (a/h= 3) for various nodal densities along the thickness direction.
Fig. 3. Convergence of the results for the central deflection of an SSSS square plate (a/h= 3) for various nodal densities.
1
ð27dÞ
For example, the SFCG plate is simply supported, free, clamped and guided along the sides 1, 2, 3 and 4, respectively (see Fig. 2).
sxx ¼ 0; v ¼ w ¼ 0 on x ¼ 0ðor x ¼ aÞ; or syy ¼ 0;
0.95
ð27cÞ
(iv) Guided edge (G)
(i) Simply supported edge (S)
u¼w¼0
ð27bÞ
-3
-4
-5
0
0.25
0.5
0.75
1
⎯z Fig. 5. Variation in non-dimensional deflection wð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
The components of the stress tensor are related to components of strain tensor as follows: n o rT ¼ sxx syy szz syz szx sxy ¼ De ð28Þ In Eq. (28), e denotes the infinitesimal strain with the following components 8 9 > > fI;x u^ I 8 9 > > > > > > exx > > > > > > > > I I > > > > ^ > > > > f v > eyy > > > ;y > > > > > > > > > > > > I I > > > > N < ezz = X < = ^ f;z w e ¼ 2e ¼ ð29Þ I I I I > > > > ^ ^ yz f þ f w v > > ;y > > ;z > > > > > > I¼1> > > > > > 2ezx > > > I I I > > > > ^I> > > > > > > f;z u^ þ f;x w > > : > > > > 2exy ; > I > I I I > : f v^ þ f u^ > ; ;x
569
In the current study two different loading conditions are considered
Uniform load q distributed on the top surface (UL) szz ¼ q;
szx ¼ szy ¼ 0 on z ¼ 0szz ¼ 0; szx ¼ szy ¼ 0 on z ¼ h
ð30aÞ
Sinusoidal load with the intensity q distributed on the bottom surface (SL)
szz ¼ 0;
szx ¼ szy ¼ 0 on z ¼ 0 px py sin ; szx ¼ szy ¼ 0 on z ¼ h
szz ¼ q sin
a
b
ð30bÞ
;y
-0.9 4
-0.8
3
-0.7
2 -0.6
⎯σxz
⎯σxx
1 0
Eh/E0: -1
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-2 -3 -4
0.25
0.5
Eh/E0:
-0.4
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-0.3 -0.2 -0.1 0
-5 0
-0.5
0.75
0
0.25
1
⎯z Fig. 6. Variation in non-dimensional in-plane normal stress s xx ð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
0.5
0.75
1
⎯z Fig. 8. Variation in non-dimensional transverse shear stress s xz ð0; 0:5b; zÞ versus z of an SSSS square plate.
2
0
1
-0.25
⎯σzz
⎯σxy
0 Eh/E0: -1
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-2
-3 0
0.25
0.5
0.75
-0.5 Eh/E0: 0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-0.75
-1
1
⎯z Fig. 7. Variation in non-dimensional in-plane shear stress s xy ð0; 0; zÞ versus z of an SSSS square plate.
-1.25
0
0.25
0.5
0.75
1
⎯z Fig. 9. Variation in non-dimensional out-of-plane normal stress s zz ð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
ARTICLE IN PRESS 570
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
4. Numerical results and discussion
4.1. Comparison studies
Several examples have been studied to examine the accuracy of the presented method in 3D analysis of thick FG rectangular plates with various boundary conditions. It is assumed that the plate is subjected to uniform or sinusoidal distributed loads (UL and SL loading conditions) and the body forces are assumed to be 0. Also, the Poisson’s ratio is considered to be 0.3. In the following analysis, dimensionless size of the sub-domain and support domain are considered as as = 0.75 and ai = 2.6. However, the effect of the value of the parameters as and ai on the accuracy of the solution is studied by Vaghefi et al. [42]. It should be noted that the number of gauss points is considered to be 216 in the local sub-domain integration. To present the obtained results the following non-dimensional parameters are introduced.
In order to examine the efficiency and accuracy of the presented approach, the obtained results are compared with the available exact 3D solution given by Kashtalyan [10]. First, the effects of nodal density on the obtained results for an FG square plate with simply supported boundary conditions (SSSS) are examined and presented in Fig. 3. Various nodal densities containing 5 5 5, 7 7 7, y, 15 15 15 regularly placed along the x-, y- and z-directions have been considered. In Fig. 3, w* refers to the exact solution obtained by Kashtalyan [10] for a SSSS FG square plate with a/h =3. As shown in the figure, it is observed that by choosing 13 13 13 nodes (2197 nodes), the solutions are converged to the exact ones for all the cases. Also, the effect of nodal density on the plate thickness direction, the direction in which the material properties change, is examined. The results for a square FG plate with SSSS boundary conditions with a/h= 3 are depicted in Fig. 4. In this figure, various nodal densities containing 13 13 5, 13 13 7, y, 13 13 15 nodes have been considered. According to the figure it can be observed that considering 13 and more nodes along the
w¼
sij Gh w z ; z¼ ; s ij ¼ q hq h
ð31Þ
In Eq. (31), the parameter Gh is the shear modulus on the edge z= h of the plate.
Table 2 Numerical results for FG square plates with various boundary conditions and UL loading condition (h/a= 0.2). Eh/E0
s xx (a/2, b/2, 0)
w max MLPG1
MLPG5
FEM
s yy (a/2, b/2, 0)
MLPG1
MLPG5
FEM
MLPG1
MLPG5
FEM
(a) SSSS 0.1 0.5 1 2 10
4.0916 8.9751 12.5997 17.6640 39.0605
4.1596 8.9357 12.6375 17.8397 39.0385
4.1215 9.0047 12.6134 17.7118 39.1558
15.3564 9.2902 7.4462 5.9410 3.4665
15.4390 9.3279 7.4199 5.9711 3.4944
15.403 9.2995 7.4588 5.9591 3.4805
15.3564 9.2902 7.4462 5.9410 3.4665
15.4390 9.3279 7.4199 5.9711 3.4944
15.4032 9.2995 7.4588 5.9591 3.4805
(b) SCSC 0.1 0.5 1 2 10
2.4311 5.4661 7.7215 10.8228 23.5185
2.4270 5.5147 7.7419 10.8902 23.6686
2.4373 5.4851 7.7336 10.8517 23.6058
9.1729 5.6761 4.6394 3.7090 2.2461
9.0727 5.7029 4.6203 3.7128 2.2267
9.1343 5.6959 4.6342 3.7334 2.2364
11.1523 6.7495 5.4580 4.3708 2.5759
11.2305 6.7348 5.4618 4.3291 2.5517
11.1814 6.7647 5.4505 4.3464 2.5650
(c) SSSF 0.1 0.5 1 2 10
12.4349 27.2606 38.4889 54.3348 122.7609
12.4062 27.5442 38.3326 54.9850 124.3991
12.5132 27.4253 38.6582 54.6430 123.4156
25.1382 15.1983 12.2387 9.6918 5.5863
25.4982 15.4213 12.3027 9.6697 5.6536
25.3471 15.2947 12.2592 9.7428 5.6276
11.9221 7.1349 5.7356 4.5566 2.7440
11.7441 7.0864 5.6989 4.5886 2.7426
11.8456 7.1255 5.7244 4.5778 2.7295
(d) SSSG 0.1 0.5 1 2 10
9.5911 20.9661 29.4798 41.8456 94.0184
9.7573 21.1375 29.4401 41.5796 92.9495
9.6649 21.0676 29.5869 41.7467 93.5501
25.3203 15.4382 12.3969 9.7842 5.7340
25.7049 15.2773 12.2797 9.7465 5.6173
25.5384 15.3930 12.3336 9.8204 5.6623
14.7404 8.8260 7.1261 5.7314 3.3589
14.5421 8.8355 7.0901 5.8132 3.3819
14.6202 8.8759 7.1448 5.7790 3.3731
(e) CFCF 0.1 0.5 1 2 10
3.5627 8.1319 11.6473 16.3762 36.9346
3.5879 8.1895 11.5694 16.3885 36.4419
3.5720 8.1516 11.6298 16.4360 36.7001
13.5523 8.1803 6.5891 5.2458 3.0279
13.3829 8.1169 6.5536 5.1909 3.0725
13.4943 8.1593 6.5822 5.2223 3.0486
3.3407 2.0525 1.6935 1.3969 0.94076
3.2949 2.0471 1.6802 1.3825 0.94204
3.3156 2.0398 1.6874 1.3884 0.93952
(f) CGCG 0.1 0.5 1 2 10
3.2034 7.2736 10.2880 14.6290 31.5918
3.2327 7.3329 10.3664 14.3906 31.8806
3.2145 7.2853 10.3357 14.5306 31.9116
13.9612 8.3649 6.8187 5.4403 3.1722
13.9381 8.4862 6.7709 5.3926 3.1434
13.8783 8.4172 6.8003 5.4060 3.1662
4.4350 2.8342 2.3295 1.9301 1.2416
4.4308 2.8105 2.3406 1.9084 1.2555
4.4510 2.8174 2.3374 1.9167 1.2468
(g) CCCC 0.1 0.5 1 2 10
1.7485 3.9547 5.6094 7.8809 16.9989
1.7457 3.9251 5.6221 7.8188 17.1946
1.7340 3.9499 5.5903 7.8589 17.1054
8.0702 5.0016 4.0694 3.2773 1.9772
8.1498 4.9610 4.0265 3.2334 1.9517
8.1096 4.9837 4.0535 3.2548 1.9706
8.0702 5.0016 4.0694 3.2773 1.9772
8.1498 4.9610 4.0265 3.2334 1.9517
8.1096 4.9837 4.0535 3.2548 1.9706
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
571
plate thickness direction can led to an acceptable solution. Therefore, to obtain the solutions in the following problems, the nodal distribution is considered to be 13 13 13. In Table 1, the non-dimensional out-of-plane displacement w of an SSSS FG square plate is given. Assuming a/h= 3 and a sinusoidal loading (see Eq. (30b)), various inhomogeneity parameters g are considered and the results are calculated in the plate center (x =0.5a, y=0.5b, z =0.5h). It should be mentioned that the load is assumed to be applied on the bottom (z= h) surface of the plate. It is clear while considering g = 0, the obtained solutions can be related to the same homogeneous plate. The results obtained by using both the methods, MLPG1 and MLPG5, are presented in the table. In Table 1 the percent of error is calculated using the following equation Present solutionExact solution ð32Þ 100 %Error ¼ Exact solution
Also, the variations in the non-dimensional stresses s xx , s xy , s xz and s zz are depicted through the thickness of the plate in Figs. 6–9, respectively. Different values of Eh/E0 are considered in plotting the figures and both the methods MLPG1 and MLPG5 are applied to obtain the results, which are compared with the reference solution [10]. In plotting Figs. 6–9, the plate is assumed to be simply supported with a/h= 3, under a sinusoidal loading at z= h. Numerical results in Figs. 6–9 are given for three different values of Eh/E0. The chosen values of the inhomogeneity ratio Eh/E0 = 10, 1.0, and 0.1 do not necessarily represent a certain material, but are rather used to show the effect of inhomogeneity on stress and displacement fields. It can be seen that the present results agree well with those reported by Kashtalyan [10].
where the exact solution refers to the analytic solution w* for an SSSS thick FG square plate with a/h= 3 given in Ref. [10]. In Fig. 5, variation in the non-dimensional out-of-plane displacement w is plotted through the thickness of the plate.
The FG plates with various boundary conditions are analyzed using the MLPG1 and MLPG5 methods. In Table 2, the maximum non-dimensional deflection w max and the non-dimensional normal stresses s xx and s yy at the center of the top surface
4.2. Numerical examples with various boundary conditions
Table 3 Numerical results for FG square plates with various boundary conditions and UL loading condition. (h/a =0.3). Eh/E0
s xx (a/2, b/2, 0)
w max MLPG1
MLPG5
FEM
MLPG1
s yy (a/2, b/2, 0) MLPG5
FEM
MLPG1
MLPG5
FEM
(a) SSSS 0.1 0.5 1 2 10
0.9707 2.1378 2.9853 4.1208 8.7134
0.9688 2.1498 2.9603 4.1098 8.6923
0.9732 2.1407 2.9792 4.1333 8.7293
7.2230 4.3084 3.4496 2.7499 1.6449
7.2034 4.2943 3.4959 2.7556 1.6566
7.2639 4.3378 3.4681 2.7673 1.6499
7.2230 4.3084 3.4496 2.7499 1.6449
7.2034 4.2943 3.4959 2.7556 1.6566
7.2639 4.3378 3.4681 2.7673 1.6499
(b) SCSC 0.1 0.5 1 2 10
0.6637 1.5302 2.1434 2.9696 6.3437
0.6675 1.5286 2.1555 3.0066 6.3506
0.6655 1.5269 2.1481 2.9890 6.3713
4.8701 3.0584 2.4868 2.0094 1.2493
4.8784 3.0311 2.4692 2.0085 -1.2412
4.8640 3.0484 2.4794 2.0125 1.2466
5.5628 3.3449 2.6511 2.1526 1.3128
5.5672 3.3437 2.6555 2.1543 1.3103
5.5558 3.3347 2.6608 2.1435 1.3073
(c) SSSF 0.1 0.5 1 2 10
2.7629 6.0992 8.5891 12.0999 26.9454
2.7597 6.1412 8.5866 12.1041 26.7625
2.7648 6.1154 8.6239 12.1192 26.8886
11.7326 7.0170 5.6075 4.4438 2.5999
11.7563 7.0088 5.6392 4.4358 2.5618
11.8060 7.0433 5.6250 4.4592 2.5836
5.3477 3.1655 2.5304 2.0512 1.2726
5.3210 -3.1502 2.5256 2.0381 1.2649
5.3603 3.1771 2.5496 2.0470 1.2677
(d) SSSG 0.1 0.5 1 2 10
2.1435 4.6698 6.5596 9.1481 20.0863
2.1407 4.6802 6.5399 9.1548 19.9522
2.1500 4.7085 6.5882 9.2033 20.0389
11.7572 6.9951 5.5844 4.4223 2.5761
11.7181 7.0800 5.6418 4.4690 2.5881
11.8310 7.0385 5.6143 4.4504 2.5790
6.8423 4.1385 3.3613 2.6698 1.6293
6.8262 4.1289 3.3419 2.6712 1.6251
6.8781 4.1569 3.3496 2.6892 1.6316
(e) CFCF 0.1 0.5 1 2 10
0.9624 2.2818 3.2947 4.7104 10.6373
0.9765 2.2871 3.3346 4.7115 10.6248
0.9686 2.2986 3.3229 4.7246 10.6851
6.6913 4.0416 3.2433 2.5593 1.5300
6.7873 4.0746 3.2478 2.5907 1.5260
6.7453 4.0482 3.2368 2.5790 1.5344
1.5581 0.94903 0.78921 0.68142 0.52378
1.5632 0.95112 0.78725 0.67291 0.52526
1.5515 0.95299 0.79313 0.67653 0.52173
(f) CGCG 0.1 0.5 1 2 10
0.8824 2.0743 2.9257 4.1037 8.8568
0.8798 2.0583 2.9246 4.1344 8.8598
0.8845 2.0650 2.9422 4.1138 8.8953
6.8910 4.1806 3.3383 2.6878 1.5976
6.9751 4.2011 3.3311 2.6746 1.5993
6.9342 4.1858 3.3568 2.6842 1.6062
2.3574 1.5543 1.2956 1.0974 0.77678
2.3489 1.5453 1.3128 1.1092 0.78193
2.3724 1.5510 1.3056 1.1022 0.78012
(g) CCCC 0.1 0.5 1 2 10
0.4948 1.1649 1.6554 2.3312 5.0918
0.5016 1.1753 1.6521 2.3469 5.1350
0.4985 1.1702 1.6626 2.3368 5.1065
4.2097 2.6062 2.0965 1.7194 1.0917
4.2029 2.5950 2.0944 1.7318 1.1040
4.2293 2.6120 2.1097 1.7277 1.0974
4.2097 2.6062 2.0965 1.7194 1.0917
4.2029 2.5950 2.0944 1.7318 1.1040
4.2293 2.6120 2.1097 1.7277 1.0974
ARTICLE IN PRESS 572
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
Fig. 10, it can be seen that the through-thickness variation in the non-dimensional out-of-plane displacement w in the FG plates is not constant and it has the maximum value near the top surface for all values of stiffness gradient Eh/E0. It should be mentioned that with increase in the ratio Eh/E0, the out-of-plane displacement w of the FG plates decreases. Also in Fig. 11, variation in the in-plane normal stress s xx through the thickness of a clamped FG plate under sinusoidal loading is shown. The plate has the same thickness ratio as in Fig. 10. According to Fig. 11, it can be seen that in FG plates, the maximum of the in-plane normal stress s xx is located near the surface with lower Young’s modulus, i.e. the maximum compressive stress happened near the top surface of the FG plate with Eh/E0 = 0.1 and the maximum of the tensile stress happened near the bottom surface of the FG plate with Eh/E0 =10.
0
⎯w
-1
Eh/E0: 0.1 0.5 1 2 10
-2
5. Conclusions
-3
0
0.25
0.5
0.75
1
⎯z Fig. 10. Variation in non-dimensional deflection wð0:5a; 0:5b; zÞ versus z of a CCCC square plate.
3
2
⎯σxx
1
0
-1
Eh/E0: 0.1 0.5 1 2 10
-2
-3 0
0.25
A 3D and completely meshless approach have been developed for the static analysis of thick FG plates with various boundary conditions. The 3D equilibrium equations have been considered and the Young’s modulus of the plate was assumed to vary exponentially through the thickness, while the other material properties were assumed to be constant. The two meshless methods including MLPG1 and MLPG5 have been extended to solve the bending problem of FG plates. Brick-shaped local domains were considered in the solution procedure. The application of brick-shaped local domains facilitates the integration and the approximation procedure. The proposed approach for construction of the shape and the test functions makes it possible to introduce more nodes in the direction of material variation. Consequently, more precise solutions can be obtained easily and efficiently. It is observed that the approach works well in analysis of FG plates. Several numerical examples have been presented. The obtained results have been compared with the available analytical and numerical solutions in the literature and an excellent agreement was seen. Although the approach has been applied for static analysis, it can be used for the dynamics, buckling and fracture analyses of thick FG plates as well.
References
0.5 ⎯z
0.75
1
Fig. 11. Variation in non-dimensional in-plane normal stress s xx ð0:5a; 0:5b; zÞ versus z of a CCCC square plate.
(z= 0) are given for square FG plates with various boundary conditions including SSSS, SCSC, SSSF, SSSG, CFCF, CGCG and CCCC. For the results presented in Table 2, the UL loading condition is assumed and the plates with the thickness ratios of h/a =0.2 and Eh/E0 = 0.1, 0.5, 1, 2, and 10 are considered. To perform a better comparison, the FEM results obtained by the authors are also presented. The FEM results are obtained by using the ANSYS software with 36 36 36 solid 45 brick elements. Table 3 presents the same results as Table 2 while considering the thickness ratio of h/a= 0.3. In Fig. 10, the through-thickness variation in the nondimensional out-of-plane displacement w of a CCCC FG plate under sinusoidal loading is shown. The results are obtained considering a/h=3 and different values of Eh/E0. According to
[1] Igari T, Notomi A, Tsunoda H, Hida K, Kotoh T, Kunishima S. Material properties of functionally gradient material for fast breeder reactor. In: Yamanouochi M, Koizumi M, Hirai T, Shiota I, editors. Proceedings of the first international symposium on functionally gradient materials, Japan: Sendai; 1990. p. 209–214. [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 Appl Biomater 1994;5(1):23–7. [3] Tani J, Liu GR. Surface waves in functionally gradient piezoelectric plates. JSME Int J Ser A (Mech Mater Eng) 1993;36(2):152–5. [4] Osaka T, Matsubara H, Homma T, Mitamura S, Noda K. Microstructural study of electroless-plated CoNiReP/NiMoP doublelayered media for perpendicular magnetic recording. Jpn J Appl Phys 1990;29(10):1939–43. [5] Watanabe Y, Nakamura Y, Fukui Y, Nakanishi K. Amagnetic-functionally graded material manufactured with deformation induced martensitic transformation. J Mater Sci Lett 1993;12(5):326–8. [6] Getto H, Ishihara S. Development of the fire retardant door with functionally gradient wood. In: Proceedings of the fourth international symposium on functionally graded materials, Japan: Tsukuba City; 1996. [7] Levinson M. The simply supported rectangular plate: an exact, three dimensional, linear elasticity solution. J Elasticity 1985;15:283–91. [8] Batra RC, Vel SS. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2001;40(7):1421–33. [9] Zhong Z, Shang ET. Three-dimensional exact analysis of a simply supported functionally gradient piezoelectric plate. Int. J. Solids Struct 2003;40(20):5335–52. [10] Kashtalyan M. Three-dimensional elasticity solution for bending of functionally graded rectangular plates. Eur J Mech A—Solid 2004;23:853–64.
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
[11] Liew KM, Teo TM. Modeling via differential quadrature method: threedimensional solutions for rectangular plates. Comput Methods Appl Mech Eng 1998;159:369–81. [12] Liew KM, Teo TM, Han JB. Three-dimensional solutions of rectangular plates by variant differential quadrature method. Int J Mech Sci 2001;43:1611–28. ¨ . Three-dimensional vibration, buckling and bending analysis of [13] Civalek O thick rectangular plates based on discrete singular convolution method. Int J Mech Sci 2007;49:752–65. [14] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [15] Liu WK, Jun S, Zhang Y. Reproducing kernel particle methods. Int J Numer Methods Fluids 1995;20:1081–106. [16] Duarte CA, Oden JT. Hp-cloud-a meshless method to solve boundary-value problems. Comput Methods Appl Mech Eng 1996;139:237–62. [17] Babuska I, Melenk J. The partition of unity method. Int J Numer Methods Eng 1997;40:727–58. [18] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [19] Atluri SN, Kim HG, Cho JY. A critical assessment of the truly meshless local Petrov–Galerkin (MLPG), and local boundary integral equation (LBIE) methods. Comput Mech 1999;24:348–72. [20] Atluri SN, Shen S. The meshless local Petrov–Galerkin (MLPG) method: a simple and less-costly alternative to the finite element and boundary element methods. CMES: Comput Modeling Eng Sci 2002;3(1):11–51. [21] Atluri SN, Shen S. The meshless LOCAL Petrov–Galerkin (MLPG) method. USA: Tech Science Press; 2002. [22] Atluri SN, Zhu T. The meshless local Petrov–Galerkin (MLPG) approach for solving problems in elasto-statics. Comput Mech 2000;25:169–79. [23] Batra RC, Ching HK. Analysis of elastodynamic deformation near a crack/ notch tip by the meshless local Petrov–Galerkin (MLPG) Method. CMES: Comput Modeling Eng Sci 2002;3(6):717–30. [24] Lin H, Atluri SN. The meshless local Petrov–Galerkin (MLPG) method for solving incompressible Navier–Stokes equations. CMES: Comput Modeling Eng Sci 2001;2(2):117–42. [25] Lin H, Atluri SN. Meshless local Petrov–Galerkin (MLPG) method for convection-diffusion problems. CMES: Comput Modeling Eng Sci 2000;1(2):45–60. [26] Sladek J, Sladek V, Atluri SN. A pure contour formulation for the meshless local boundary intergral equation method in thermoelasticity. CMES: Comput Modeling Eng Sci 2001;2(4):423–33. [27] Atluri SN, Cho JY, Kim HG. Analysis of thin beams, using the meshless local Petrov–Galerkin method, with generalized moving least squares interpolations. Comput Mech 1999;24:334–47. [28] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) formulation for static and free vibration analyses of thin plates. CMES: Comput Modeling Eng Sci 2001;2(4):463–76.
573
[29] Long S, Atluri SN. A meshless local Petrov–Galerkin (MLPG) method for solving the bending problem of a thin plate. CMES: Comput Modeling Eng Sci 2002;3(1):53–64. [30] 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. [31] Xia P, Long SY, Cui HX, Li GY. The static and free vibration analysis of a nonhomogeneous moderately thick plate using the meshless local radial point interpolation method. Eng Anal Boundary Elem 2009;33(6):770–7. [32] Sladek J, Sladek V, Zhang Ch. Stress analysis in anisotropic functionally graded materials by the MLPG method. Eng Anal Boundary Elem 2005;29:597–609. [33] Wang H, Qin QH. Meshless approach for thermo-mechanical analysis of functionally graded materials. Eng Anal Boundary Elem 2008;32(9): 704–12. [34] Kim HG, Atluri SN. Arbitrary placement of secondary nodes, and error control, in the meshless local Petrov–Galerkin (MLPG) method. CMES: Comput Modeling Eng Sci 2000;1(3):11–32. [35] Ching HK, Batra RC. Determination of crack tip fields in linear elastostatics by the meshless local Petrov–Galerkin (MLPG) method. CMES: Comput Model Eng Sci 2001;2(2):273–90. [36] Koohkan H, Baradaran GH, Vaghefi R. A completely meshless analysis of cracks in isotropic functionally graded materials. In: Proceedings of the IMechE, Part C: J Mech Eng Sci 2009, in press. [37] Tang Z, Shen S, Atluri SN. Analysis of materials with strain-gradient effects: A meshless local Petrov–Galerkin (MLPG) approach, with nodal displacements only. CMES: Comput Modeling Eng Sci 2003;4(1): 177–96. [38] Li Q, Shen S, Han ZD, Atluri SN. Application of meshless local Petrov–Galerkin (MLPG) to problems with singularities, and material discontinuities, in 3-D Elasticity. CMES: Comput Modeling Eng Sci 2003;4(5):571–85. [39] Han ZD, Atluri SN. Truly meshless local Petrov–Galerkin (MLPG) solutions of traction and displacement BIEs. CMES: Comput Modeling Eng Sci 2003;4(6):665–78. [40] Han ZD, Atluri SN. Meshless local Petrov–Galerkin (MLPG) approaches for solving 3-D problems in elasto-statics. CMES: Comput Modeling Eng & Sci 2004;6(2):169–88. [41] Han ZD, Atluri SN. A meshless local Petrov–Galerkin (MLPG) approach for 3dimensional elasto-dynamics. CMC: Comput Mater Continua 2004;1(2): 129–40. [42] Vaghefi R., Baradaran G.H., Koohkan H., Three-dimensional static analysis of rectangular thick plates by using meshless local Petrov–Galerkin method. In: Proceedings of the IMechE, Part C: J. Mech Eng Sci, vol 223, C9; 2009. p. 1983–1996
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Three-dimensional static analysis of thick functionally graded plates by using meshless local Petrov–Galerkin (MLPG) method R. Vaghefi, G.H. Baradaran , H. Koohkan Department of Mechanical Engineering, Shahid Bahonar University of Kerman, Kerman, Iran
a r t i c l e in f o
a b s t r a c t
Article history: Received 30 July 2009 Accepted 20 January 2010 Available online 20 February 2010
In this paper, a version of meshless local Petrov–Galerkin (MLPG) method is developed to obtain threedimensional (3D) static solutions for thick functionally graded (FG) plates. The Young’s modulus is considered to be graded through the thickness of plates by an exponential function while the Poisson’s ratio is assumed to be constant. The local symmetric weak formulation is derived using the 3D equilibrium equations of elasticity. Moreover, the field variables are approximated using the 3D moving least squares (MLS) approximation. Brick-shaped domains are considered as the local sub-domains and support domains. In this way, the integrations in the weak form and approximation of the solution variables are done more easily and accurately. The proposed approach to construct the shape and the test functions make it possible to introduce more nodes in the direction of material variation. Consequently, more precise solutions can be obtained easily and efficiently. Several numerical examples containing the stress and deformation analysis of thick FG plates with various boundary conditions under different loading conditions are presented. The obtained results have been compared with the available analytical and numerical solutions in the literature and an excellent consensus is seen. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Three-dimensional meshless local Petrov–Galerkin method Functionally graded material Three-dimensional theory of elasticity Thick plate Moving least squares approximation
1. Introduction Functionally graded materials (FGMs) are advanced composites in which the material properties vary smoothly and continuously by a predetermined function. Generally, those are composed of ceramics and metals so that the material properties such as the Young’s modulus vary from the metallic side to the ceramic one. Initially, they were designed as thermal barrier coatings in space applications. However, recently FGMs have gained lots of applications in nuclear reactors [1], dental and medical implants [2], piezoelectric and thermo electric devices [3–5] and fire retardant doors [6]. Investigation of the mechanical behavior of new materials such as FGMs is an attractive research area in mechanics. Many researchers have studied the mechanical behavior of FGMs having various geometries and loading conditions. Among the geometries, which are considered to analyze, the plates are the most significant due to numerous applications in engineering structures. It should be mentioned that the plate is a three-dimensional (3D) structure in which one dimension is relatively much smaller than the other dimensions. In 2D theories such as the classical plate theory (CPT), in the first order shear deformation plate
Corresponding author. Tel.: + 98 341 2111763; fax: + 98 341 2120964.
E-mail address: [email protected] (G.H. Baradaran). 0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.01.005
theory (FSDT) and the higher order shear deformation plate theory (HSDT), various assumptions are made to obtain a 2D formulation for plates. Obviously, it is easier to find a solution for the plate problems in 2D formulation. However, some errors occur in solutions due to the assumptions considered in 2D theories. With the increase in the thickness of plates, the errors amplify. If no simplifying assumptions are considered, the 3D elasticity equations must be used for elastic analysis of plates. The 3D solution does not involve any limitation of 2D solutions. It is obvious that if 3D solutions are attainable, it will be more accurate than the solutions obtained by the mentioned 2D theories. Analytical and numerical methods have been used for 3D analyses of plates [7–13]. Analytical methods are applicable for some simple cases and generally, the analytical solutions cannot be found for plates with complicated geometries and boundary conditions. Consequently, numerical methods such as the finite element method (FEM) and differential quadrature (DQ) method are the alternatives that have been employed for 3D analysis of plates. Recently, meshless methods, which do not need burdensome effort of mesh generation have gained lots of attention. A variety of these methods have been developed, such as the element-free Galerkin (EFG) method [14], the reproducing kernel particle method (RKPM) [15], hp-clouds [16], and the partition of unity method (PUM) [17], which have been successfully applied for the solutions of some engineering problems. However, most of
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
these methods are not truly meshless, because, they use background cells for the evaluation of integrals in the weak formulation of problems. The meshless local Petrov–Galerkin (MLPG) method proposed by Atluri et al. [18–21] considers local weak formulation instead of the global weak formulation, and does not require a background mesh for evaluation of the integrals. In the MLPG method, the integrals of the weak form are evaluated over local subdomains that partly cover each other. The trial and test functions are chosen from totally different functional spaces. Furthermore, the physical size of the test and trial domains is not necessary to be the same, which makes the MLPG a very flexible method. Based on the concept of the MLPG, six different methods have been introduced, which are labeled as MLPG1–MLPG6 [21]. These six methods differ due to the type of test function considered in the weak formulation. Among them, the MLPG1 and MLPG5 seem to be the most promising formulations; so they are used in the present study. In the MLPG1 a namely fourth-order spline function is considered as test function, while the Heaviside step function is employed as test function in MLPG5. It is noticeable that the MLPG5 does not involve any domain integration or singular integrals. The MLPG methods have been employed in a wide range of applications, for example elasto-statics [22], elastodynamics [23], fluid mechanics [24], convection–diffusion problems [25], thermoelasticity [26], beam problems [27,28], plate problems [29–31], FGM problems [32,33], fracture mechanics [34–36], and strain gradient theory [37]. So far, the applications of the MLPG method are mostly limited to 2D problems. The major reason for the 3D analysis is considered less, related to difficulties in evaluating the integrals over the local sub-domains, especially when a sub-domain intersects the global boundary of problem. However, significant efforts have been carried out to develop the MLPG method to solve 3D problems. Li et al. [38] applied a combination of MLPG2 and MLPG5 for the solution of two classical 3D problems, viz., the Boussinesq problem and the Eshelby’s inclusion problem. They used MLPG5 for nodes inside the domain and MLPG2 for nodes on the boundary of the problem domain. Han and Atluri [39] used the MLPG method to solve 3D elastic fracture problems. Also, Han and Atluri [40] developed three kinds of the MLPG methods using different test functions for analysis of thick beams and spheres and examined the performance of the methods. They used Delaunay triangulation algorithm and introduced a suitable mapping procedure for evaluating the integrals over the spherical sub-domains. The integrals on the volume of the sub-domains were calculated by dividing the domains into cone shapes. The MLPG domain discretization method has been employed to 3D elasto-dynamic problems of impact and fragmentation by Han and Atluri [41]. Vaghefi et al. [42] developed two different 3D MLPG procedures including MLPG1 and MLPG5 for the elasto-static analysis of thick rectangular plates with various boundary conditions. Brick-shaped domains are considered as sub-domains and support domains. The purpose of the present paper is to develop a complete 3D MLPG procedure for the elasto-static analysis of thick FG plates with various boundary conditions. The Young’s modulus of the plate is assumed to vary exponentially through the thickness, and the Poisson’s ratio is assumed to be constant. The local symmetric weak form (LSWF) is used to formulate the problem and 3D MLS approximation is used to approximate the field variables. To impose the essential boundary conditions, the penalty method is adopted. Since the global domain of the problem is cubic, brickshaped domains are considered as sub-domains and support domains for evaluating the integrals of the weak form and approximating the solution variables, respectively. Moreover, the approach used for the construction of the shape and the test functions makes it possible to add nodes in any direction of the plate (through length, width and thickness) depending upon the
565
required accuracy. Therefore, sufficient number of nodes has been added in the thickness direction of the plate to represent the normal and shear stress variations through the thickness with good accuracy. Consequently, more precise solutions can be obtained easily and efficiently. In order to illustrate the accuracy of the proposed approach, several numerical examples are presented and the results are compared with the known solutions in the literature.
2. MLPG formulation for 3D elasticity The 3D elasto-static equilibrium equations in a domain of the volume O, which is bounded by the surface G, are given by
sij;j þ bi ¼ 0; in O
ð1Þ
where sij is the stress tensor and bi is the body force vector. The indices i, j which take the values of 1, 2, 3 refer to the Cartesian coordinates x, y, z, respectively. The boundary conditions are assumed to be ui ¼ u i on Gu ;
ð2aÞ
sij nj ¼ t i on Gt
ð2bÞ
where ui and ti represent the displacement components and the surface traction components, respectively. Gu is the boundary with the prescribed displacement u i while Gt is the boundary with the prescribed traction t i . nj denotes the unit outward vector normal to the boundary G. The weak formulation is constructed over the local subdomains, which are located inside the global domain O. The local sub-domains may overlap with each other and must cover the whole global domain. Various shapes with different sizes can be chosen as sub-domains, but appropriate ones should be considered to obtain precise outcomes. Since the global domain of the rectangular plate is hexahedral, brick-shaped local domains are considered as sub-domain and support domain (see Fig. 1). Selecting the brick-shaped local domain makes the mapping procedure easy and no special treatment is needed when the local sub-domain intersects the global boundary. 2.1. Local symmetric weak form for 3D elastic body The generalized local weak form of the equilibrium equations over a local sub-domain around node I is written as follows: Z Z ðsij;j þ bi ÞnI dOs a ðui u i ÞnI dG ¼ 0 ð3Þ OIs
GIsu
where ui, nI are the trial and the test functions, respectively and can be chosen from different functional spaces. In Eq. (3), the term
local boundary ∂Ωs=Γsi for an internal node
sub-domain Ωs of node I
x
Γu
Γ
Γt
I Γsi
h
z
Γs ∂Ωs=ΓsiUΓs
Ω
y
Fig. 1. Brick-shaped sub-domains in the global domain of a rectangular plate.
ARTICLE IN PRESS 566
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
GIsu illustrates a part of the boundary @OIs over which the essential boundary conditions are prescribed and a denotes a penalty parameter, which is introduced in order to impose the essential boundary conditions (a b1). Using the relation sij,jnI = (sijnI),i sijnI,j and the divergence theorem, Eq. (3) leads to Z Z Z sij nj nI dG ðsij nI;j bi nI Þ dOs a ðui u i ÞnI dG ¼ 0 ð4Þ
where fI(x) is called the MLS approximation shape function and is defined as
It should be noted that the boundary @OIs of the local subdomain is divided into three parts, @OIs ¼ GIsi [ GIsu [ GIst , where GIsi is a part of the local boundary located inside the global domain and has no contact with the global boundary of the problem. Also, GIst is a part of the local boundary that coincides with the global traction boundary and GIsu is a part of the local boundary that coincides with the global displacement boundary. Writing Eq. (4) in terms of GIsi , GIst and GIsu and applying the natural boundary condition ti ¼ sij nj ¼ t i on GIst , the following local symmetric weak form (LSWF) for the linear 3D elasticity problem can be obtained Z Z Z Z sij nI;j dO ti nI dG ti nI dG þ a ui nI dG
uh;i ¼
@OIs
OIs
OIs
GIsi
Z
¼ GIst
GIsu
GIsu
Z
t j nI dG þ a
GIsu
u i nI dG þ
GIsu
Z OIs
bi nI dO
ð5Þ
The test function nI(x) is chosen such that it is positive inside local domain Os and vanishes outside Os.
2.2. 3D approximation using moving least square (MLS) approximation The MLS approximation is considered such that the unknown trial approximant uh(x) of the function u(x) is defined by uh ðxÞ ¼ pT ðxÞaðxÞ
8x A Ox
ð6Þ
T
where p ¼ ½p1 ðxÞ; p2 ðxÞ; . . . ; pm ðxÞ is a complete monomial basis of order m and Ox is the domain of definition of the MLS approximation for the trial function at point x. The quadratic basis vector pT for 3D problems are defined as pT ðxÞ ¼ ½1; x; y; z; x2 ; y2 ; z2 ; xy; yz; zx
quadratic basis
ð7Þ
The coefficient vector a(x) is determined by minimizing a weighted discrete L2-norm, which is defined as JðaðxÞÞ ¼
N X
I
^ T GðxÞ½PaðxÞu ^ gI ðxÞ½pðxI ÞaðxÞu^ 2 ¼ ½PaðxÞu
ð8Þ
fI ðxÞ ¼
m X
pj ðxÞ½A1 ðxÞBðxÞjI
ð11Þ
j¼1
The partial derivatives of the trial function are introduced as follows N X
fI;i ðxÞu^ I
ð12Þ
I¼1 I
where (.),i denotes q(.)/qxi and f;i ðxÞ are derivatives of the MLS shape function and can be obtained as
fI;i ðxÞ ¼
m X
½pj;i ðA1 BÞjI þ pj ðA1 B;i þ A1 ;i BÞjI
ð13Þ
j¼1 1 where A1 Þ;i represents the derivative of the inverse of the ;i ¼ ðA matrix A with respect to xI. A fourth-order spline weight function for a one-dimensional domain is defined as follows: 8 2 3 4 > < 16 dI þ8 dI 3 dI 0 rdI r rI dI rI rI rI ¼ ð14Þ gI ðxÞ ¼ f > rI :0 d I Z rI
In Eq. (14), dI =jx xIj while rI is the size of the support for the I weight function gI defined as rI ¼ ai :d . ai is the dimensionless size I of the support domain and d is the average nodal spacing between two neighboring nodes in the vicinity of node I. Using the brick-shaped support domain, the weight function for the 3D problem can be obtained by a simple extension of the onedimensional function of (14) as follows: ! I dIy d dI ð15Þ gI ðxÞ ¼ f Ix f I f Iz ¼ gx gy gz rx ry rz where the functions gx, gy and gz are obtained by replacing dI/rI with dIx =rxI , dIy =ryI and dIz =rzI in Eq. (14), respectively. The parameters dIx ; rxI ; dIy ; ryI ; dIz and rzI are defined as I
I
I
dIx ¼ jxxI j; dIy ¼ jyyI j; dIz ¼ jzzI j; rxI ¼ ai d x ; ryI ¼ ai :d y ; rzI ¼ ai d z ð16Þ I
I
I
In Eq. (16), d x , d y and d z are the average nodal spacings between two neighboring nodes in the vicinity of the node I in the x-, y- and z-direction, respectively. Then, the dimensions of the support domain used for constructing the MLS shape functions I I I become ð2ai d x Þ ð2ai d y Þ ð2ai d z Þ. 2.3. Test function
I¼1
where gI(x) is the weight function related to node I and gI(x) 40 for all points placed in the support of the weight function. In Eq. I (8), xI is the position of the Ith node, u^ ðI ¼ 1; 2; . . . ; N A Ox Þ is the fictitious nodal value, which may not be equal to uh(xI) and P, G and u^ have the following form 2
2
3 pT ðx1 Þ 6 T 7 6 p ðx2 Þ 7 7 P¼6 6 ^ 7 4 5 pT ðxN Þ
6 6 ; G¼6 6 4
Nm
g1 ðxÞ 0 ^ 0
0 &
& 0
0 ^
3
7 7 7 0 7 5 gN ðxÞ
; u^ ¼ ½u^ 1 ; u^ 2 ; . . . ; u^ N 1N
NN
ð9Þ The coefficient vector a(x) is determined and substituted into I Eq. (6) to obtain the approximation function uh ðu^ Þ as follows: uh ðxÞ ¼ Uu^ ¼
N X I¼1
fI ðxÞu^ I
ð10Þ
The test function is considered to be either a fourth-order spline function or a Heaviside step function and the corresponding MLPG formulations are called MLPG1 and MLPG5, respectively. Brick-shaped sub-domain is also chosen for the support of the test function. In a same manner discussed in the previous section, the dimensions of the sub-domain for node I can be I I I defined as ð2as d x Þ ð2as d y Þ ð2as d z Þ, Where as is a constant, greater than 0 and less than 1, chosen as 0.75 in the present study. Considering the MLPG1 method, the weight function of Eq. (14) has been used to define the test function. The test function for the brick-shaped sub-domain is constructed following the same procedure as mentioned in the previous part for constructing the weight function. Therefore, the test function for MLPG1 is defined as 0 1 0 1 0 1 dIx A @ dIy A @ dIz A I @ nðxÞ ¼ c ¼ f ð17Þ f f ¼ nx ny nz I I I as d x as :d y as :d z
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
The node I is located at the center of the local sub-domain, therefore, the support of the test function equals to the size of the brick-shaped sub-domain. The LSWF stated in Eq. (5) can be written as Z
Z
sij nI;j dO I
Os
Z
GIsu
Z ti nI dG þ a Z
t j nI dG þ a
¼ GIst
GIsu
GIsu
u i nI dG þ
ui nI dG Z OIs
2
0
with (
bi nI dO
ð18Þ
Si ¼
8 9 8 J 9 8 9 J 3 > u^ > u^ > > > > 0 > > > < b1 > = < = = > < x> J J J 07 v^ ; b ¼ b2 ; 5; u^ ¼ u^ y ¼ > > > > > > > :b > ; > > :w Sz ; > : u^ J > 3 ^J;
0 Sy
Sx 6 S¼4 0
567
0
ð23aÞ
z
1
if ui is prescribed on Gu
0
if ui is not prescribed on Gu
i ¼ x; y; z
;
and 2 In the concept of MLPG5, the Heaviside step function (see Eq. (19)) is chosen as the test function. (
nðxÞ ¼
1
at x A Os
0
ð19Þ
at x= 2Os
D11
6 6 D12 6 EðzÞð1nÞ 6 6 D12 D ¼ DðzÞ ¼ 6 ð1 þ nÞð12nÞ 6 0 6 6 0 4 0
D12
D12
0
0
D11
D12
0
0
D12
D11
0
0
0
0
D22
0
0
0
0
D22
0
0
0
0
0
3
7 0 7 7 0 7 7 7; 0 7 7 0 7 5 D22
n
12n with D11 ¼ 1; D12 ¼ : ; D22 ¼ 1n 2ð1nÞ 2.4. Discretization of the weak form and numerical implementation for the 3D plate problem Substituting the MLS appproximation function, Eq. (10) into Eq. (5) and summing up of all nodes, the following discretized system of linear equations is obtained. The discretized system for the MLPG1 method can be written as Z OIs
M X
BI
J
DBJ u^ dO þ a
Z
J¼1
Z ¼
GIst
WI t d G þ a
Z GIst
GIsu
WI S
M X
J
UJ u^ dG
Z GIsu
J¼1
WI u dG þ
Z
WI SND
M X
J
BJ u^ dG
J¼1
WI bdO
OIst
ð20Þ
Simplifying Eq. (20), the following system of linear algebraic J equations in terms of u^ is obtained. M X
J KIJ u^ ¼ f I ;
I ¼ 1; 2; . . . ; M:
ð21Þ
J¼1
where M is the total number of nodes. The so-called ‘‘stiffness’’ matrix K and ‘‘load’’ vector f are defined as Z Z Z BI DBJ dO þ a WI SUJ dG WI SNDBJ dG ð22aÞ KIJ ¼ OIs
fI ¼
Z
GIsu
Z
GIst
WI t dG þ a
GIsu
GIsu
WI u dG þ
Z OIs
WI b dO
where nx, ny and nz represent the component of outward normal vector on the boundary of the local sub-domain. E and n are the Young’s modulus and the Poisson’s ratio, respectively. The Young’s modulus is assumed to vary exponentially through the thickness according to EðzÞ ¼ Eh egðz=h1Þ ;
fJ;x
6 6 0 6 6 6 0 6 J B ¼6 6 fJ 6 ;y 6 6 6 0 4
0
fJ;x fJ;z
fJ;z
2
7 2 I 0 7 7 7 c;x fJ;z 7 7 I 6 7; B ¼ 6 6 0 4 0 7 7 7 0 fJ;y 7 7 5
0
g ¼ lnðE0 =Eh Þ
ð25Þ
where E0 = E(0) and Eh = E(h). Also, u, v and w are displacement components of a point along the x, y and z axes, respectively and cI is the test function for the MLPG1, which was defined in Eq. (17). By using relations (10) and (12), the unknown variables and their derivatives are obtained. Considering the MLPG5 method, the stiffness matrix and the load vector are given by Z KIJ ¼
ð22bÞ fI ¼
Z GIst
GIsi
Z NDBJ dG
Z t dG þ a
GIsu
GIsu
Z SNDBJ dG þ a
u dG þ
Z OIs
GIsu
SUJ dG
ð26aÞ
b dO
ð26bÞ
O 0
0
cI;y
cI;y
0
cI;x
cI;z
0
cI;z
0
cI;y
0
cI;z
3
nx
0
0
ny
0
6 N¼4 0 0
ny
0
nx
nz
0
nz
0
ny
nz
3
a
b
cI;x
z
y Side 4
2
fJ
6 J 07 5;U ¼ 6 4 0 nx 0
0
fJ 0
0
3
2
Side 2
Side 1
7 0 7 7; 5
fJ;x
0
ð24Þ
3
0
fJ;y
0 r z rh
where g is the inhomogeneity parameter. Poisson’s ratio is considered to be constant. Since a FGM is assumed to be isotropic and the Poisson’s ratio to be constant, it is possible to express the inhomogeneity parameter g in terms of the values of the Young’s modulus E= E(z) on the top (z= 0) and the bottom (z =h) surfaces of the plate as
where in 3D space, the matrices and the vectors in Eqs. (22) can be expressed as 2
ð23bÞ
cI
6 7 I 6 0 7 5; W ¼ 4 0 J f 0
0
cI 0
0
c
E0
3
7 0 7 5; I
x Side 3
h
z
Fig. 2. 3D coordinate system and problem dimensions.
Eh
ARTICLE IN PRESS 568
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
3. Defining the problem domain and boundary conditions A rectangular Cartesian coordinate system, shown in Fig. 2, is used to describe the geometry of a rectangular elastic plate with the volume O defined by 0oxoa, 0oyob and 0ozoh. The axes of x and y are positioned in the top surface while, z axis is directed along the thickness. Nodes are located in all three directions to discretize the plate continuum. By the present method, unequal distribution of nodes in different directions is possible. Therefore, nodes can be added in each direction easily to obtain desirable accuracy. The edges x=0,a and y=0,b are boundaries in which the essential boundary conditions or the natural boundary conditions are defined. The top and bottom surfaces of the plate (z=0,h) are loading surfaces. Four types of edge boundary conditions are considered, which are nominated as the simply supported (S), clamped (C), free (F) and guided (G) boundary conditions. These boundary conditions on the edges of the plate are defined as follows:
(ii) Clamped edge (C) u ¼ v ¼ w ¼ 0 on x ¼ 0ðor x ¼ aÞ; or y ¼ 0ðor y ¼ bÞ: (iii) Free edge (F)
sxx ¼ sxy ¼ sxz ¼ 0 on x ¼ 0ðor x ¼ aÞ; or syy ¼ syx ¼ syz ¼ 0 on y ¼ 0ðor y ¼ bÞ: sxy ¼ sxz ¼ 0; u ¼ 0; on x ¼ 0ðor x ¼ aÞ; or syx ¼ syz ¼ 0; v ¼ 0; on y ¼ 0ðor y ¼ bÞ:
1.04 1.03 1.02
on y ¼ 0ðor y ¼ bÞ:
ð27aÞ
1.03
w/w*
1.01 w/w*
1.04
1 0.99
1.02
0.98
1.01
0.97
1
0.96
0.99
0.95 γ =10-1 (MLPG1) γ =10-1 (MLPG5) γ =-10-1 (MLPG1) γ =-10-1 (MLPG5) γ =0 (MLPG1) γ =0 (MLPG5)
0.98 0.97 0.96 5x5x5
7x7x7
g
3D elasticity solution [10]
3D MLPG1
Error (%)
3D MLPG5
Error (%)
10 1 10 2 10 3 10 4 10 5 10 6 10 10 10 10 10 10 0
1.4146 1.3496 1.3433 -1.3426 1.3426 1.3426 1.2740 1.3355 1.3419 1.3425 1.3425 1.3425 1.3430
1.4096 1.3448 1.3384 1.3378 1.3377 1.3377 1.2693 1.3312 1.3370 1.3377 1.3377 1.3377 1.3377
0.35 0.36 0.36 0.36 0.36 0.36 0.37 0.32 0.37 0.36 0.36 0.36 0.39
1.4229 1.3573 1.3509 1.3503 1.3502 1.3502 1.2811 1.3432 1.3495 1.3502 1.3502 1.3502 1.3502
0.58 0.57 0.57 0.57 0.57 0.57 0.56 0.58 0.57 0.57 0.57 0.57 0.54
2 3 4 5 6
13x13x5
13x13x7 13x13x9 13x13x11 13x13x13 13x13x15 Nodal density
0
9x9x9 11x11x11 13x13x13 15x15x15 Nodal density
-1
Eh/E0:
-2
0.1 Analytical[10] 0.1 Present(MLPG1) 0.1 Present(MLPG5) 10 Analytical [10] 10 Present(MLPG1) 10 Present(MLPG5) 1 Analytical [10] 1 Present(MLPG1) 1 Present(MLPG5)
⎯w
Table 1 Non-dimensional out-of-plane displacement wða=2; b=2; h=2Þ in FG square plates with SSSS boundary conditions and SL loading condition. (a/h= 3).
γ= 10-1 (MLPG1) γ= 10-1 (MLPG5) γ= -10-1 (MLPG1) γ= -10-1 (MLPG5) γ= 0 (MLPG1) γ= 0 (MLPG5)
Fig. 4. Convergence of the results for the central deflection of an SSSS square plate (a/h= 3) for various nodal densities along the thickness direction.
Fig. 3. Convergence of the results for the central deflection of an SSSS square plate (a/h= 3) for various nodal densities.
1
ð27dÞ
For example, the SFCG plate is simply supported, free, clamped and guided along the sides 1, 2, 3 and 4, respectively (see Fig. 2).
sxx ¼ 0; v ¼ w ¼ 0 on x ¼ 0ðor x ¼ aÞ; or syy ¼ 0;
0.95
ð27cÞ
(iv) Guided edge (G)
(i) Simply supported edge (S)
u¼w¼0
ð27bÞ
-3
-4
-5
0
0.25
0.5
0.75
1
⎯z Fig. 5. Variation in non-dimensional deflection wð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
The components of the stress tensor are related to components of strain tensor as follows: n o rT ¼ sxx syy szz syz szx sxy ¼ De ð28Þ In Eq. (28), e denotes the infinitesimal strain with the following components 8 9 > > fI;x u^ I 8 9 > > > > > > exx > > > > > > > > I I > > > > ^ > > > > f v > eyy > > > ;y > > > > > > > > > > > > I I > > > > N < ezz = X < = ^ f;z w e ¼ 2e ¼ ð29Þ I I I I > > > > ^ ^ yz f þ f w v > > ;y > > ;z > > > > > > I¼1> > > > > > 2ezx > > > I I I > > > > ^I> > > > > > > f;z u^ þ f;x w > > : > > > > 2exy ; > I > I I I > : f v^ þ f u^ > ; ;x
569
In the current study two different loading conditions are considered
Uniform load q distributed on the top surface (UL) szz ¼ q;
szx ¼ szy ¼ 0 on z ¼ 0szz ¼ 0; szx ¼ szy ¼ 0 on z ¼ h
ð30aÞ
Sinusoidal load with the intensity q distributed on the bottom surface (SL)
szz ¼ 0;
szx ¼ szy ¼ 0 on z ¼ 0 px py sin ; szx ¼ szy ¼ 0 on z ¼ h
szz ¼ q sin
a
b
ð30bÞ
;y
-0.9 4
-0.8
3
-0.7
2 -0.6
⎯σxz
⎯σxx
1 0
Eh/E0: -1
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-2 -3 -4
0.25
0.5
Eh/E0:
-0.4
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-0.3 -0.2 -0.1 0
-5 0
-0.5
0.75
0
0.25
1
⎯z Fig. 6. Variation in non-dimensional in-plane normal stress s xx ð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
0.5
0.75
1
⎯z Fig. 8. Variation in non-dimensional transverse shear stress s xz ð0; 0:5b; zÞ versus z of an SSSS square plate.
2
0
1
-0.25
⎯σzz
⎯σxy
0 Eh/E0: -1
0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-2
-3 0
0.25
0.5
0.75
-0.5 Eh/E0: 0.1 Analytical [10] 0.1 Present (MLPG1) 0.1 Present (MLPG5) 10 Analytical [10] 10 Present (MLPG1) 10 Present (MLPG5) 1 Analytical [10] 1 Present (MLPG1) 1 Present (MLPG5)
-0.75
-1
1
⎯z Fig. 7. Variation in non-dimensional in-plane shear stress s xy ð0; 0; zÞ versus z of an SSSS square plate.
-1.25
0
0.25
0.5
0.75
1
⎯z Fig. 9. Variation in non-dimensional out-of-plane normal stress s zz ð0:5a; 0:5b; zÞ versus z of an SSSS square plate.
ARTICLE IN PRESS 570
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
4. Numerical results and discussion
4.1. Comparison studies
Several examples have been studied to examine the accuracy of the presented method in 3D analysis of thick FG rectangular plates with various boundary conditions. It is assumed that the plate is subjected to uniform or sinusoidal distributed loads (UL and SL loading conditions) and the body forces are assumed to be 0. Also, the Poisson’s ratio is considered to be 0.3. In the following analysis, dimensionless size of the sub-domain and support domain are considered as as = 0.75 and ai = 2.6. However, the effect of the value of the parameters as and ai on the accuracy of the solution is studied by Vaghefi et al. [42]. It should be noted that the number of gauss points is considered to be 216 in the local sub-domain integration. To present the obtained results the following non-dimensional parameters are introduced.
In order to examine the efficiency and accuracy of the presented approach, the obtained results are compared with the available exact 3D solution given by Kashtalyan [10]. First, the effects of nodal density on the obtained results for an FG square plate with simply supported boundary conditions (SSSS) are examined and presented in Fig. 3. Various nodal densities containing 5 5 5, 7 7 7, y, 15 15 15 regularly placed along the x-, y- and z-directions have been considered. In Fig. 3, w* refers to the exact solution obtained by Kashtalyan [10] for a SSSS FG square plate with a/h =3. As shown in the figure, it is observed that by choosing 13 13 13 nodes (2197 nodes), the solutions are converged to the exact ones for all the cases. Also, the effect of nodal density on the plate thickness direction, the direction in which the material properties change, is examined. The results for a square FG plate with SSSS boundary conditions with a/h= 3 are depicted in Fig. 4. In this figure, various nodal densities containing 13 13 5, 13 13 7, y, 13 13 15 nodes have been considered. According to the figure it can be observed that considering 13 and more nodes along the
w¼
sij Gh w z ; z¼ ; s ij ¼ q hq h
ð31Þ
In Eq. (31), the parameter Gh is the shear modulus on the edge z= h of the plate.
Table 2 Numerical results for FG square plates with various boundary conditions and UL loading condition (h/a= 0.2). Eh/E0
s xx (a/2, b/2, 0)
w max MLPG1
MLPG5
FEM
s yy (a/2, b/2, 0)
MLPG1
MLPG5
FEM
MLPG1
MLPG5
FEM
(a) SSSS 0.1 0.5 1 2 10
4.0916 8.9751 12.5997 17.6640 39.0605
4.1596 8.9357 12.6375 17.8397 39.0385
4.1215 9.0047 12.6134 17.7118 39.1558
15.3564 9.2902 7.4462 5.9410 3.4665
15.4390 9.3279 7.4199 5.9711 3.4944
15.403 9.2995 7.4588 5.9591 3.4805
15.3564 9.2902 7.4462 5.9410 3.4665
15.4390 9.3279 7.4199 5.9711 3.4944
15.4032 9.2995 7.4588 5.9591 3.4805
(b) SCSC 0.1 0.5 1 2 10
2.4311 5.4661 7.7215 10.8228 23.5185
2.4270 5.5147 7.7419 10.8902 23.6686
2.4373 5.4851 7.7336 10.8517 23.6058
9.1729 5.6761 4.6394 3.7090 2.2461
9.0727 5.7029 4.6203 3.7128 2.2267
9.1343 5.6959 4.6342 3.7334 2.2364
11.1523 6.7495 5.4580 4.3708 2.5759
11.2305 6.7348 5.4618 4.3291 2.5517
11.1814 6.7647 5.4505 4.3464 2.5650
(c) SSSF 0.1 0.5 1 2 10
12.4349 27.2606 38.4889 54.3348 122.7609
12.4062 27.5442 38.3326 54.9850 124.3991
12.5132 27.4253 38.6582 54.6430 123.4156
25.1382 15.1983 12.2387 9.6918 5.5863
25.4982 15.4213 12.3027 9.6697 5.6536
25.3471 15.2947 12.2592 9.7428 5.6276
11.9221 7.1349 5.7356 4.5566 2.7440
11.7441 7.0864 5.6989 4.5886 2.7426
11.8456 7.1255 5.7244 4.5778 2.7295
(d) SSSG 0.1 0.5 1 2 10
9.5911 20.9661 29.4798 41.8456 94.0184
9.7573 21.1375 29.4401 41.5796 92.9495
9.6649 21.0676 29.5869 41.7467 93.5501
25.3203 15.4382 12.3969 9.7842 5.7340
25.7049 15.2773 12.2797 9.7465 5.6173
25.5384 15.3930 12.3336 9.8204 5.6623
14.7404 8.8260 7.1261 5.7314 3.3589
14.5421 8.8355 7.0901 5.8132 3.3819
14.6202 8.8759 7.1448 5.7790 3.3731
(e) CFCF 0.1 0.5 1 2 10
3.5627 8.1319 11.6473 16.3762 36.9346
3.5879 8.1895 11.5694 16.3885 36.4419
3.5720 8.1516 11.6298 16.4360 36.7001
13.5523 8.1803 6.5891 5.2458 3.0279
13.3829 8.1169 6.5536 5.1909 3.0725
13.4943 8.1593 6.5822 5.2223 3.0486
3.3407 2.0525 1.6935 1.3969 0.94076
3.2949 2.0471 1.6802 1.3825 0.94204
3.3156 2.0398 1.6874 1.3884 0.93952
(f) CGCG 0.1 0.5 1 2 10
3.2034 7.2736 10.2880 14.6290 31.5918
3.2327 7.3329 10.3664 14.3906 31.8806
3.2145 7.2853 10.3357 14.5306 31.9116
13.9612 8.3649 6.8187 5.4403 3.1722
13.9381 8.4862 6.7709 5.3926 3.1434
13.8783 8.4172 6.8003 5.4060 3.1662
4.4350 2.8342 2.3295 1.9301 1.2416
4.4308 2.8105 2.3406 1.9084 1.2555
4.4510 2.8174 2.3374 1.9167 1.2468
(g) CCCC 0.1 0.5 1 2 10
1.7485 3.9547 5.6094 7.8809 16.9989
1.7457 3.9251 5.6221 7.8188 17.1946
1.7340 3.9499 5.5903 7.8589 17.1054
8.0702 5.0016 4.0694 3.2773 1.9772
8.1498 4.9610 4.0265 3.2334 1.9517
8.1096 4.9837 4.0535 3.2548 1.9706
8.0702 5.0016 4.0694 3.2773 1.9772
8.1498 4.9610 4.0265 3.2334 1.9517
8.1096 4.9837 4.0535 3.2548 1.9706
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
571
plate thickness direction can led to an acceptable solution. Therefore, to obtain the solutions in the following problems, the nodal distribution is considered to be 13 13 13. In Table 1, the non-dimensional out-of-plane displacement w of an SSSS FG square plate is given. Assuming a/h= 3 and a sinusoidal loading (see Eq. (30b)), various inhomogeneity parameters g are considered and the results are calculated in the plate center (x =0.5a, y=0.5b, z =0.5h). It should be mentioned that the load is assumed to be applied on the bottom (z= h) surface of the plate. It is clear while considering g = 0, the obtained solutions can be related to the same homogeneous plate. The results obtained by using both the methods, MLPG1 and MLPG5, are presented in the table. In Table 1 the percent of error is calculated using the following equation Present solutionExact solution ð32Þ 100 %Error ¼ Exact solution
Also, the variations in the non-dimensional stresses s xx , s xy , s xz and s zz are depicted through the thickness of the plate in Figs. 6–9, respectively. Different values of Eh/E0 are considered in plotting the figures and both the methods MLPG1 and MLPG5 are applied to obtain the results, which are compared with the reference solution [10]. In plotting Figs. 6–9, the plate is assumed to be simply supported with a/h= 3, under a sinusoidal loading at z= h. Numerical results in Figs. 6–9 are given for three different values of Eh/E0. The chosen values of the inhomogeneity ratio Eh/E0 = 10, 1.0, and 0.1 do not necessarily represent a certain material, but are rather used to show the effect of inhomogeneity on stress and displacement fields. It can be seen that the present results agree well with those reported by Kashtalyan [10].
where the exact solution refers to the analytic solution w* for an SSSS thick FG square plate with a/h= 3 given in Ref. [10]. In Fig. 5, variation in the non-dimensional out-of-plane displacement w is plotted through the thickness of the plate.
The FG plates with various boundary conditions are analyzed using the MLPG1 and MLPG5 methods. In Table 2, the maximum non-dimensional deflection w max and the non-dimensional normal stresses s xx and s yy at the center of the top surface
4.2. Numerical examples with various boundary conditions
Table 3 Numerical results for FG square plates with various boundary conditions and UL loading condition. (h/a =0.3). Eh/E0
s xx (a/2, b/2, 0)
w max MLPG1
MLPG5
FEM
MLPG1
s yy (a/2, b/2, 0) MLPG5
FEM
MLPG1
MLPG5
FEM
(a) SSSS 0.1 0.5 1 2 10
0.9707 2.1378 2.9853 4.1208 8.7134
0.9688 2.1498 2.9603 4.1098 8.6923
0.9732 2.1407 2.9792 4.1333 8.7293
7.2230 4.3084 3.4496 2.7499 1.6449
7.2034 4.2943 3.4959 2.7556 1.6566
7.2639 4.3378 3.4681 2.7673 1.6499
7.2230 4.3084 3.4496 2.7499 1.6449
7.2034 4.2943 3.4959 2.7556 1.6566
7.2639 4.3378 3.4681 2.7673 1.6499
(b) SCSC 0.1 0.5 1 2 10
0.6637 1.5302 2.1434 2.9696 6.3437
0.6675 1.5286 2.1555 3.0066 6.3506
0.6655 1.5269 2.1481 2.9890 6.3713
4.8701 3.0584 2.4868 2.0094 1.2493
4.8784 3.0311 2.4692 2.0085 -1.2412
4.8640 3.0484 2.4794 2.0125 1.2466
5.5628 3.3449 2.6511 2.1526 1.3128
5.5672 3.3437 2.6555 2.1543 1.3103
5.5558 3.3347 2.6608 2.1435 1.3073
(c) SSSF 0.1 0.5 1 2 10
2.7629 6.0992 8.5891 12.0999 26.9454
2.7597 6.1412 8.5866 12.1041 26.7625
2.7648 6.1154 8.6239 12.1192 26.8886
11.7326 7.0170 5.6075 4.4438 2.5999
11.7563 7.0088 5.6392 4.4358 2.5618
11.8060 7.0433 5.6250 4.4592 2.5836
5.3477 3.1655 2.5304 2.0512 1.2726
5.3210 -3.1502 2.5256 2.0381 1.2649
5.3603 3.1771 2.5496 2.0470 1.2677
(d) SSSG 0.1 0.5 1 2 10
2.1435 4.6698 6.5596 9.1481 20.0863
2.1407 4.6802 6.5399 9.1548 19.9522
2.1500 4.7085 6.5882 9.2033 20.0389
11.7572 6.9951 5.5844 4.4223 2.5761
11.7181 7.0800 5.6418 4.4690 2.5881
11.8310 7.0385 5.6143 4.4504 2.5790
6.8423 4.1385 3.3613 2.6698 1.6293
6.8262 4.1289 3.3419 2.6712 1.6251
6.8781 4.1569 3.3496 2.6892 1.6316
(e) CFCF 0.1 0.5 1 2 10
0.9624 2.2818 3.2947 4.7104 10.6373
0.9765 2.2871 3.3346 4.7115 10.6248
0.9686 2.2986 3.3229 4.7246 10.6851
6.6913 4.0416 3.2433 2.5593 1.5300
6.7873 4.0746 3.2478 2.5907 1.5260
6.7453 4.0482 3.2368 2.5790 1.5344
1.5581 0.94903 0.78921 0.68142 0.52378
1.5632 0.95112 0.78725 0.67291 0.52526
1.5515 0.95299 0.79313 0.67653 0.52173
(f) CGCG 0.1 0.5 1 2 10
0.8824 2.0743 2.9257 4.1037 8.8568
0.8798 2.0583 2.9246 4.1344 8.8598
0.8845 2.0650 2.9422 4.1138 8.8953
6.8910 4.1806 3.3383 2.6878 1.5976
6.9751 4.2011 3.3311 2.6746 1.5993
6.9342 4.1858 3.3568 2.6842 1.6062
2.3574 1.5543 1.2956 1.0974 0.77678
2.3489 1.5453 1.3128 1.1092 0.78193
2.3724 1.5510 1.3056 1.1022 0.78012
(g) CCCC 0.1 0.5 1 2 10
0.4948 1.1649 1.6554 2.3312 5.0918
0.5016 1.1753 1.6521 2.3469 5.1350
0.4985 1.1702 1.6626 2.3368 5.1065
4.2097 2.6062 2.0965 1.7194 1.0917
4.2029 2.5950 2.0944 1.7318 1.1040
4.2293 2.6120 2.1097 1.7277 1.0974
4.2097 2.6062 2.0965 1.7194 1.0917
4.2029 2.5950 2.0944 1.7318 1.1040
4.2293 2.6120 2.1097 1.7277 1.0974
ARTICLE IN PRESS 572
R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
Fig. 10, it can be seen that the through-thickness variation in the non-dimensional out-of-plane displacement w in the FG plates is not constant and it has the maximum value near the top surface for all values of stiffness gradient Eh/E0. It should be mentioned that with increase in the ratio Eh/E0, the out-of-plane displacement w of the FG plates decreases. Also in Fig. 11, variation in the in-plane normal stress s xx through the thickness of a clamped FG plate under sinusoidal loading is shown. The plate has the same thickness ratio as in Fig. 10. According to Fig. 11, it can be seen that in FG plates, the maximum of the in-plane normal stress s xx is located near the surface with lower Young’s modulus, i.e. the maximum compressive stress happened near the top surface of the FG plate with Eh/E0 = 0.1 and the maximum of the tensile stress happened near the bottom surface of the FG plate with Eh/E0 =10.
0
⎯w
-1
Eh/E0: 0.1 0.5 1 2 10
-2
5. Conclusions
-3
0
0.25
0.5
0.75
1
⎯z Fig. 10. Variation in non-dimensional deflection wð0:5a; 0:5b; zÞ versus z of a CCCC square plate.
3
2
⎯σxx
1
0
-1
Eh/E0: 0.1 0.5 1 2 10
-2
-3 0
0.25
A 3D and completely meshless approach have been developed for the static analysis of thick FG plates with various boundary conditions. The 3D equilibrium equations have been considered and the Young’s modulus of the plate was assumed to vary exponentially through the thickness, while the other material properties were assumed to be constant. The two meshless methods including MLPG1 and MLPG5 have been extended to solve the bending problem of FG plates. Brick-shaped local domains were considered in the solution procedure. The application of brick-shaped local domains facilitates the integration and the approximation procedure. The proposed approach for construction of the shape and the test functions makes it possible to introduce more nodes in the direction of material variation. Consequently, more precise solutions can be obtained easily and efficiently. It is observed that the approach works well in analysis of FG plates. Several numerical examples have been presented. The obtained results have been compared with the available analytical and numerical solutions in the literature and an excellent agreement was seen. Although the approach has been applied for static analysis, it can be used for the dynamics, buckling and fracture analyses of thick FG plates as well.
References
0.5 ⎯z
0.75
1
Fig. 11. Variation in non-dimensional in-plane normal stress s xx ð0:5a; 0:5b; zÞ versus z of a CCCC square plate.
(z= 0) are given for square FG plates with various boundary conditions including SSSS, SCSC, SSSF, SSSG, CFCF, CGCG and CCCC. For the results presented in Table 2, the UL loading condition is assumed and the plates with the thickness ratios of h/a =0.2 and Eh/E0 = 0.1, 0.5, 1, 2, and 10 are considered. To perform a better comparison, the FEM results obtained by the authors are also presented. The FEM results are obtained by using the ANSYS software with 36 36 36 solid 45 brick elements. Table 3 presents the same results as Table 2 while considering the thickness ratio of h/a= 0.3. In Fig. 10, the through-thickness variation in the nondimensional out-of-plane displacement w of a CCCC FG plate under sinusoidal loading is shown. The results are obtained considering a/h=3 and different values of Eh/E0. According to
[1] Igari T, Notomi A, Tsunoda H, Hida K, Kotoh T, Kunishima S. Material properties of functionally gradient material for fast breeder reactor. In: Yamanouochi M, Koizumi M, Hirai T, Shiota I, editors. Proceedings of the first international symposium on functionally gradient materials, Japan: Sendai; 1990. p. 209–214. [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 Appl Biomater 1994;5(1):23–7. [3] Tani J, Liu GR. Surface waves in functionally gradient piezoelectric plates. JSME Int J Ser A (Mech Mater Eng) 1993;36(2):152–5. [4] Osaka T, Matsubara H, Homma T, Mitamura S, Noda K. Microstructural study of electroless-plated CoNiReP/NiMoP doublelayered media for perpendicular magnetic recording. Jpn J Appl Phys 1990;29(10):1939–43. [5] Watanabe Y, Nakamura Y, Fukui Y, Nakanishi K. Amagnetic-functionally graded material manufactured with deformation induced martensitic transformation. J Mater Sci Lett 1993;12(5):326–8. [6] Getto H, Ishihara S. Development of the fire retardant door with functionally gradient wood. In: Proceedings of the fourth international symposium on functionally graded materials, Japan: Tsukuba City; 1996. [7] Levinson M. The simply supported rectangular plate: an exact, three dimensional, linear elasticity solution. J Elasticity 1985;15:283–91. [8] Batra RC, Vel SS. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2001;40(7):1421–33. [9] Zhong Z, Shang ET. Three-dimensional exact analysis of a simply supported functionally gradient piezoelectric plate. Int. J. Solids Struct 2003;40(20):5335–52. [10] Kashtalyan M. Three-dimensional elasticity solution for bending of functionally graded rectangular plates. Eur J Mech A—Solid 2004;23:853–64.
ARTICLE IN PRESS R. Vaghefi et al. / Engineering Analysis with Boundary Elements 34 (2010) 564–573
[11] Liew KM, Teo TM. Modeling via differential quadrature method: threedimensional solutions for rectangular plates. Comput Methods Appl Mech Eng 1998;159:369–81. [12] Liew KM, Teo TM, Han JB. Three-dimensional solutions of rectangular plates by variant differential quadrature method. Int J Mech Sci 2001;43:1611–28. ¨ . Three-dimensional vibration, buckling and bending analysis of [13] Civalek O thick rectangular plates based on discrete singular convolution method. Int J Mech Sci 2007;49:752–65. [14] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Eng 1994;37:229–56. [15] Liu WK, Jun S, Zhang Y. Reproducing kernel particle methods. Int J Numer Methods Fluids 1995;20:1081–106. [16] Duarte CA, Oden JT. Hp-cloud-a meshless method to solve boundary-value problems. Comput Methods Appl Mech Eng 1996;139:237–62. [17] Babuska I, Melenk J. The partition of unity method. Int J Numer Methods Eng 1997;40:727–58. [18] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [19] Atluri SN, Kim HG, Cho JY. A critical assessment of the truly meshless local Petrov–Galerkin (MLPG), and local boundary integral equation (LBIE) methods. Comput Mech 1999;24:348–72. [20] Atluri SN, Shen S. The meshless local Petrov–Galerkin (MLPG) method: a simple and less-costly alternative to the finite element and boundary element methods. CMES: Comput Modeling Eng Sci 2002;3(1):11–51. [21] Atluri SN, Shen S. The meshless LOCAL Petrov–Galerkin (MLPG) method. USA: Tech Science Press; 2002. [22] Atluri SN, Zhu T. The meshless local Petrov–Galerkin (MLPG) approach for solving problems in elasto-statics. Comput Mech 2000;25:169–79. [23] Batra RC, Ching HK. Analysis of elastodynamic deformation near a crack/ notch tip by the meshless local Petrov–Galerkin (MLPG) Method. CMES: Comput Modeling Eng Sci 2002;3(6):717–30. [24] Lin H, Atluri SN. The meshless local Petrov–Galerkin (MLPG) method for solving incompressible Navier–Stokes equations. CMES: Comput Modeling Eng Sci 2001;2(2):117–42. [25] Lin H, Atluri SN. Meshless local Petrov–Galerkin (MLPG) method for convection-diffusion problems. CMES: Comput Modeling Eng Sci 2000;1(2):45–60. [26] Sladek J, Sladek V, Atluri SN. A pure contour formulation for the meshless local boundary intergral equation method in thermoelasticity. CMES: Comput Modeling Eng Sci 2001;2(4):423–33. [27] Atluri SN, Cho JY, Kim HG. Analysis of thin beams, using the meshless local Petrov–Galerkin method, with generalized moving least squares interpolations. Comput Mech 1999;24:334–47. [28] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) formulation for static and free vibration analyses of thin plates. CMES: Comput Modeling Eng Sci 2001;2(4):463–76.
573
[29] Long S, Atluri SN. A meshless local Petrov–Galerkin (MLPG) method for solving the bending problem of a thin plate. CMES: Comput Modeling Eng Sci 2002;3(1):53–64. [30] 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. [31] Xia P, Long SY, Cui HX, Li GY. The static and free vibration analysis of a nonhomogeneous moderately thick plate using the meshless local radial point interpolation method. Eng Anal Boundary Elem 2009;33(6):770–7. [32] Sladek J, Sladek V, Zhang Ch. Stress analysis in anisotropic functionally graded materials by the MLPG method. Eng Anal Boundary Elem 2005;29:597–609. [33] Wang H, Qin QH. Meshless approach for thermo-mechanical analysis of functionally graded materials. Eng Anal Boundary Elem 2008;32(9): 704–12. [34] Kim HG, Atluri SN. Arbitrary placement of secondary nodes, and error control, in the meshless local Petrov–Galerkin (MLPG) method. CMES: Comput Modeling Eng Sci 2000;1(3):11–32. [35] Ching HK, Batra RC. Determination of crack tip fields in linear elastostatics by the meshless local Petrov–Galerkin (MLPG) method. CMES: Comput Model Eng Sci 2001;2(2):273–90. [36] Koohkan H, Baradaran GH, Vaghefi R. A completely meshless analysis of cracks in isotropic functionally graded materials. In: Proceedings of the IMechE, Part C: J Mech Eng Sci 2009, in press. [37] Tang Z, Shen S, Atluri SN. Analysis of materials with strain-gradient effects: A meshless local Petrov–Galerkin (MLPG) approach, with nodal displacements only. CMES: Comput Modeling Eng Sci 2003;4(1): 177–96. [38] Li Q, Shen S, Han ZD, Atluri SN. Application of meshless local Petrov–Galerkin (MLPG) to problems with singularities, and material discontinuities, in 3-D Elasticity. CMES: Comput Modeling Eng Sci 2003;4(5):571–85. [39] Han ZD, Atluri SN. Truly meshless local Petrov–Galerkin (MLPG) solutions of traction and displacement BIEs. CMES: Comput Modeling Eng Sci 2003;4(6):665–78. [40] Han ZD, Atluri SN. Meshless local Petrov–Galerkin (MLPG) approaches for solving 3-D problems in elasto-statics. CMES: Comput Modeling Eng & Sci 2004;6(2):169–88. [41] Han ZD, Atluri SN. A meshless local Petrov–Galerkin (MLPG) approach for 3dimensional elasto-dynamics. CMC: Comput Mater Continua 2004;1(2): 129–40. [42] Vaghefi R., Baradaran G.H., Koohkan H., Three-dimensional static analysis of rectangular thick plates by using meshless local Petrov–Galerkin method. In: Proceedings of the IMechE, Part C: J. Mech Eng Sci, vol 223, C9; 2009. p. 1983–1996