Boundary element method applied to the bending analysis of thin functionally graded plates

The present work introduces the boundary element method applied to the bending analysis of functionally graded plates. It is assumed that material properties are graded through the thickness direction of the plate according to a power law distribution. The neutral surface position for such plate is determined and the classical plate theory based on the exact neutral surface position is employed to extract the equilibrium equations. A direct approach based on the Green’s identity is used to formulate boundary element method. By introducing a novel approach, domain integrals which arise from distributed transverse loads are transformed into boundary integrals. In case studies, three geometrical shapes including, rectangular, circular and elliptic for functionally graded plates with/without hole are considered. Comparative studies are first carried out to evaluate the sufficiency of the proposed method for bending analysis of isotropic and functionally graded plates subjected to the transverse loads. Then, a series parametric study is performed to examine the influences of the power of functionally graded material, boundary conditions and geometrical parameters on the deformation and stress of functionally graded plates.


INTRODUCTION
Thin plates are light weight structures with high load-carrying capacity, and technological effectiveness.These are extensively used in engineering applications, such as aerospace, mechanical and civil and ocean engineering which manufactured in various geometrical shapes and made from various materials.In recent years, a new class of composite materials known as functionally graded materials (FGMs) has gained considerable attention as advanced structural materials.FGMs have superior thermal and mechanical performance to conventional homogeneous materials.Due to this fact, they have a wide variety of engineering applications especially for the purpose of removing mismatches of thermo-mechanical properties between two neighboring layers and reducing stress level in structures.A number of research works have been carried out in static and dynamic analysis of functionally graded (FG) plates (see e.g.Ref. [1,5,11,[19][20][21]27]).For instance, in the field of static analysis, Reddy et al [19] investigated the axisymmetric bending of functionally graded solid and annular circular plates.They presented exact solutions for deflections, force and moment resultants using the first-order shear deformation theory (FSDT).Reddy and Cheng [20] studied the three-dimensional thermo mechanical deformations of simply supported, functionally graded rectangular plates by using an asymptotic method.In general, the material properties of FG plate do not have symmetry about the middle plane of the FG plate.This condition leads to the stretchingbending coupling in constitutive equations of FG plate.Abrate [1] showed that there is no stretching-bending coupling in constitutive equations of FG plates, if the reference surface is properly selected.Zhang and Zhou [27] presented a theoretical analysis to the FG thin plates based on the physical neutral surface.They carried out the bending, buckling and free vibration analysis of simply supported rectangular and clamped circular FG plates.Recently, Bodaghi and Saidi [5] developed an exact analytical solution for buckling of functionally graded rectangular plates subjected to nonuniformly distributed in-plane loading acting on two opposite simply supported edges.
As mentioned above, analytical solutions of plate problems using classical methods is limited to relatively simple plate geometry, load configuration, and boundary supports.If these conditions are more complicated, the classical analysis methods become increasingly tedious or even impossible.In such cases, approximate methods are the only approaches that can be employed for the solution of practically important plate problems.The boundary element method (BEM) is a popular computational tool among numerical methods.The BEM is well suited for treating complicated boundaries, for discontinuous internal actions, mixed boundary conditions, etc.In particular, the BEM has been applied successfully for the solutions of plate bending problems, its advantages have been demonstrated [3,4,6,8,12,15,22,23,26].One of the main advantages of this method is a possibility of reduced dimensionality of the problem which leads to a reduced set of equations and makes smaller amount of data required for computation.For the linear problem, the BEM is especially effective because integral equations are formulated only on the boundary of the domain under consideration.The various boundary element formulations for thin plate bending can be generally categorized as the direct and the indirect formulations.The first boundary element formulation was based on the so-called indirect methods [12], where the integral equations do not relate to the natural variables of the problem, such as deflections, rotations, bending or twisting moments and shear forces.Instead, they involve some source distribution densities, apparently without any physical meaning.Some authors have followed the same bases to propose alternative indirect schemes [3,26].The natural variables appearing in the integral equations gave origin to the direct methods applied to plate bending, as in Bezine [4], Stern [22] and Tottenham [23].They established the main bases of the boundary element technique for plate bending, nowadays used as a standard tool.
In the general plate bending boundary element method, domain integrals arise in the formulation owing to the distributed load on the domain.The fact that domain integrals need to be evaluated spoils the pure boundary character of the BEM and weakens the advantages this method has over domain methods.During the past three decades, various techniques such as Galerkin tensor method [6], dual reciprocity method [15], multiple reciprocity method [16] and radial integration method [8] have been developed that successfully overcome this problem and at the same time preserve the Latin American Journal of Solids and Structures 10(2013) 549 -570 purely boundary character of the BEM.A brief comment on the popularly used transformation methods was described by Gao [8].Using the direct boundary element method based on Green's identity, Gospodinov and Ljutskanov [10] investigated the static analysis of thin rectangular plates.They discretized the plate domain to solve domain integrals in their BEM formulation.They also used the indirect boundary element method for dynamic and stability analyses of the plates.Du et al. [7] presented some fundamental aspects of the direct boundary element method of the Kirchhoff theory of thin plate.The boundary element results of plates in bending can exhibit disturbances in the vicinity of points where abrupt changes in the boundary conditions occur.Due to this event, Venturini and Paiva [24] proposed several different ways of defining the boundary element system of equations to improve their numerical responses and consequently to increase the reliability of the technique.Paris and de Leon [17] formulated the thin plate bending problem by means of two coupled Poisson equations.The domain integrals were evaluated approximating the integrands by a series of simple domain functions whose coefficients were calculated by a collocation procedure at points placed along the boundary and domain.Recently, Leonetti et al. [14] presented a symmetric model for the boundary element analysis of Kirchhoff plates.They showed that the convergence of their model is more regular than the collocation boundary model.
Although some research works dealt with bending analysis of thin isotropic plates, but there are few research works that dedicate to the bending analysis of non-homogeneous plate.For example, Paiva et al. [18] presented an approach for anisotropic thin-plate bending problems using the boundary element formulation when the source points were located on the boundary.Albuquerque et al. [2] presented a boundary element formulation without any domain integral for anisotropic plate bending problems using the radial integration method.A robust boundary element method that can be used to solve elastic problems with nonlinearly varying material parameters, such as the functionally graded material and damage mechanics problems was presented by Gao et al. [9].Recently, Ruocco and Minutolo [21] investigated two-dimensional stress analysis of multi-region functionally graded materials using a field boundary element model.To the best of authors' knowledge, there is no research work on the bending analysis of functionally graded plate based on the boundary element method.
The present work develops a pure boundary element method for bending analysis of functionally graded plates based on the classical plate theory (CPT) and physical neutral surface concept.The material properties of the FG plate are assumed to vary continuously and smoothly through the thickness according to the power-law distribution of the volume fraction of the constituents.The equilibrium equations are derived from the principle of minimum total potential energy.The direct boundary element method is employed to solve plate bending problem.By introducing a novel approach and using auxiliary potential functions, domain integrals which arise from distributed loads are transformed into boundary integrals.To certify the accuracy of the present boundary element method, the results obtained by the present analysis are compared with those available in the literature.Moreover, the effects of power of FGM and geometrical parameters together with various combinations of boundary conditions on the deformation and stress of FG plates are investigated in detail.

