Arash Shahbaztabar·Koosha Arteshyar

Abstract We extend the differential quadrature element method(DQEM)to the buckling analysis of uniformly in-plane loaded functionally graded (FG) plates fully or partially resting on the Pasternak model of elastic support. Material properties of the FG plate are graded in the thickness direction and assumed to obey a power law distribution of the volume fraction of the constituents.To set up the global eigenvalue equation,the plate is divided into sub-domains or elements and the generalized differential quadrature procedure is applied to discretize the governing,boundary and compatibility equations.By assembling discrete equations at all nodal points, the weighting coefficient and force matrices are derived. To validate the accuracy of this method,the results are compared with those of the exact solution and the finite element method.At the end,the effects of different variables and local elastic support arrangements on the buckling load factor are investigated.

Keywords Differential quadrature element method·Pasternak elastic support·Critical buckling load·Functionally graded plates

1 Introduction

Functionally graded materials (FGMs), introduced by a group of scientists working on a space plane project [1] in the mid-1980s, are a class of composites usually produced by mixing ceramic and metal materials. The smooth stress distribution property of FGMs,a product of the continuous change in the material properties from one surface to another,gives them precedence over laminated composites. A wide range of research projects related to the buckling analysis of FG structural members has been carried out. Based on the classical plate theory(CPT)and applying the Fourier series method and Stokes' transformation, Latifi et al. [2] studied buckling of thin FG plates.In the framework of the first and higher order shear deformation plate theory(FSDT/HSDT),mechanical and thermal buckling of FG plates have also been investigated by applying various numerical methods,local Kriging meshless method [3], isogeometric analysis[4], mesh-free method [5,6], finite element method [7,8],extended Kantorovich method[9],Chebyshev-Ritz method[10]and Chebyshev spectral collocation method[11].Based on the third shear deformation plate theory(TSDT),an exact solution was presented by Thai and Kim [12] for buckling analysis of FG plates on Pasternak elastic support.An analytical solution was carried out by Bodaghi and Saidi [13]for buckling behavior of FG rectangular plates based on the HSDT. By applying an inverse trigonometric shear deformation theory (ITSDT), an analytical solution for bending and buckling of functionally graded plates was developed by Kulkarni et al.[14].

It is well-known that difficulties in applying the analytical methods arise when complex geometries,patch loading,local elastic supports,mixed boundary conditions and other discontinuities are taken into account.Therefore,the use of numerical approaches is unavoidable. Differential quadrature method(DQM)is a powerful numerical method which was developed by Bellman and Casti [15] and can take precedenceovertheconventionalfiniteelementandfinitedifference, two well-known numerical methods that are often computationally expensive. The pioneer scholars adopting the DQM in structural problems were Bert et al. [16]. The early applications of DQM included many restrictions；however,investigators have made valuable efforts to rectify these limitations so that differential quadrature method has become a highly efficient and accepted numerical approach for the solution of various engineering problems. For example,restrictions such as calculating the weighting coefficients,were improved by Shu and Richards[17].Implementation of multiple boundary conditions and discontinuities in geometry and loading led to the introduction of various methods such as δ-technique [18], the modified weighting coefficient approach[19],the equation replaced approach[20],the quadrature element method(QEM)[21]and the differential quadrature element method(DQEM)[22-26].

The study of buckling behavior of plates made of functionally graded materials on elastic supports is confined to problems in which the plate is completely on the elastic foundation.Many engineering systems,however,can be considered as a plate that is partially resting on an elastic support.By using the differential quadrature method, Nobakhti and Aghdam [27] and Jahromi et al. [28] studied bending and free vibration of isotropic plates on local elastic supports.Applying this method depends strongly upon a correct distribution of the grid points in such a way that the interface of the elastic support would be covered.This drawback can be removed by applying DQEM. In the present work, the application of DQEM is developed to investigate the buckling characteristics of uniformly in-plane loaded functionally graded plates fully or locally resting on an elastic foundation. To set up the global eigenvalue equation, the plate is divided into sub-domains or elements, and the generalized differential quadrature procedure is adopted to the governing, boundary and compatibility equations. By assembling discrete equations at all nodal points, the weighting coefficient and force matrices are derived. The validity of the present formulation is confirmed via a comparison study with the analytical solution reported in the literature, and the impact of various parameters and local elastic support arrangements on the critical buckling load factor is illustrated.

2 Governing equations

2.1 Kinematics

Consider a rectangular plate made of functionally graded materials resting fully or partially on a Pasternak elastic support and subjected to uniformly distributed in-plane loads,as shown in Fig.1.According to the FSDT,the components of displacements may be defined as

