# Analysis of Sandwich Plates by Generalized Differential Quadrature Method

Authors: A. J. M. Ferreira; E. Viola; F. Tornabene; N. Fantuzzi; A. M. Zenkour
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We combine a layer-wise formulation and a generalized differential quadrature technique for predicting the static deformations and free vibration behaviour of sandwich plates. Through numerical experiments, the capability and efficiency of this strong-form technique for static and vibration problems are demonstrated, and the numerical accuracy and convergence are thoughtfully examined.


---

## Body

## 1. Introduction

The two big groups of theories for the analysis of laminated plates are based on displacements only or on displacements and stresses at interfaces of the laminate. This paper deals with the first group, where again two approaches can be considered: the equivalent-single layer (ESL) theories (where all the layers are referred to the same degrees of  freedom) or layer-wise (LW) theories (where specific degrees of freedom are linked to specific layers).The classical laminated plate theory and the first-order shear deformation theory [1–5] are the simplest ESL theories, at the cost of introducing artificial shear correction factors. ESL higher-order theories usually consider higher-order expansion of displacements with respect to the thickness coordinate of the plate, with the further advantage of disregarding shear correction factors [6–12]. Most of the previously mentioned ESL theories may ill-represent the transverse shear stresses for sandwich laminates, where the material properties of skins are typically much larger than those of the core. In such problems, LW theories that consider independent degrees of freedom for each layer should be considered instead. Several LW theories have been proposed in literature [13–21] usually considering translational degrees of freedom at the layer interfaces. Over the last decades, an automatic approach for generation of any kind of ESL or LW theories was developed by Carrera [18–20], which allows the user to freely develop new theories without much effort. Various examples of the use of Carrera’s unified formulation with meshless methods by the first author can be found in [22–27].Unlike LW theories that consider translations at each laminate interface, in this work we adopt a layer-wise theory that considers both translations and rotations, based on an expansion of Mindlin’s first-order shear deformation theory in each layer. The displacement continuity at layer’s interface is enforced. In each layer, the theory produces constant, very accurate transverse shear stress.As is well known, the analysis of composite plates is typically performed by finite element method. On the contrary, in the present study, a layer-wise shear deformation theory is implemented within generalized differential quadrature technique. The latter meshfree method in combination with this layer-wise theory carries out an accurate calculation of the natural frequencies, displacements, and stresses of sandwich plates. This study appears for the first time and serves to fill the gap of knowledge in this research topic.


## 2. A Layer-Wise Theory