Physical neutral surface concept
The material properties of FGMs vary smoothly and continuously from one surface to the other.This is achieved by gradually varying the volume fraction of the constituent materials.FGMs are usually made from a mixture of metal and ceramic, or a combination of different metals.In the present work, it is assumed that the FG plate is made of ceramic and metal.Since, the material properties of the FG plate vary through the thickness direction, the neutral surface of such plate may not coincide with its geometric middle surface.Therefore, stretching and bending deformations of FG plates are coupled.Some researchers [1,5,27] have shown that there is no stretchingbending coupling in constitutive equations if the origin of coordinate system is suitably selected in the thickness direction of the FG plate so as to be the neutral surface.To specify the position of neutral surface of FG plates, two different planes are considered for the measurement of z , namely, ms z and ns z , as depicted in Fig. 1.The volume-fraction of metal (V m ) can be expressed based on these coordinates as [5,27]: where h is the thickness of the plate, the parameter D is the distance of neutral surface from the middle surface and n denotes the power of FGM ( 0) n ≥ .The effective Young's modulus (E)   based on the Voigt model can be expressed as: where E m and E c are the Young's modulus of the metal and ceramic, respectively.The situation of the neutral surface of the FG plate is determined to satisfy the first moment with respect to Young's modulus being zero as follows [5,27]: Latin American Journal of Solids and Structures 10(2013) 549 -570 Therefore, the position of neutral surface can be obtained as: Eq. ( 4) shows that the distance of neutral surface from the middle surface (Δ) is zero for ho- mogeneous isotropic plates.