Fig.1 Geometry and coordinates of an FG plate subjected to uniform compression loads on the Pasternak foundation

where u0, v0, and w0are displacement components of an arbitrary point in the middle surface；ψxand ψyare rotation angles about the y and x axes,respectively. The strains can be expressed in terms of displacements as

with

2.2 Stress-strain relations

The Young's modulus(E)and the mass density(ρ)of the FG plate are considered to be graded in the thickness direction and to obey a power law distribution as

where subscripts“m”and“c”denote the metal and ceramic constituents,respectively,and α denotes the power law index.

The stresses can be related to the strains by the generalized Hooke's law as

The quantities Cijare given by

where υ stands for Poisson's ratio.In the present formulation,this quantity is assumed to have a constant value of υ =0.3.

2.3 Equations of motion

By applying Hamilton's principle, the following governing differential equations can be derived for vibration and buckling of FG plates on the elastic foundation

Nx0and Ny0are the normal in-plane forces per unit length and Nxy0denotes the shear in-plane force per unit length.k1represents the Winkler coefficient and k2is the shear stiffness parameter of the foundation. The force and moment resultants and the mass moments of inertia(I0,I1,I2)can be obtained as follows:

The shear correction factor (KS) is taken as 5/6 in this study；this value can present stable and accepted results for FG plates.The governing Eq.(7)can be written in terms of the displacement components by applying Eqs.(2),(5)and(8)as

where

The displacement functions of the FG plate are assumed to be:

where ω denotes the circular (natural) frequency. U(x,y),V(x,y),W(x,y),Ψx(x,y)and Ψy(x,y)are unknown functions. The uniform normal and shear forces per unit length are assumed to be Nx0= -γ1Ncr, Ny0= -γ2Ncrand Nxy0= 0.The intensity of normal forces is assumed to be negative in compression.Substituting Eq.(11)into Eq.(9),one can obtain the following dimensionless governing equations:

The non-dimensional terms in Eq.(12)are defined as follows:

In the present study, the boundary conditions along the plate edges at x=0,a and y=0,b are defined as follows:

· Simply supported(S):

· Clamped(C):

· Free(F):

3 Application of the DQEM

To employ the differential quadrature element method,the plate should be divided into elements according to the arrangements of the local elastic support. Each of these elements is further discretized into equally or non-equally distributed grid points by using the generalized differential quadrature method. To do so, the discrete form of the governing differential Eq.(12)for each element can be written as:

where nx and ny are the total number of nodal points in an element with respect to the ξ and η directions,respectively,and ne represents the number of elements.are weighting coefficients of the first and second-order partial derivatives of U,V,W,Ψxand Ψywith respect to ξ and η,respectively.The weighting coefficients can be calculated by the following equations[20,22]:

where

The following non-uniform distributed mesh points can be used in the above formulas:

It should be noted, the weighting coefficients of each element can be simply obtained by modifying Eq. (16) as follows:

where ly and lx are,respectively,the width and length of the relevant elements in the computational domain.For example,thelengthandwidthof eachelement for aplatewiththreeelements shown in Fig.2b,can be taken as lx1=lx2=lx3=1 and ly1=0.3, ly2=0.4, ly3=0.3.To complete the DQE solution of an arbitrary element,two different sets of equations,either boundaryconditions or compatibilityconditions,should be added to Eq.(15).More details can be found in the work presented by Liu and Liew [25]. For a FG plate, the compatibility conditions can be expressed:

· At common nodes of two elements(Fig.3a,b):for 0＜x ＜1 and j=1,2,…,ny,

for 0＜y ＜1 and i=1,2,…,nx,

· At common node of four elements(Fig.3c):

or

Fig.3 Connected grid points at interface boundaries of adjoining elements

The discrete forms of the stress resultants in the above equations are as follows:

To obtain the critical values of the load factor(Δ)or natural frequency parameters (β) of the plate, the weightingcoefficients,mass and force matrices of all elements should be assembled into a global system as follows:

Table 1 Material properties of functionally graded plate

where KTdenotes the total weighting coefficients matrix,BTis the total force matrix, MTis the global mass matrix,and Λ represents the global displacements vector.

4 Numerical results

Computer code in MATLAB was developed to solve the eigenvalue problem and to obtain the critical buckling load factor(Δ),given in Eq.(25).First,comparison studies were carried out to show the accuracy of the present method compared to the analytical results and the finite element solution.Then, parametric studies were performed to examine the effects of various parameters on the critical buckling load.For convenience,a 4-capital-letter combination of S,C and F is used to define the boundary conditions.S,C and F denote the simply supported, clamped and free boundary conditions,respectively.For instance,SCSF defines a plate which has simply supported boundary condition at edges ξ =0,1,clamped at η=0 and free at η=1.Two materials,aluminum(Al) and stainless steel (SUS304), were considered as the metal constituents and alumina(Al2O3)as the ceramic.The properties of these materials are shown in Table 1.