The layer-wise theory proposed in this paper is based on the assumption of a first-order shear deformation theory in each layer and the imposition of displacement continuity at layer’s interfaces. Given that, a layer-wise formulation links the degrees of freedom at each layer. Nevertheless, the problem may become large if a large number of layers are considered, although the present formulation is adequate for a general laminate. In this paper, we concentrate on sandwich plates built from two external, stiff layers and one central, soft core. The three-layer laminate is shown schematically in Figure1.Figure 1
1D representation of the layer-wise kinematics.In order to simplify the imposition of displacement continuity at the core-skin interfaces, we first define the displacement field for the middle layer (the core in the sandwich plate) as(1)u(2)(x,y,z)=u0(x,y)+z(2)θx(2)(x,y),(2)v(2)(x,y,z)=v0(x,y)+z(2)θy(2)(x,y),(3)w(2)(x,y,z)=w0(x,y),
where u and v are the in-plane displacements at any point (x,y,z), u0 and v0 are the in-plane displacements of the points (x,y,0) on the midplane, w is the deflection, and θx(2) and θy(2) are the rotations of the normals to the midplane about the y- and x-axes, respectively.The corresponding displacement fields for the (skins) upper layer (3) and lower layer (1) are given, respectively, as
(4)u(3)(x,y,z)=u0(x,y)+h22θx(2)(x,y)+h32θx(3)(x,y)+z(3)θx(3)(x,y),v(3)(x,y,z)=v0(x,y)+h22θy(2)(x,y)+h32θy(3)(x,y)+z(3)θy(3)(x,y),w(3)(x,y,z)=w0(x,y),u(1)(x,y,z)=u0(x,y)-h22θx(2)(x,y)-h12θx(1)(x,y)+z(1)θx(1)(x,y),v(1)(x,y,z)=v0(x,y)-h22θy(2)(x,y)-h12θy(1)(x,y)+z(1)θy(1)(x,y),w(1)(x,y,z)=w0(x,y),
where hk are the kth layer thickness and z(k)∈[-hk/2,hk/2] are the kth layer z coordinates.Deformations for a generic layerk are given by
(5){ϵxx(k)ϵyy(k)γxy(k)γxz(k)γyz(k)}={∂u(k)∂x∂v(k)∂y∂u(k)∂y+∂v(k)∂x∂u(k)∂z+∂w(k)∂x∂v(k)∂z+∂w(k)∂y}.
It is common that the ESL theories decouple the in-plane deformations with the out-of-plane deformations. The in-plane deformations can be expressed as
(6){ϵxx(k)ϵyy(k)γxy(k)}={ϵxxm(k)ϵyym(k)γxym(k)}+z(k){ϵxxf(k)ϵyyf(k)γxyf(k)}+{ϵxxmf(k)ϵyymf(k)γxymf(k)}
and shear deformations as
(7){γxz(k)γyz(k)}={∂w0∂x+θx(k)∂w0∂y+θy(k)}.
The membrane components are given by
(8){ϵxxm(k)ϵyym(k)γxym(k)}={∂u0∂x∂v0∂y∂u0∂y+∂v0∂x},
where superscript m denotes membrane components.The bending components can be expressed as(9){ϵxxf(k)ϵyyf(k)γxyf(k)}={∂θx(k)∂x∂θy(k)∂y∂θx(k)∂y+∂θy(k)∂x},
where superscript f denotes flexural or bending components and the membrane-bending coupling components for layers 2, 3, and 1 are, respectively, given as
(10){ϵxxmf(2)ϵyymf(2)γxymf(2)}={000},{ϵxxmf(3)ϵyymf(3)γxymf(3)}={h22∂θx(2)∂x+h32∂θx(3)∂xh22∂θy(2)∂y+h32∂θy(3)∂yh22(∂θx(2)∂y+∂θy(2)∂x)+h32(∂θx(3)∂y+∂θy(3)∂x)},{ϵxxmf(1)ϵyymf(1)γxymf(1)}={-h22∂θx(2)∂x-h12∂θx(1)∂x-h22∂θy(2)∂y-h12∂θy(1)∂y-h22(∂θx(2)∂y+∂θy(2)∂x)-h12(∂θx(1)∂y+∂θy(1)∂x)},
where superscript mf denotes membrane-flexural components.The present layer-wise formulation neglects the transverse normal deformations, by assuming a constant transverse displacement in every layer. Neglectingσz(k) for each orthotropic layer, the stress-strain relations in the fiber local coordinate system can be expressed as
(11){σ1(k)σ2(k)τ12(k)τ23(k)τ31(k)}=[Q11(k)Q12(k)000Q12(k)Q22(k)00000Q33(k)00000Q44(k)00000Q55(k)]{ε1(k)ε2(k)γ12(k)γ23(k)γ31(k)},
where subscripts 1 and 2 are, respectively, the principal middle surface directions, 3 is the direction normal to the middle plane of the plate, and the reduced stiffness components, Qij(k), are given by
(12)Q11(k)=E1(k)1-ν12(k)ν21(k),Q22(k)=E2(k)1-ν12(k)ν21(k),Q12(k)=ν12(k)Q22(k),Q33(k)=G12(k),Q44(k)=G23(k),Q55(k)=G31(k),ν21(k)=ν12(k)E2(k)E1(k),
in which E1(k), E2(k), ν12(k), G12(k), G23(k), and G31(k) are material properties of the lamina k. Note that shear correction factors are not considered in the stress-strain laws, as they must be done in the first-order shear deformation theories [21].The stress-strain relations in the globalx-y-z coordinate system can be obtained by coordinate transformation as
(13){σxx(k)σyy(k)τxy(k)τyz(k)τzx(k)}=[Q¯11(k)Q¯12(k)Q¯16(k)00Q¯12(k)Q¯22(k)Q¯26(k)00Q¯16(k)Q¯26(k)Q¯66(k)00000Q¯44(k)Q¯45(k)000Q¯45(k)Q¯55(k)]{εxx(k)εyy(k)γxy(k)γyz(k)γzx(k)}.
By considering θ as the angle between x-axis and 1-axis, with 1-axis being the first principal material axis, connected usually with fiber direction, the components Q¯ij(k) can be computed by coordinate transformation (see [13] for details).The equations of motion of this layer-wise theory are derived from the dynamic version of the principle of virtual displacements. In the present work, only symmetric laminates are considered; therefore,u0, v0, and related stress resultants can be discarded. The virtual strain energy (δU), the virtual kinetic energy (δK), and the virtual work done by applied forces (δV), assuming a three-layer laminate, are given by
(14)δU=∫Ω0‍∑k=13‍{∫-hk/2hk/2‍[σxx(zδϵxxf(k)+δϵxxmf(k))+σyy(zδϵyyf(k)+δϵyymf(k))+τxy(zδγxyf(k)+δγxymf(k))+τxzδγxz(k)+τyzδγyz(k)(zδεxxf(k)+δεxxmf(k))]∫-hk/2hk/2‍[σxx(zδεxxf(k)+δεxxmf(k))dz}dxdy=∫Ω0‍∑k=13‍(Nxx(k)δϵxxmf(k)+Mxx(k)δϵxxf(k)+Nyy(k)δϵyymf(k)+Myy(k)δϵyyf(k)+Nxy(k)δγxymf(k)+Mxy(k)δγxyf(k)+Qx(k)δγxz(k)+Qy(k)δγyz(k))dxdy.
It should be noted that the integrals in the thickness direction are now in the layer k domain, from -hk/2 to hk/2. Consider
(15)δK=∫Ω0‍∑k=13‍∫-hk/2hk/2‍ρ(k)×(u˙(k)δu˙(k)+v˙(k)δv˙(k)+w˙(k)δw˙(k))dzdxdy,
where boldface quantities are vectors and (·) denotes time differentiation. Finally, the virtual work is given as
(16)δV=-∫Ω0‍qδw0dxdy,
where δ quantities are virtual, Ω0 denotes the midplane of the laminate, q is the external distributed load, and
(17)(Nαβ(k),Mαβ(k))=∫-hk/2hk/2‍σαβ(k)(1,z)dzk,Qα(k)=∫-hk/2hk/2‍ταz(k)dzk,
where α, β take the symbols x,y.SubstitutingδU, δK, and δV into the virtual work statement, using the fundamental lemma of the calculus of variations, and integrating by parts to relieve from any derivatives of the generalized displacements (see [13] for details), the equations of motion with respect to 7 degrees of freedom (w0, θx(1), θy(1), θx(2), θy(2), θx(3), and θy(3)) are obtained. Consider
(18)δw0:∑k=13‍(∂Qx(k)∂x+∂Qy(k)∂y)-q=∑k=13‍I0(k)w¨0,δθx(1):h12∂Nxx(1)∂x-∂Mxx(1)∂x+h12∂Nxy(1)∂y-∂Mxy(1)∂y+Qx(1)=I0(1)(h1h24θ¨x2+h124θ¨x1)+I2(1)θ¨x1,δθy(1):h12∂Nyy(1)∂y-∂Myy(1)∂y+h12∂Nxy(1)∂x-∂Mxy(1)∂x+Qy(1)=I0(1)(h1h24θ¨y2+h124θ¨y1)+I2(1)θ¨y1,δθx(2):h22∂Nxx(1)∂x-h22∂Nxx(3)∂x-∂Mxx(2)∂x+h22∂Nxy(1)∂y-h22∂Nxy(3)∂y-∂Mxy(2)∂y+Qx(2)=I0(1)(h224θ¨x2+h1h24θ¨x1)+I0(3)(h224θ¨x2+h2h34θ¨x3)+I2(2)θ¨x2,δθy(2):h22∂Nyy(1)∂y-h22∂Nyy(3)∂y-∂Myy(2)∂y+h22∂Nxy(1)∂x-h22∂Nxy(3)∂x-∂Mxy(2)∂x+Qy(2)=I0(1)(h224θ¨y2+h1h24θ¨y1)+I0(3)(h224θ¨y2+h2h34θ¨y3)+I2(2)θ¨y2,δθx(3):-h32∂Nxx(3)∂x-∂Mxx(3)∂x-h32∂Nxy(3)∂y-∂Mxy(3)∂y+Qx(3)=I0(3)(h2h34θ¨x2+h324θ¨x3)+I2(3)θ¨x3,δθy(3):h32∂Nyy(3)∂y-∂Myy(3)∂y-h32∂Nxy(3)∂x-∂Mxy(3)∂x+Qy(3)=I0(3)(h2h34θ¨y2+h324θ¨y3)+I2(3)θ¨y3,
where
(19)(I0(k),I2(k))=∫-hk/2hk/2‍ρ(k)(1,z2)dz,
where ρk, the specific mass of the material of the k is layer and hk is the thickness of the kth layer.


## 3. The Generalized Differential Quadrature Method

The differential quadrature (DQ) method was presented by Bellman and his associates in the early 1970s [28]. It is a numerical discretization technique for the approximation of derivatives. The DQ method was initiated from the idea of conventional integral quadrature. In fact, one problem which arises frequently in structural mechanics and in many other engineering problems is the evaluation of the integral ∫ab‍f(x)dx. If F is a function such that dF/dx=f, then the value of the given integral is F(b)-F(a). Unfortunately, in practical problems, it is extremely difficult, and most of the times impossible, to obtain an explicit expression for F. The values of f, perhaps, can be known at a discrete set of points, and, in this situation, a numerical approach is essential.Following the idea of integral quadrature, Bellman et al. [28] suggested that the first-order derivative of the function f(x) with respect to x at a grid point xi is approximated by a linear sum of all the functional values in the whole domain as
(20)df(x)dx|x=xi=∑j=1N‍aijx,(1)f(xj)fori=1,2,…,N,
where aijx,(1) represent the weighting coefficients, N is the total number of grid points x1,x2,…,xN in the whole domain, and f(xj) is the calculated value of f(x) at the point x=xj. However, the weighting coefficients which were used depend, in the Bellman approach, on the number of grid points and the type of discretization that was used. Some years later the DQ approach has been extended by Shu [29], who introduced the generalized differential quadrature (GDQ) method. In [29], the weighting coefficients calculation is performed through recursive formulae, and they are independent of the number of grid points N and the discretization type. The interested reader can find a brief review on GDQ applications in [30–40].To compute the first-order weighting coefficientsaijx,(1), the following algebraic formulae are derived:
(21)aijx,(1)=L(1)(xi)(xi-xj)L(1)(xj)fori,j=1,2,…,N,i≠jaiix,(1)=-∑k=1,k≠iN‍aik(1)fori=j,
where L(1)(xk)=∏j=1,j≠kN‍(xk-xj), which are the derivatives of the Lagrange polynomials L(x)=(x-x1)(x-x2)⋯(x-xN)=∏j=1N‍(x-xj). As reported above, Shu [29] found out a recurrence formulation for the weighting coefficients calculation of the second- and higher-order derivatives. Consider
(22)aijx,(n)=n(aiix,(n-1)aijx,(1)-aijx,(n-1)xi-xj)fori≠j,n=2,3,…,N-1,aiix,(n)=-∑k=1,k≠iN‍aikx,(n)fori=j.
So any derivative of an order greater than the first can be written as
(23)d(n)f(x)dxn|x=xi=∑j=1N‍aijx,(n)f(xj)fori=1,2,…,N,n=2,3,…,N-1.
The previous derivative discretization was presented for the one-dimension problems only. Nevertheless, the simple one-dimensional case can be directly extended to the multidimensional one for any regular shape, such as a rectangular or circular. Given a two-dimensional physical system, in Cartesian coordinates, described by a function f(x,y), the problem values depend on the nodal coordinates x and y. Since GDQ works with regular grid along the main axes, N and M will indicate the points along x and y, respectively. Using the discretization of the derivatives in the one-dimensional case (23), a generic higher-order derivative along x and y for the two-dimensional case can be reported. Consider
(24)∂(n)f(x,y)∂xn|x=xi,y=yj=∑k=1N‍aikx,(n)f(xk,yj)fori=1,2,…,N,n=1,2,…,N-1,∂(m)f(x,y)∂ym|x=xi,y=yj=∑l=1M‍ajly,(m)f(xi,yl)forj=1,2,…,M,m=1,2,…,M-1,
where aikx,(n) and ajly,(m) are the weighting coefficients of order n and m along x and y, respectively. Moreover, the mixed derivative can be written in the same way as (24). Consider
(25)∂(n+m)f(x,y)∂xn∂ym|x=xi,y=yj=∑k=1N‍aikx,(n)(∑l=1M‍ajly,(m)f(xk,yl))fori=1,2,…,N,j=1,2,…,M,n=1,2,…,N-1,m=1,2,…,M-1,
where aikx,(n) and ajly,(m) have the same meaning in (24).


### 3.1. Grid Distributions

As it has been widely proven [34–37], the GDQ numerical accuracy strongly depends on the grid distribution choice for a certain numerical problem. On the one hand, the natural and simplest choice of uniform grids does not lead to stable and accurate results for any GDQ computation; on the other hand, the Chebyshev-Gauss-Lobatto (C-G-L) grid point distribution gives accurate results in most of cases. In the numerical results proposed in this paper a C-G-L grid distribution is considered for both directions x and y. Consider
(26)xi=12(1-cos(i-1N-1)π),yj=12(1-cos(j-1M-1)π),
where N, M are the total number of points in the given directions x, y, respectively.


### 3.2. The Static Problem

The partial differential system of equations developed in Section2 has been discretised and numerically solved using the GDQ method. Since GDQ allows turning any derivative in an algebraic expression, the fundamental system of equations is solved in its strong form directly, within its boundary conditions. A static boundary problem can be indicated as follows:
(27)ℓ(u(x))+q(x)=0inΩ,ℓB(u(x))=gon∂Ω,
where ℓ is a partial differential operator and contains all the spatial derivatives of the problem with respect to a Cartesian coordinate system x. The model degrees of freedom are u(x), and the external forces acting on the domain Ω are indicated by q(x). It is obvious that the differential problem cannot be solved without its boundary conditions, where ℓB is the differential operator related to the boundary conditions of the given system of equations and g are the boundary loads.Once the boundary-value problem (27) is discretised by GDQ the following set of algebraic equations is obtained:
(28)ℒU+Q=0,ℒBU=G,
where U is a vector representing a set of unknown functional values at all the interior points, ℒ is a matrix carried out from the GDQ method, and Q, G are the known vectors arising from the functions q(x) and g.In order to obtain the global stiffness matrix of the physical problem under study, the discretised boundary conditions must be embedded into the discrete fundamental system of equations. For the sake of conciseness, the global static system (28) is reported in matrix form as follows:
(29)KU=Q-,
where Q- is the external boundary and domain forces and K is the global stiffness matrix which can be inverted leading to the problem solution. The interested reader can refer to the paper by Viola et al. [39] for further details about static problem resolutions of plates and shells.


### 3.3. The Eigenproblem

Given a generic set of equations of motion, the linear eigenvalue problem can be achieved eliminating the external forces, as it is shown in the following:(30)ℓ(u(x,t))=ρu¨(x,t),
where the right-hand term represents the inertia term, given by the product of density ρ and the acceleration term u¨(x,t). It is noted that the symbol ¨ stands for the second derivative respect to time of a given variable. In general, an eigenproblem seeks particular numbers, called eigenvalues (λ), and certain vectors, called eigenvectors U(x). Mathematically speaking, the eigenvalue problem is obtained from (30) setting as follows:
(31)u(x,t)=U(x)eiλt.
Equation (31) indicates that variable separation technique is employed. In fact, the problem degree of freedom u which is function of space and time is given by the product of two functions, which are functions of space and time separately. Hence, substituting relation (31) into (30), a generalized eigenvalue problem can be achieved. For the sake of completeness, it is presented with its boundary conditions as
(32)ℓ(U(x))+ρλ2U(x)=0inΩ,ℓBU(x)=0on∂Ω.
Applying GDQ rule, the eigenvalue problem (32) can be written in discretised form as follows:
(33)ℒU+λ2MU=0,ℒBU=0,
where M indicates the mass matrix of the structural system under study. Solving the discretised system (33), the eigenvalues λ and eigenvectors U of the given system can be obtained. Finally, the reader can refer to [38] for solutions of eigenvalue problems regarding plates and shells.


## 3.1. Grid Distributions

As it has been widely proven [34–37], the GDQ numerical accuracy strongly depends on the grid distribution choice for a certain numerical problem. On the one hand, the natural and simplest choice of uniform grids does not lead to stable and accurate results for any GDQ computation; on the other hand, the Chebyshev-Gauss-Lobatto (C-G-L) grid point distribution gives accurate results in most of cases. In the numerical results proposed in this paper a C-G-L grid distribution is considered for both directions x and y. Consider
(26)xi=12(1-cos(i-1N-1)π),yj=12(1-cos(j-1M-1)π),
where N, M are the total number of points in the given directions x, y, respectively.


## 3.2. The Static Problem

The partial differential system of equations developed in Section2 has been discretised and numerically solved using the GDQ method. Since GDQ allows turning any derivative in an algebraic expression, the fundamental system of equations is solved in its strong form directly, within its boundary conditions. A static boundary problem can be indicated as follows:
(27)ℓ(u(x))+q(x)=0inΩ,ℓB(u(x))=gon∂Ω,
where ℓ is a partial differential operator and contains all the spatial derivatives of the problem with respect to a Cartesian coordinate system x. The model degrees of freedom are u(x), and the external forces acting on the domain Ω are indicated by q(x). It is obvious that the differential problem cannot be solved without its boundary conditions, where ℓB is the differential operator related to the boundary conditions of the given system of equations and g are the boundary loads.Once the boundary-value problem (27) is discretised by GDQ the following set of algebraic equations is obtained:
(28)ℒU+Q=0,ℒBU=G,
where U is a vector representing a set of unknown functional values at all the interior points, ℒ is a matrix carried out from the GDQ method, and Q, G are the known vectors arising from the functions q(x) and g.In order to obtain the global stiffness matrix of the physical problem under study, the discretised boundary conditions must be embedded into the discrete fundamental system of equations. For the sake of conciseness, the global static system (28) is reported in matrix form as follows:
(29)KU=Q-,
where Q- is the external boundary and domain forces and K is the global stiffness matrix which can be inverted leading to the problem solution. The interested reader can refer to the paper by Viola et al. [39] for further details about static problem resolutions of plates and shells.


## 3.3. The Eigenproblem

Given a generic set of equations of motion, the linear eigenvalue problem can be achieved eliminating the external forces, as it is shown in the following:(30)ℓ(u(x,t))=ρu¨(x,t),
where the right-hand term represents the inertia term, given by the product of density ρ and the acceleration term u¨(x,t). It is noted that the symbol ¨ stands for the second derivative respect to time of a given variable. In general, an eigenproblem seeks particular numbers, called eigenvalues (λ), and certain vectors, called eigenvectors U(x). Mathematically speaking, the eigenvalue problem is obtained from (30) setting as follows:
(31)u(x,t)=U(x)eiλt.
Equation (31) indicates that variable separation technique is employed. In fact, the problem degree of freedom u which is function of space and time is given by the product of two functions, which are functions of space and time separately. Hence, substituting relation (31) into (30), a generalized eigenvalue problem can be achieved. For the sake of completeness, it is presented with its boundary conditions as
(32)ℓ(U(x))+ρλ2U(x)=0inΩ,ℓBU(x)=0on∂Ω.
Applying GDQ rule, the eigenvalue problem (32) can be written in discretised form as follows:
(33)ℒU+λ2MU=0,ℒBU=0,
where M indicates the mass matrix of the structural system under study. Solving the discretised system (33), the eigenvalues λ and eigenvectors U of the given system can be obtained. Finally, the reader can refer to [38] for solutions of eigenvalue problems regarding plates and shells.


## 4. Numerical examples

### 4.1. Static Problems of Cross-Ply Laminated Plates

A simply supported square laminated plate of sidea and thickness h is composed of four equal layers oriented at [0°/90°/90°/0°]. The plate is subjected to a sinusoidal vertical pressure of the form pz=Psin(πx/a)sin(πy/a) with the origin of the coordinate system located at the lower left corner on the midplane and P the maximum load (at the plate center).The orthotropic material properties are given byE1=25.0E2, G12=G13=0.5E2, G23=0.2E2, and ν12=0.25.The in-plane displacements, the transverse displacements, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(34)w¯=102w(a/2,a/2,0)h3E2Pa4,σ¯xx=σxx(a/2,a/2,h/2)h2Pa2,σ¯yy=σyy(a/2,a/2,h/4)h2Pa2,τ¯xz=τxz(0,a/2,0)hPa.In Table1 for the present layer-wise theory, using 7×7 up to 19×19 grid points, the results with higher-order solutions by Akhras et al. [41] and Reddy [11], FSDT solutions by Reddy and Chao [42], and an exact solution by Pagano [43] are compared. Furthermore, the numerical computations by the authors using radial basis functions (RBFs) collocation, by the authors with Reddy’s theory [44], and a layer-wise theory [45] are compared too. The present layer-wise displacements are in excellent agreement for thinner or thicker plates. Highly accurate normal stresses and transverse shear stresses are also obtained. It should be noted that the transverse shear stresses are directly obtained by the constitutive relations, not from the equations of equilibrium, so no recovery procedure has been employed. The thickness-stretching approach is important for thicker plates, as our approach produces the best results for transverse displacements, without shear correction.Table 1
[
0
°
/
9
0
°
/
9
0
°
/
0
°
] square laminated plate under sinusoidal load, applied at (z=h/2).
a
/
h
Method
w
¯
σ
¯
x
x
σ
¯
y
y
τ
¯
z
x
4
HSDT finite strip method [41]
1.8939
0.6806
0.6463
0.2109
HSDT [11]
1.8937
0.6651
0.6322
0.2064
FSDT [42]
1.7100
0.4059
0.5765
0.1398
Elasticity [43]
1.954
0.720
0.666
0.270
Ferreira et al. [44] (N=21)
1.8864
0.6659
0.6313
0.1352
Ferreira (layer-wise) [45] (N=21)
1.9075
0.6432
0.6228
0.2166
Present (7×7 grid)
1.9087
0.6423
0.6258
0.2173
Present (11×11 grid)
1.9091
0.6429
0.6265
0.2173
Present (15×15 grid)
1.9091
0.6429
0.6265
0.2173
Present (19×19 grid)
1.9091
0.6429
0.6265
0.2173
10
HSDT finite strip method [41]
0.7149
0.5589
0.3974
0.2697
HSDT [11]
0.7147
0.5456
0.3888
0.2640
FSDT [42]
0.6628
0.4989
0.3615
0.1667
Elasticity [43]
0.743
0.559
0.403
0.301
Ferreira et al. [44] (N=21)
0.7153
0.5466
0.4383
0.3347
Ferreira (layer-wise) [45] (N=21)
0.7309
0.5496
0.3956
0.2888
Present (7×7 grid)
0.7300
0.5481
0.3963
0.2993
Present (11×11 grid)
0.7303
0.5487
0.3966
0.2993
Present (15×15 grid)
0.7303
0.5487
0.3966
0.2993
Present (19×19 grid)
0.7303
0.5487
0.3966
0.2993
100
HSDT finite strip method [41]
0.4343
0.5507
0.2769
0.2948
HSDT [11]
0.4343
0.5387
0.2708
0.2897
FSDT [42]
0.4337
0.5382
0.2705
0.1780
Elasticity [43]
0.4347
0.539
0.271
0.339
Ferreira et al. [44] (N=21)
0.4365
0.5413
0.3359
0.4106
Ferreira (layer-wise) [45] (N=21)
0.4374
0.5420
0.2697
0.3232
Present (7×7 grid)
0.4345
0.5385
0.2708
0.3359
Present (11×11 grid)
0.4348
0.5391
0.2711
0.3359
Present (15×15 grid)
0.4348
0.5391
0.2711
0.3359
Present (19×19 grid)
0.4348
0.5391
0.2711
0.3359


### 4.2. Three-Layer Square Sandwich Plate in Bending, under Uniform Load

A simply supported sandwich square plate, under a uniform transverse load, is considered. This is the classical sandwich example of Srinivas [46]. The material properties of the sandwich core are expressed in the stiffness matrix, Q¯core, as
(35)Q¯core=[0.9997810.2311920000.2311920.524886000000.262931000000.266810000000.159914].
Skins material properties are related to core properties by a factor R as Q¯skin=RQ¯core.Transverse displacement and stresses are normalized through factors(36)w¯=w(a2,a2,0)0.999781hq,σ¯x1=σx(1)(a/2,a/2,-h/2)q,σ¯x2=σx(1)(a/2,a/2,-2h/5)q,σ¯x3=σx(2)(a/2,a/2,-2h/5)q,σ¯y1=σy(1)(a/2,a/2,-h/2)q,σ¯y2=σy(1)(a/2,a/2,-2h/5)q,σ¯y3=σy(2)(a/2,a/2,-2h/5)q,τ¯xz1=τxz(2)(0,a/2,0)q,τ¯xz2=τxz(2)(0,a/2,-2h/5)q.
Transverse displacement and stresses for a sandwich plate are indicated in Tables 2, 3, and 4 and compared with various sources, for three R ratios. This present approach results are in excellent agreement both in transverse displacement and transverse shear stresses. In particular for higher values of R, the present formulation computes more accurate results than those of Pandya and kant [47]. Some other sources such as laminated shell finite elements [48], multiquadrics [21], and third-order formulation presented by Ferreira et al. [44] are compared. In all cases, the lesson learned is that layer-wise approach should be used instead of ESL theories in the analysis of sandwich soft-core plates, where the skin properties are much higher than core properties. The present technique shows that the thickness-stretching approach is important for thicker plates.Table 2
Square laminated plate under uniform load forR=5.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
5
HSDT [47]
256.13
62.38
46.91
9.382
3.089
2.566
FSDT [47]
236.10
61.87
49.50
9.899
3.313
2.444
CLT
216.94
61.141
48.623
9.783
4.5899
3.386
Ferreira and Barbosa [48]
258.74
59.21
45.61
9.122
3.593
3.593
Ferreira (N=15) [21]
257.38
58.725
46.980
9.396
3.848
2.839
Exact [46]
258.97
60.353
46.623
9.340
4.3641
3.2675
HSDT [44] (N=11)
253.6710
59.6447
46.4292
9.2858
3.8449
1.9650
HSDT [44] (N=15)
256.2387
60.1834
46.8581
9.3716
4.2768
2.2227
HSDT [44] (N=21)
257.1100
60.3660
47.0028
9.4006
4.5481
2.3910
Present (7×7 grid)
257.3663
59.8307
46.4128
9.2826
3.8840
2.6959
Present (11×11 grid)
258.1305
60.0447
46.4237
9.2847
4.0603
2.1679
Present (15×15 grid)
258.1720
60.0772
46.4008
9.2802
4.0916
2.0593
Present (19×19 grid)
258.1794
60.0889
46.3925
9.2785
4.0979
2.0386Table 3
Square laminated plate under uniform load forR=10.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
10
HSDT [47]
152.33
64.65
51.31
5.131
3.147
2.587
FSDT [47]
131.095
67.80
54.24
4.424
3.152
2.676
CLT
118.87
65.332
48.857
5.356
4.3666
3.7075
Ferreira and Barbosa [48]
159.402
64.16
47.72
4.772
3.518
3.518
Ferreira (N=15) [21]
158.55
62.723
50.16
5.01
3.596
3.053
Exact [46]
159.38
65.332
48.857
4.903
4.0959
3.5154
Third-order [44] (N=11)
153.0084
64.7415
49.4716
4.9472
2.7780
1.8207
Third-order [44] (N=15)
154.2490
65.2223
49.8488
4.9849
3.1925
2.1360
Third-order [44] (N=21)
154.6581
65.3809
49.9729
4.9973
3.5280
2.3984
Present (7×7 grid)
158.1023
64.5969
48.7204
4.8720
3.7227
3.1570
Present (11×11 grid)
158.8400
64.9304
48.6523
4.8652
3.9394
2.4751
Present (15×15 grid)
158.8945
64.9907
48.6005
4.8601
3.9772
2.3470
Present (19×19 grid)
158.9036
65.0106
48.5834
4.8583
3.9841
2.3263Table 4
Square laminated plate under uniform load forR=15.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
15
HSDT [47]
110.43
66.62
51.97
3.465
3.035
2.691
FSDT [47]
90.85
70.04
56.03
3.753
3.091
2.764
CLT
81.768
69.135
55.308
3.687
4.2825
3.8287
Ferreira and Barbosa [48]
121.821
65.650
47.09
3.140
3.466
3.466
Ferreira (N=15) [21]
121.184
63.214
50.571
3.371
3.466
3.099
Exact [46]
121.72
66.787
48.299
3.238
3.9638
3.5768
Third-order [44] (N=11)
113.5941
66.3646
49.8957
3.3264
2.1686
1.5578
Third-order [44] (N=15)
114.3874
66.7830
50.2175
3.3478
2.6115
1.9271
Third-order [44] (N=21)
114.6442
66.9196
50.3230
3.3549
3.0213
2.2750
Present (7×7 grid)
120.5292
65.8690
48.2620
3.2175
3.6285
3.3561
Present (11×11 grid)
121.2654
66.3110
48.1025
3.2068
3.8616
2.6168
Present (15×15 grid)
121.3247
66.3973
48.0221
3.2015
3.9001
2.4901
Present (19×19 grid)
121.3340
66.4233
47.9984
3.1999
3.9063
2.4731


### 4.3. Free Vibration Problems of Cross-Ply Laminated Plates

In this example, all layers of the laminate are assumed to be of the same thickness and density and made of the same linearly elastic composite material. The following material parameters of a layer are used:E1/E2=10, 20, 30, or 40; G12=G13=0.6E2; G3=0.5E2; ν12=0.25. The subscripts 1 and 2 denote the directions normal and transverse to the fiber direction in a lamina, which may be oriented at an angle to the plate axes. The ply angle of each layer is measured from the global x-axis to the fiber direction.The example under consideration is a simply supported square plate with a cross-ply lamination scheme[0°/90°/90°/0°]. The thickness and length of the plate are denoted by h and a, respectively. The thickness-to-span ratio h/a=0.2 is employed in the computation. Table 5 lists the fundamental frequency of the simply supported laminate made of various modulus ratios of E1/E2. Figure 2 illustrates the modes of vibration for E1/E2=40 with 15×15 grid points. It is found that the present meshless results are in very close agreement with the values of [49] and the meshfree results of Liew et al. [50] based on the FSDT. The small differences may be due to the consideration of the through-the-thickness deformations in the present formulation.Table 5
The normalized fundamental frequency of the simply supported cross-ply laminated square plate[0°/90°/90°/0°] (w-=(wa2/h)ρ/E2,h/a=0.2).
Method
Grid
E
1
/
E
2
10
20
30
40
Liew et al. [50]
8.2924
9.5613
10.320
10.849
Exact (Khdeir and Librescu) [49]
8.2982
9.5671
10.326
10.854
Present
11
×
11
8.5846
9.8384
10.5695
11.0649
15
×
15
8.5846
9.8384
10.5696
10.0649
19
×
19
8.5846
9.8384
10.5695
10.0649First 4 modes of vibration of four-layer[0°/90°/90°/0°] simply supported laminated plate (w-=(wa2/h)ρ/E2,h/a=0.2), E1/E2=30, grid 19×19 points.
(a)
(b)
(c)
(d)


## 4.1. Static Problems of Cross-Ply Laminated Plates

A simply supported square laminated plate of sidea and thickness h is composed of four equal layers oriented at [0°/90°/90°/0°]. The plate is subjected to a sinusoidal vertical pressure of the form pz=Psin(πx/a)sin(πy/a) with the origin of the coordinate system located at the lower left corner on the midplane and P the maximum load (at the plate center).The orthotropic material properties are given byE1=25.0E2, G12=G13=0.5E2, G23=0.2E2, and ν12=0.25.The in-plane displacements, the transverse displacements, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(34)w¯=102w(a/2,a/2,0)h3E2Pa4,σ¯xx=σxx(a/2,a/2,h/2)h2Pa2,σ¯yy=σyy(a/2,a/2,h/4)h2Pa2,τ¯xz=τxz(0,a/2,0)hPa.In Table1 for the present layer-wise theory, using 7×7 up to 19×19 grid points, the results with higher-order solutions by Akhras et al. [41] and Reddy [11], FSDT solutions by Reddy and Chao [42], and an exact solution by Pagano [43] are compared. Furthermore, the numerical computations by the authors using radial basis functions (RBFs) collocation, by the authors with Reddy’s theory [44], and a layer-wise theory [45] are compared too. The present layer-wise displacements are in excellent agreement for thinner or thicker plates. Highly accurate normal stresses and transverse shear stresses are also obtained. It should be noted that the transverse shear stresses are directly obtained by the constitutive relations, not from the equations of equilibrium, so no recovery procedure has been employed. The thickness-stretching approach is important for thicker plates, as our approach produces the best results for transverse displacements, without shear correction.Table 1
[
0
°
/
9
0
°
/
9
0
°
/
0
°
] square laminated plate under sinusoidal load, applied at (z=h/2).
a
/
h
Method
w
¯
σ
¯
x
x
σ
¯
y
y
τ
¯
z
x
4
HSDT finite strip method [41]
1.8939
0.6806
0.6463
0.2109
HSDT [11]
1.8937
0.6651
0.6322
0.2064
FSDT [42]
1.7100
0.4059
0.5765
0.1398
Elasticity [43]
1.954
0.720
0.666
0.270
Ferreira et al. [44] (N=21)
1.8864
0.6659
0.6313
0.1352
Ferreira (layer-wise) [45] (N=21)
1.9075
0.6432
0.6228
0.2166
Present (7×7 grid)
1.9087
0.6423
0.6258
0.2173
Present (11×11 grid)
1.9091
0.6429
0.6265
0.2173
Present (15×15 grid)
1.9091
0.6429
0.6265
0.2173
Present (19×19 grid)
1.9091
0.6429
0.6265
0.2173
10
HSDT finite strip method [41]
0.7149
0.5589
0.3974
0.2697
HSDT [11]
0.7147
0.5456
0.3888
0.2640
FSDT [42]
0.6628
0.4989
0.3615
0.1667
Elasticity [43]
0.743
0.559
0.403
0.301
Ferreira et al. [44] (N=21)
0.7153
0.5466
0.4383
0.3347
Ferreira (layer-wise) [45] (N=21)
0.7309
0.5496
0.3956
0.2888
Present (7×7 grid)
0.7300
0.5481
0.3963
0.2993
Present (11×11 grid)
0.7303
0.5487
0.3966
0.2993
Present (15×15 grid)
0.7303
0.5487
0.3966
0.2993
Present (19×19 grid)
0.7303
0.5487
0.3966
0.2993
100
HSDT finite strip method [41]
0.4343
0.5507
0.2769
0.2948
HSDT [11]
0.4343
0.5387
0.2708
0.2897
FSDT [42]
0.4337
0.5382
0.2705
0.1780
Elasticity [43]
0.4347
0.539
0.271
0.339
Ferreira et al. [44] (N=21)
0.4365
0.5413
0.3359
0.4106
Ferreira (layer-wise) [45] (N=21)
0.4374
0.5420
0.2697
0.3232
Present (7×7 grid)
0.4345
0.5385
0.2708
0.3359
Present (11×11 grid)
0.4348
0.5391
0.2711
0.3359
Present (15×15 grid)
0.4348
0.5391
0.2711
0.3359
Present (19×19 grid)
0.4348
0.5391
0.2711
0.3359


## 4.2. Three-Layer Square Sandwich Plate in Bending, under Uniform Load

A simply supported sandwich square plate, under a uniform transverse load, is considered. This is the classical sandwich example of Srinivas [46]. The material properties of the sandwich core are expressed in the stiffness matrix, Q¯core, as
(35)Q¯core=[0.9997810.2311920000.2311920.524886000000.262931000000.266810000000.159914].
Skins material properties are related to core properties by a factor R as Q¯skin=RQ¯core.Transverse displacement and stresses are normalized through factors(36)w¯=w(a2,a2,0)0.999781hq,σ¯x1=σx(1)(a/2,a/2,-h/2)q,σ¯x2=σx(1)(a/2,a/2,-2h/5)q,σ¯x3=σx(2)(a/2,a/2,-2h/5)q,σ¯y1=σy(1)(a/2,a/2,-h/2)q,σ¯y2=σy(1)(a/2,a/2,-2h/5)q,σ¯y3=σy(2)(a/2,a/2,-2h/5)q,τ¯xz1=τxz(2)(0,a/2,0)q,τ¯xz2=τxz(2)(0,a/2,-2h/5)q.
Transverse displacement and stresses for a sandwich plate are indicated in Tables 2, 3, and 4 and compared with various sources, for three R ratios. This present approach results are in excellent agreement both in transverse displacement and transverse shear stresses. In particular for higher values of R, the present formulation computes more accurate results than those of Pandya and kant [47]. Some other sources such as laminated shell finite elements [48], multiquadrics [21], and third-order formulation presented by Ferreira et al. [44] are compared. In all cases, the lesson learned is that layer-wise approach should be used instead of ESL theories in the analysis of sandwich soft-core plates, where the skin properties are much higher than core properties. The present technique shows that the thickness-stretching approach is important for thicker plates.Table 2
Square laminated plate under uniform load forR=5.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
5
HSDT [47]
256.13
62.38
46.91
9.382
3.089
2.566
FSDT [47]
236.10
61.87
49.50
9.899
3.313
2.444
CLT
216.94
61.141
48.623
9.783
4.5899
3.386
Ferreira and Barbosa [48]
258.74
59.21
45.61
9.122
3.593
3.593
Ferreira (N=15) [21]
257.38
58.725
46.980
9.396
3.848
2.839
Exact [46]
258.97
60.353
46.623
9.340
4.3641
3.2675
HSDT [44] (N=11)
253.6710
59.6447
46.4292
9.2858
3.8449
1.9650
HSDT [44] (N=15)
256.2387
60.1834
46.8581
9.3716
4.2768
2.2227
HSDT [44] (N=21)
257.1100
60.3660
47.0028
9.4006
4.5481
2.3910
Present (7×7 grid)
257.3663
59.8307
46.4128
9.2826
3.8840
2.6959
Present (11×11 grid)
258.1305
60.0447
46.4237
9.2847
4.0603
2.1679
Present (15×15 grid)
258.1720
60.0772
46.4008
9.2802
4.0916
2.0593
Present (19×19 grid)
258.1794
60.0889
46.3925
9.2785
4.0979
2.0386Table 3
Square laminated plate under uniform load forR=10.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
10
HSDT [47]
152.33
64.65
51.31
5.131
3.147
2.587
FSDT [47]
131.095
67.80
54.24
4.424
3.152
2.676
CLT
118.87
65.332
48.857
5.356
4.3666
3.7075
Ferreira and Barbosa [48]
159.402
64.16
47.72
4.772
3.518
3.518
Ferreira (N=15) [21]
158.55
62.723
50.16
5.01
3.596
3.053
Exact [46]
159.38
65.332
48.857
4.903
4.0959
3.5154
Third-order [44] (N=11)
153.0084
64.7415
49.4716
4.9472
2.7780
1.8207
Third-order [44] (N=15)
154.2490
65.2223
49.8488
4.9849
3.1925
2.1360
Third-order [44] (N=21)
154.6581
65.3809
49.9729
4.9973
3.5280
2.3984
Present (7×7 grid)
158.1023
64.5969
48.7204
4.8720
3.7227
3.1570
Present (11×11 grid)
158.8400
64.9304
48.6523
4.8652
3.9394
2.4751
Present (15×15 grid)
158.8945
64.9907
48.6005
4.8601
3.9772
2.3470
Present (19×19 grid)
158.9036
65.0106
48.5834
4.8583
3.9841
2.3263Table 4
Square laminated plate under uniform load forR=15.
R
Method
w
¯
σ
¯
x
1
σ
¯
x
2
σ
¯
x
3
τ
¯
x
z
1
τ
¯
x
z
2
15
HSDT [47]
110.43
66.62
51.97
3.465
3.035
2.691
FSDT [47]
90.85
70.04
56.03
3.753
3.091
2.764
CLT
81.768
69.135
55.308
3.687
4.2825
3.8287
Ferreira and Barbosa [48]
121.821
65.650
47.09
3.140
3.466
3.466
Ferreira (N=15) [21]
121.184
63.214
50.571
3.371
3.466
3.099
Exact [46]
121.72
66.787
48.299
3.238
3.9638
3.5768
Third-order [44] (N=11)
113.5941
66.3646
49.8957
3.3264
2.1686
1.5578
Third-order [44] (N=15)
114.3874
66.7830
50.2175
3.3478
2.6115
1.9271
Third-order [44] (N=21)
114.6442
66.9196
50.3230
3.3549
3.0213
2.2750
Present (7×7 grid)
120.5292
65.8690
48.2620
3.2175
3.6285
3.3561
Present (11×11 grid)
121.2654
66.3110
48.1025
3.2068
3.8616
2.6168
Present (15×15 grid)
121.3247
66.3973
48.0221
3.2015
3.9001
2.4901
Present (19×19 grid)
121.3340
66.4233
47.9984
3.1999
3.9063
2.4731


## 4.3. Free Vibration Problems of Cross-Ply Laminated Plates

In this example, all layers of the laminate are assumed to be of the same thickness and density and made of the same linearly elastic composite material. The following material parameters of a layer are used:E1/E2=10, 20, 30, or 40; G12=G13=0.6E2; G3=0.5E2; ν12=0.25. The subscripts 1 and 2 denote the directions normal and transverse to the fiber direction in a lamina, which may be oriented at an angle to the plate axes. The ply angle of each layer is measured from the global x-axis to the fiber direction.The example under consideration is a simply supported square plate with a cross-ply lamination scheme[0°/90°/90°/0°]. The thickness and length of the plate are denoted by h and a, respectively. The thickness-to-span ratio h/a=0.2 is employed in the computation. Table 5 lists the fundamental frequency of the simply supported laminate made of various modulus ratios of E1/E2. Figure 2 illustrates the modes of vibration for E1/E2=40 with 15×15 grid points. It is found that the present meshless results are in very close agreement with the values of [49] and the meshfree results of Liew et al. [50] based on the FSDT. The small differences may be due to the consideration of the through-the-thickness deformations in the present formulation.Table 5
The normalized fundamental frequency of the simply supported cross-ply laminated square plate[0°/90°/90°/0°] (w-=(wa2/h)ρ/E2,h/a=0.2).
Method
Grid
E
1
/
E
2
10
20
30
40
Liew et al. [50]
8.2924
9.5613
10.320
10.849
Exact (Khdeir and Librescu) [49]
8.2982
9.5671
10.326
10.854
Present
11
×
11
8.5846
9.8384
10.5695
11.0649
15
×
15
8.5846
9.8384
10.5696
10.0649
19
×
19
8.5846
9.8384
10.5695
10.0649First 4 modes of vibration of four-layer[0°/90°/90°/0°] simply supported laminated plate (w-=(wa2/h)ρ/E2,h/a=0.2), E1/E2=30, grid 19×19 points.
(a)
(b)
(c)
(d)


## 5.   Conclusions

The first-order and the third-order shear deformation theories are laminate-wise, with laminate degrees of freedom, where all layers have the same rotations. Layer-wise formulations can better accommodate the deformation behaviour of some laminates, in particular for sandwich laminates, where core and skin materials are so different.In this paper the static and free vibration analysis of composite laminated plates by the use of generalized differential quadrature and using a layer-wise theory with independent rotations in each layer is performed.The equations of motion were derived and interpolated; moreover, boundary conditions interpolation was schematically formulated.Finally, composite laminated plate and sandwich plate were considered for testing the present methodology. The obtained results show an excellent accuracy for all cases.



---

*Source: 964367-2013-09-25.xml*