Governing equations
According to the classical plate theory and physical neutral surface concept, the displacement components of a material point within the plate domain in Cartesian coordinates system can be written as: where ψ 1 , ψ 2 and w are the displacements of neutral surface of FG plate along the x,y and z ns coordinate directions, respectively.Substituting Eqs. ( 5) into linear strain-displacement relationship, kinematic equations are obtained as: where The constitutive relations of the plate in the plane stress state are expressed as [25]:

) where
Latin American Journal of Solids and Structures 10(2013) 549 -570 The parameter n is Poisson ratio and it is assumed to be constant through the thickness of the FG plate [1,5,19,20,27].By using the principle of minimum total potential energy, the equilibrium equations of plate are derived as follows: Substituting Eq. ( 6) into Eq.( 9) and the subsequent results into Eqs.(11), gives the stress resultants in term of the displacements as: N = Aε 0 ; M = Dε 1 (12) where and The matrices A and D are extensional and bending stiffness matrices.Eq. ( 12) reveals the fact that there is no stretching-bending coupling matrix, and subsequently, the coupling in the constitutive equations when the physical neutral surface concept is utilized.Finally, the bending equation of the FG plate may be obtained upon substitution of Eqs.(12) into Eq.( 10) as below: where ∇ 2 is known as the Laplace operator, and the constant D is the equivalent flexural rigidity of the FG plate which is simplified as and D c = for a homogeneous fully metallic and ceramic plates, respectively.Utilizing the boundary element method, the bending equation ( 15) will be solved for stress analysis of FG plates in the next section.

Fundamental solution
Considering Eq. ( 15), the boundary element formulation is now well established [13].Let us consider a point source placed at point P(x,y) of the xy -plane as shown in Fig. 2a.Its density at Q(ξ, η) can be expressed by Dirac delta function as follows: The deflection v(Q,P) produced at point Q satisfies the governing equation as follows: A singular particular solution of Eq. ( 17) is called the fundamental solution of the governing equation (15).Applying Green's second identity for w and ∇ 2 v yields: where n is normal to the boundary Γ .Using Eq. ( 18) and following Ref.[17], the fundamen- tal solution can be obtained as v = 1 8π r 2 ln r .
The governing equation ( 15) can be written in the new form by introducing a potential function as: where By considering Eqs. ( 17), ( 19) and ( 20), and also using Green's second identity for ϕ and ν , the representation of integral equation ( 18) can be rewritten as follows: where ε P is a coefficient which depends on the position of point P .This coefficient for the smooth boundary is defined as follows [13]: It should be noted that the integral representation of the solution (21) involves the following domain integral: Although the integrand v F D is known, the fact that domain integral need to be evaluated spoils the pure boundary character of the method.In the present work, to overcome this drawback, a novel procedure is introduced by converting the domain integral (23) to a boundary integral.To this end, an auxiliary function called v is assumed by satisfying the following equation: Then, Green's second identity for v and ∇ 2 v is written as: Now, by substituting Eq. ( 24) into Eq.( 25) and using of Green's second identity for ∇ 2 v and v , the following equation can be obtained: It is seen that the domain integral (23) has been successfully converted to the above boundary integral.Finally, substituting Eq. ( 26) into Eq.( 21) yields: In Eq. ( 27), the quantity ϕ and its derivation ϕ n are unknown.To specify these quantities, a new potential function called u and satisfied the following Laplace equation is assumed.
In the well-known procedure as expressed in Ref. [13], the fundamental solution u can be de- rived as u = 1 2π lnr .By virtue of Eqs. ( 19) and (28), Green's second identity for j and u can be written as: It can be seen that the integral representation of the solution (29) involves the domain integral.To convert this domain integral to a boundary integral, similar process that was presented for Eq. ( 23) were carried out.
An auxiliary function by satisfying the following Poisson equation is considered as: Applying Green's identity for the functions u and u that satisfies Eq. (30) yields: Substituting Eq. (31) into Eq.( 29) gives: Thus, the bending equation of FG plates ( 15) is reformulated into two dependent boundary integral solutions ( 27) and (32).