Table 2 Comparison of buckling load parameter Δ= Ncra2/(Emh3)for Al/Al2O3 square plates(λ=1)

Table 3 Comparison of buckling load parameter Δ = Ncra2/(Emh3) for SSSS Al/Al2O3 square plates on Pasternak elastic foundation under uniaxial compression(λ=1, γ1 =-1, γ2 =0)

4.1 Comparison studies

In the first example,non-dimensional critical buckling load factors Δ = Ncra2/(Emh3)of Al/Al2O3square plates(λ=1) subjected to uniaxial (γ1= -1, γ2= 0) and biaxial(γ1=-1, γ2=-1)compressive loads were obtained for three different types of boundary conditions.Different values of δ and power law index(α)were selected.The numerical data were achieved by only eleven grid points in each direction (nx =ny =11) for a plate element and are reported in Table 2.This table also encompasses the results of the exact solution based on the third shear deformation plate theory,which does not need shear correction factor,obtained by Thai and Kim[12].As can be observed from Table 2,the results are in excellent agreement, even for thicker FG plates with large values of thickness to length ratio.

In the next case, the critical values of the buckling load factor Δ = Ncra2/(Emh3) for simply supported Al/Al2O3square plates(λ=1)with full elastic support and subjected to uniaxial (γ1= -1, γ2= 0) compressive load are compared in Table 3 with those of the exact solution presented by Thai and Kim[12]based on the TSDT.As can be inferred, the present results are very close to the analytical solutions.

Assessment of the critical values of buckling load (Ncr)of a simply supported isotropic plate partially on the Winkler elastic foundation (= 100,= 0) and subjected to uniaxial compressive load (γ1= -1, γ2= 0) is the last example to check the accuracy of the present formulation. Because the solution related to this problem is not reported in the open literature,the present numerical results were compared with those obtained by finite element software ABAQUS(Release 6.13)and summarized in Table 4.The dimensions and material properties of the plate in this analysis were taken as a = 1 m, b = 1 m, h = 0.01 m,E =200 GPa,υ=0.3.For finite element analysis,a quadrilateral element (S4R) was used. It is evident from Table 4 that good agreement can be observed for all arrangements of the local elastic supports.

Table 4 Comparison of critical buckling load Ncr for isotropic plates partially resting on Winkler elastic foundation and subjected to uniaxial compression (a = b = 1, h = 0.01, tx1 = tx2 = ty1 = ty2 =0.2, γ1 =-1, γ2 =0, ~k1 =100, =0)

Table 4 Comparison of critical buckling load Ncr for isotropic plates partially resting on Winkler elastic foundation and subjected to uniaxial compression (a = b = 1, h = 0.01, tx1 = tx2 = ty1 = ty2 =0.2, γ1 =-1, γ2 =0, ~k1 =100, =0)

4.2 Parametric survey

In the previous section,it was confirmed that the differential quadrature element method,in spite of applying a small number of nodal points,can provide accurate solutions to the structural problems.In the following section,we investigate the effects of various quantities including aspect ratio (λ),power law index(α),foundation stiffness parameters,thickness to length ratio(δ)and different types of locally elastic support arrangements (Fig. 2), on the critical load factor.Figure 2 depicts three different configurations of the locally supported foundations including boundary strip elastic supports along the edges of the FG plate(Fig.2a,Case A),local elastic support along two parallel edges(Fig.2b,Case B)and local elastic supports at four corners of the FG plate(Fig.2c,Case C). Manhole covers, ship deck hatches and offshore platforms are some typical examples of these cases；there are many engineering systems that could be considered as a plate partially resting on elastic supports.These particular arrangements of partially supported foundations can illustrate how the differential quadrature element method can be used for different cases. For example, selecting the number of elements depends on the local elastic foundation arrangements.For a plate with partially elastic supports along two parallel edges,the number of elements that can be used is three(ne=3),while,for a plate with elastic supports at four corners,nine elements(ne=9)should be employed.

Fig.4 Variation of critical buckling load factor for FG plates fully and partially supported on Pasternak foundation versus the gradient index under uniaxial compression (SSSS, δ = 0.1, λ = 1, ~k1 =1000, 100, Al/Al2O3, γ1 =-1, γ2 =0)

Fig.5 Variation of critical buckling load factor for FG plates fully and partially supported on Pasternak foundation and subjected to biaxial compression versus δ ratio(CCCC, λ=1, α =1, ~k1 =1000, ~k2 =100, SUS304/Al2O3)