The BEM with superparametric boundary elements
In the present work, superparametric elements are employed to discretize the plate boundary [13].The boundary Γ is discretized into N superparametric elements (see Fig. 2b), which are num- bered in the counter-clockwise sense.The geometry is modeled by a parabolic arc (two degree polynomial), whereas the boundary quantity is assumed to be constant along the element and equal to its value at the nodal middle point.The discretized form of Eqs. ( 27) and (32) are expressed for a given point p i on Γ as: where Γ j is the segment on which the j th node is located and over which integration is carried out.Also, p i is the nodal point of the i th element which remains constant, while the point q varies over the j th element (see Fig. 2b).Definition of all parameters presented in Eq. ( 33) and (34) are expressed as follows: Latin American Journal of Solids and Structures 10(2013) 549 -570 ( ) ,n q = ( ) ,x q n x + ( ) ,y q n y , ( ) ,n iq = ( ) ,r iq iq ∂r iq ∂n = ( ) ,r iq iq cos α iq (35a) where α iq is the angle between the vector r iq and the unit vector n which is normal to the boundary at point q and the S, x and y are the shape function and elemental position vectors, respectively.Equations ( 33) and (34) are applied consecutively for all the nodes p i (i = 1...N ) yielding two systems of N linear algebraic equations, which are arranged in matrix form as follows: These two systems of equations are coupled and should be simultaneously solved.Prior to solving them, it is necessary to separate the unknown from the known quantities which is possible by applying the plate boundary conditions.

Application of boundary conditions
The functionally graded plate boundary can be any combination of clamped, simply supported and free boundary conditions.These boundary conditions which are developed from the principle of minimum total potential energy can be expressed in terms of normal and tangential coordinates as follows: Clamped: w = w ,n = 0 (38a) where M n and Q eff are the resultant normal moment and the effective shear force, respectively, which are defined as: where the Laplacian has the form: in which R is the radius of curvature of boundary.By using Eq. ( 20), these boundary conditions can be simplified for j th node as follows: Clamped: Simply supported: w j = 0 and The known boundary conditions are applied at nodal points, and then Eqs. ( 36) and ( 37) are rearranged according to unknowns nodal parameters.The prescribed quantities are multiplied by their influence vectors and added to the right-hand side, and the unknowns and their influence vectors are placed on the left-hand side, thus leading to two standard sets of linear algebraic equations in terms of unknown nodal values as: which can be solved for unknown boundary values at each boundary point.Once these unknown boundary values have been obtained, it is possible to calculate any of the required internal values, using Eqs.( 27) and (32).Moreover, the partial derivatives of w and ϕ can be evaluated at points within domain by direct differentiation of Eqs. ( 27) and (32) for ε P = 1 .

NUMERICAL RESULTS AND DISCUSSION
In order to present the results of developed boundary element method, the three geometrical shapes for a functionally graded plate including, rectangular, circular and elliptic shapes with/without hole are considered.The geometry and dimensions of mentioned plates are given in Fig. 3.

Comparative and convergence studies
To validate the present formulation and ensuring the accuracy and convergence of the proposed boundary element method, several bending analyses of thin plates are first solved and results are compared with those that are available in the literature for isotropic homogenous annular and rectangular plates [25] and functionally graded circular plates [19] subjected to transverse uniformly distributed load.Also, the of convergence is implemented by increasing the boundary elements for the plate.
For the transverse uniformly distributed load, the auxiliary functions v and u presented in Eqs. ( 26) and ( 34) are given as follows: In Table 1, the non-dimensional maximum deflection and resultant moments obtained from the present BEM are compared with the exact closed form solution obtained by Ventsel and Krauthammer [25] for bending analysis of isotropic simply supported rectangular plates with different aspect ratios (L / W ) .The plate has been modeled in two cases with 20 and 40 boundary elements.From this table, it is observed that the present BEM converges exactly when the number of boundary elements increases.Furthermore, in Table 2, the non-dimensional maximum deflection of an isotropic annular plate with different boundary conditions obtained from the present BEM are compared with those obtained from the exact solution given by Ventsel and Krauthammer [25].The annular plate with a clamped (or simply) support at the inner edge (r = R i ) and free at the outer edge (r = R o ) is named as clamped (or simply supported)-free an- nular plate and vise versa.It is seen from Table 2 that the results of proposed BEM in good correlation with the exact results as maximum discrepancy is about 0.59 % when 300 boundary elements have been used.Another comparative study for evaluation of maximum deflection between the present BEM and exact solution developed by Reddy et al. [19] is carried out in Table 3 for FG circular plates with different power of FGMs and boundary conditions.The material properties of titanium-Zirconica FG circular plate (E m / E c = 0.396, ν = 0.288) are taken from Ref. [19].Also, the non-dimensional deflection has been considered as w * = 64wD c FR o 4 which was defined by Reddy et al. [19].It is seen from this table that both results for simply supported FG circular plate are the same for all power of FGMs when 100 boundary elements are used, although errors below 1.4% are observed be-tween results of clamped FG circular plates in this case.This table also shows that a better convergence for clamped FG circular plate is achieved when 200 boundary elements is employed as maximum discrepancy reduces to 0.34 % in this case.This comparison study verifies the validity of functionally graded classical plate theory based on the physical neutral surface concept.Moreover, the stress distribution (σ rr * = σ rr D c FD ) related to the present problem for clamped FG plate has been modeled by the present boundary element method for checking convergency that is presented in Fig. 4.This figure shows the good accuracy for stress distribution when a few number of boundary elements have been employed.Finally, the three comparative studies show that the results obtained from the proposed BEM agree well with exact results of plates with straight and curved simply support edges.The numerical results presented in the next section have been calculated with maximum number of boundary elements after performing the convergence study.

Case studies
In order to obtain new results, it is assumed that the functionally graded plate is made of a mixture of Aluminum (E m = 70GPa) and Alumina (E c = 380GPa) with constant Poisson ratio (0.3) [11].The mechanical load is assumed to be uniform in spatial domain.
) are plotted along the radial direction of FG annular plate, respectively.It is assumed that the annular FG plate with different power of FGMs can have two mixtures of clamped and free boundary conditions at inner and outer edges.Fig. 5 reveals that fully metallic and fully ceramic annular plates have the maximum and minimum value of transverse deflection, respectively.Also, the deflection of functionally graded annular plates occurs between those of full-metal and full-ceramic plates.This event is due to the fact that the stiffness of the metal is lower than that of ceramic so the stiffness of FG plate increases by increasing the power of FGM from zero (full-metal) up to the maximum value for n → ∞ (fullceramic).Related with the variation of non-dimensional radial and tangential resultant moments along the radial direction of FG annular plate, Fig. 6 and 7 show a similar trend as in Fig. 5, namely that for the considered FG plate deflection, the resultant moments of FG annular plates (in absolute sense) decrease with the increase of the power of FGMs.It is worth noting that the tangential resultant moment of free-clamped FG annular plates with any power of FGMs are zero around r / R o = 0.85 .Finally, Figs.5-7 show that free-clamped FG annular plate has smaller deflection and absolute resultant moments in comparison with clamped-free FG annular plate.In fact, the more length clamped constraint at the edges increases the more stiffness of the annular plate, which results in a smaller deflection and resultant moments.
Through the thickness distribution of non-dimensional normal stress (σ rr * = σ rr D c FD ) of free- clamped FG annular plate related to the present example at r / R o = 0.75 is shown in Fig. 8.As it can be seen from this figure, unlike the homogeneous plates (n = 0) , the normal stress of FG plates has not vanished at the middle In fact, the mentioned stress is equal to zero at neutral surface which is located at distance of Δ from the middle surface.The variation of non-dimensional transverse deflection, radial and tangential resultant moments along semi-minor and semi-major axes of FG elliptic plates with free-clamped boundary conditions are illustrated in Fig. 9.The primary inference drawn from Fig. 9 is that the deflection and absolute resultant moments of FG elliptic plates decrease with the increase of the power of FGM.Moreover, it is found that the behavior of resultant moments along semi-minor and semi-major axes of FG elliptic plates is not similar.By focusing on values of the resultant moments along semi-major axis of free-clamped FG elliptic plates presented in Figs 9c and 9e, it can be found that regardless of the power of FGM, the resultant moments are zero near clamped edge.