Fig.6 Variation of critical buckling load factor for FG plates partially on Pasternak foundation under biaxial compression versus the aspect ratio(δ = 0.1, α = 1,= 1000, = 100, tx1 = tx2 = ty1 =ty2 =0.3, SUS304/Al2O3)

Figure 4 depicts the variation of the critical buckling load factor Δ = Ncra2/(Emh3) versus the gradient index for a simply supported Al/Al2O3square plate (λ =1) subjected to uniaxial compressive (γ1= -1, γ2= 0) load and fully or partially resting on Pasternak elastic strips along the edges. The values of the stiffness parameters of the elastic foundation are taken as ~k1= 1000, ~k2= 100 and tx1=tx2= ty1= ty2= tp.From Fig.4,it can be inferred that by increasing the power law index,the critical buckling load factor of FG plates supported fully(tp=0.5)or locally on an elastic foundation is decreased.

Variation of the critical buckling load factor Δ =Ncra2/(Emh3)as a function of thickness to length ratio(δ)is plotted in Fig.5 for a CCCC SUS304/Al2O3square plate(λ = 1) fully (ty = 0.5) and partially resting on Pasternak elastic supports along two parallel edges (ty1= ty2=tp, tx1= tx2= 0) and subjected to biaxial compressive load(γ1=-1, γ2=-1).The foundation parameters are=1000,=100 and α=1.It is concluded that by taking large values of δ,critical buckling load factors decrease.Furthermore, it can be found that by increasing the contact area between the plate and foundation from tp=0.1 to tp=0.5,larger critical buckling load factors are obtained.

Table 5 Buckling load factor Δ= Ncra2/(Emh3)for SSSS Al/Al2O3 plates fully on elastic foundation

Table 5 continued

Figure 6 illustrates the variation of the critical buckling load factor Δ = Ncra2/(Emh3) with respect to the aspect ratio for a SUS304/Al2O3square plate (λ = 1) fully (ty= 0.5) and partially resting on Pasternak elastic supports at four corners and subjected to biaxial compressive load(γ1=-1, γ2=-1),when α=1 and tx1=tx2=ty1=ty2=0.3.Six combinations of boundary conditions are considered.It is concluded that the critical values of load factor Δ=Ncra2/(Emh3)for boundary conditions SSSS,SCSC and CCCC can be increased by enhancing the aspect ratio.However,an unstable behavior(increasing-decreasing trend)for boundary conditions having a free edge can be observed.It is also found that the buckling load factor(Δ)for a FG plate with one or more free edges is less sensitive to the values of aspect ratio, and by imposing high restraining boundary conditions a dramatic change can be observed.

The critical load factors Δ=Ncra2/(Emh3)for FG plates completely laid on the elastic supports under uniaxial and biaxial compressive loads are reported in Tables 5 and 6.The values of critical load parameter for two configurations of local elastic supports are summarized in Tables 7 and 8.A wide range of values for the power law index,aspect ratios,foundation parameters and thickness to length ratios have been considered.It can be seen from these tables that extending the foundation stiffness coefficients and the contact area between the plate and foundation, regardless of boundary conditions,strongly influence on the buckling load parameter.

5 Conclusions

In this paper, we investigated the buckling behavior of uniaxial and biaxial in-plane loaded functionally graded rectangular plates fully or partly on the Pasternak elastic support by applying the DQ element method. Comparison studies proved that this method can provide accurate solution of the structural problems in spite of applying relatively few nodal points.After demonstrating the precision of the DQ element formulation,the effects of different quantities on the buckling load parameter were studied in detail.The findings reveal the following conclusions.

(1) With increasing values of power law index(α)and thickness to length ratios,the critical buckling load parameter decreases.

(2) The variation of the critical load parameter with respect to the aspect ratios for a plate having one or more free edges is small, and dramatic changes can be observed by imposing high restraining boundary conditions.

(3) In addition to enhanced values of stiffness parameters of elastic support,as well as increasing the aspect ratios and interaction area between the plate and foundation,the critical buckling load parameters increase.

Table 6 Buckling load factor Δ= Ncra2/(Emh3)for CCCC Al/Al2O3 plates fully on elastic foundation

Table 6 continued

Table 7 Buckling load factor Δ= Ncra2/(Emh3)for SSSS Al/Al2O3 square plates partially on elastic supports(Case A,λ=1,δ=0.1)

Table 8 Buckling load factor Δ= Ncra2/(Emh3)for SFSF SUS304/Al2O3 square plates partially on elastic supports(Case C,λ=1,δ=0.1)