CONCLUDING REMARKS
In the present article, a pure boundary element method has been developed for bending analysis of functionally graded plates subjected to the transverse load.The material properties have been assumed to vary through the thickness as a power law distribution.Based on the classical plate theory and physical neutral surface concept, the equilibrium equations have been derived.A direct approach based on Green's identity has been used to formulate boundary element method.
By introducing a novel approach and using auxiliary functions domain integrals which arose from distributed loads are transformed into boundary integrals.The bending results have been obtained for a functionally graded plate with various geometrical shapes including, straight and curved boundary domain.Several comparison studies are investigated to reveal accuracy of the present formulation and procedure.Furthermore, the effects of power of FGM and geometric parameters together with various combinations of boundary conditions on the deformation and stress of functionally graded plates have been discussed in detail.Finally the following main results can be concluded: 1) The obtained results from the present BEM agree well with exact results of plate with straight and curved simply support edges.
2) The functionally graded classical plate theory based on the physical neutral surface concept leads to accurate results.3) Fully metallic and fully ceramic annular plates, respectively, have the maximum and minimum value of non-dimensional deflection and absolute resultants moments.Furthermore, these quantities decrease with the increase of power of FGM. 4) The free-clamped FG annular plate has smaller deflection and absolute resultant moments in comparison with clamped-free FG annular plate.5) Unlike the homogeneous plates, for which their normal stress is vanished at the middle surface, normal and stress of functionally graded plate is equal to zero at neutral surface.6) The resultant moments along semi-major axis of free-clamped FG elliptic plates are zero near clamped edge.

Figure 1
Figure 1 The position of middle surface and neutral surface for a functionally graded plate.

δψ 1 :
N xx,x + N xy,y = 0 δψ 2 : N xy,x + N yy,y = 0 δw : M xx,xx + 2M xy,xy + M yy,yy + F = 0 (10) where d represents the variational symbol; F = F(x,y) is the mechanical load per unit area, and N and M are resultant forces and resultant moments, respectively, which are defined by the following expressions

Figure 2
Figure 2 Geometry and modeling of boundary domain: (a) Coordinate system, (b) Modeling of the boundary with superparametric elements, (c) Definition of nodal-point location, relative distance and relative angle.

Figure 4
Figure 4 Convergence test of the non-dimensional maximum stress ) ( * FD D rr c rr σ σ =

Figure 5 = 2 (
Figure 5 Non-dimensional deflection * w along the radial direction of FG annular plate with different boundary conditions for various powers of FGM R 0 R i = 2 ( ) : (a) clamped-free, (b) free-clamped.

Figure 6 2 (Figure 7 2 (
Figure 6 Non-dimensional radial resultant moment * rr M along the radial direction of FG annular plate with different boundary conditions for various powers of FGM R 0 R i = 2 ( ) : (a) clamped-free, (b) free-clamped.

Figure 8 Figure 9
Figure 8 Non-dimensional normal stress σ rr * through the thickness of free-clamped FG annular plate for various powers of FGM R 0 R i = 2, r R 0 = 3 4 () .

Table 1
Convergence test of the non-dimensional maximum deflection and resultant moments for isotropic simply supported rectangular plates.
* The number of boundary elements used in the present BEM † The discrepancy Table 2 Convergence test of the non-dimensional maximum deflection (w * = 64 w D