# Boundary Conditions in 2D Numerical and 3D Exact Models for Cylindrical Bending Analysis of Functionally Graded Structures

Authors: F. Tornabene; S. Brischetto; N. Fantuzzi; M. Bacciocchi
Journal: Shock and Vibration
 (2016)

---

## Abstract

The cylindrical bending condition for structural models is very common in the literature because it allows an incisive and simple verification of the proposed plate and shell models. In the present paper, 2D numerical approaches (the Generalized Differential Quadrature (GDQ) and the finite element (FE) methods) are compared with an exact 3D shell solution in the case of free vibrations of functionally graded material (FGM) plates and shells. The first 18 vibration modes carried out through the 3D exact model are compared with the frequencies obtained via the 2D numerical models. All the 18 frequencies obtained via the 3D exact model are computed when the structures have simply supported boundary conditions for all the edges. If the same boundary conditions are used in the 2D numerical models, some modes are missed. Some of these missed modes can be obtained modifying the boundary conditions imposing free edges through the direction perpendicular to the direction of cylindrical bending. However, some modes cannot be calculated via the 2D numerical models even when the boundary conditions are modified because the cylindrical bending requirements cannot be imposed for numerical solutions in the curvilinear edges by definition. These features are investigated in the present paper for different geometries (plates, cylinders, and cylindrical shells), types of FGM law, lamination sequences, and thickness ratios.


---

## Body

## 1. Introduction

The cylindrical bending conditions have a great diffusion in the open literature because they allow an immediate and simple verification of several plate and shell models. Some of these plate models which use the cylindrical bending hypotheses for the verifications are discussed in the following. Chen et al. [1] used the state-space method for the investigation of a simply supported cross-ply laminated plate embedding viscous interfaces in cylindrical bending. Chen and Lee [2, 3] also used the state-space method for the cylindrical bending analysis of simply supported angle-ply laminated plates with interfacial damage and simply supported angle-ply laminated plates subjected to a static load. An experimental three-roller cylindrical bending investigation for plates was proposed by Gandhi and Raval [4]; the comparison with analytical solutions was also conducted. Oral and Darendeliler [5] proposed a methodology for the design of plate-forming dies in cylindrical bending using optimization techniques which allow the cost reductions of die production. Considering the large deflection of a thin beam, under certain conditions, the solution of the plate problem is not unique [6]. The bending-gradient plate theory, seen as an extension of the Reissner-Mindlin plate model, was used in [7] for the cylindrical bending analysis of laminated plates. Nimbolkar and Jain [8] used the same Reissner-Mindlin plate theory for the cylindrical bending investigation of composite and elastic plates subjected to the mechanical transverse load under plain strain conditions. Sayyad and Ghugal [9] investigated laminated plates using a nth order shear deformation theory under cylindrical bending requirements. Some of the most important three-dimensional exact plate solutions were proposed by Pagano [10] and Pagano and Wang [11] for the cylindrical bending analysis of multilayered composite plates. Saeedi et al. [12] developed a 2D plate layer-wise model for the analysis of delamination growth in multilayered plates subjected to cylindrical bending loadings. This method is alternative to the finite element method. Shu and Soldatos [13] developed a new stress analysis method to investigate the stress distributions in angle-ply laminated plates when subjected to cylindrical bending. Further cylindrical bending analyses for plates consider the inclusion of functionally graded material (FGM) layers; typical examples are [14, 15]. In [16–21], the cylindrical bending analysis was proposed for multilayered plates embedding piezoelectric and/or FGM layers. The cylindrical bending analysis for multilayered composite and/or piezoelectric plates proposed in [22–24] includes the application of thermal loads. The cylindrical bending analysis of shells is not so diffused in the literature. Few cases can be found, and they concern the cylindrical bending conditions imposed to cylindrical shells. Simply supported angle-ply laminated cylindrical shells in cylindrical bending were investigated in [25] using the state-space formulation for the bending and free vibration analyses. A three-dimensional piezoelectric model based on the perturbation method for the cylindrical bending electromechanical load applications to cylindrical shells was proposed in [26]. Yan et al. [27] analyzed the cylindrical bending behavior of a simply supported angle-ply laminated cylindrical shell embedding viscoelastic interfaces.The authors proposed the free vibration analysis of simply supported FGM plates and shells in [28, 29] where the 3D exact shell model was compared with several 2D numerical models such as the classical finite element (FE) one and the classical and refined Generalized Differential Quadrature (GDQ) methods. Similar comparisons were also proposed in [30, 31] for one-layered and multilayered isotropic, composite, and sandwich plates and shells. In such analyses, the considered geometries are plates, cylinders, and cylindrical and spherical shell panels. Low and high frequencies were calculated for thick and thin simply supported plates and shells with constant radii of curvature. The comparison between the 2D numerical methods (FE and GDQ solutions) and the 3D exact solution is possible only if an exhaustive vibration mode investigation is conducted. This investigation allows understanding how to make the comparison between numerical and analytical methods. This comparison is not easy because the exact 3D solution gives infinite vibration modes for all the possible combinations of half-wave numbers (m,n) and 2D numerical models propose a finite number of vibration modes because they use a finite number of degrees of freedom in the plane and in the thickness direction. It should be underlined that the 3D exact model is able to investigate any kind of geometric ratio (from extremely thick to very thin structures); nevertheless it can model only simply supported boundary conditions and it can investigate structures wherein constant radii of curvature are present. The 2D numerical approaches are more versatile (in terms of boundary conditions and type of structure analyzed); however they can be used only for moderately thick structures. It is obvious that thick structures can be analyzed only when higher order displacement fields are set in such 2D models. Since commercial codes implement only through the thickness linear theories, they can investigate accurately only moderately thick shells. The method proposed in [28–31] to compare the 3D exact shell model and the 2D numerical models is based on the calculation of the frequencies via the 2D numerical code and their appropriate visualization of the relative vibration modes via the GDQ and/or FE method. After this visualization, the evaluation of the 3D exact frequencies is possible by means of the appropriate imposition of half-wave numbers obtained from the vibration mode study. From this procedure, it is clear how the 3D analysis could give some frequencies that cannot be calculated by the 2D numerical codes. This feature was not the aim of the papers [28–31], and for this reason these missed frequencies are investigated in the present new work for FGM structures.This new analysis is conducted calculating the first 18 frequencies for one-layered and sandwich FGM plates and shells by means of the 3D exact model. The half-wave numbers are combined along the two directionsα and β in the plane considering the simply supported boundary conditions for the four edges. In these first 18 frequencies there are some modes that cannot be calculated by means of the 2D numerical models because they are not solutions of the employed mathematical model with the four simply supported edges. In general, these missed modes by the 2D numerical models correspond to the cylindrical bending conditions. These conditions in the 3D exact solution mean that one of the two half-wave numbers is equal to zero through an edge direction with a transverse displacement different from zero along the perpendicular edge. These particular frequencies are not solutions of the 2D numerical models built with the four edges which are simply supported, but a different mathematical model must be defined where two edges are simply supported and the other two are free. These conditions allow the cylindrical bending, with the transverse displacement through two parallel edges different from zero, in numerical models in order to obtain the missed frequencies and modes. This condition is possible only when the zero half-wave number is imposed through a rectilinear side because the half-wave number equals zero in the curvilinear edge does not mean that all the derivatives in the that direction are zero because of the presence of the curvature. This last condition is in conflict with the definition of cylindrical bending.The 3D exact shell model is described and validated in Section2. It uses the equilibrium equations written using the generic curvilinear orthogonal coordinates. Such equations are solved in closed form by means of simply supported boundary conditions. The differential equations along z, written in layer-wise form, are solved using the exponential matrix method (also known as transfer matrix method or the state-space approach). The 2D FE and GDQ models are proposed and opportunely validated in Section 3. The 2D FE model makes use of a very common commercial finite element code. The 2D GDQ models give results by means of an in-house academic software and they are based on classical and refined shell models. The results about the comparisons between 3D exact and 2D numerical models, with the remark about the boundary conditions, are shown Section 4. Section 5 gives the main conclusions.


## 2. Exact Solution of a 3D Shell Model

The three-dimensional exact shell model employed for the comparisons proposed in Section4 has been elaborated in [32–40] for the free frequency analyses of single-layered and multilayered isotropic, orthotropic, composite, FGM, and sandwich plates and shells with constant radii of curvature and for the free vibrations of single- and double-walled carbon nanotubes including van der Waals interactions. The equilibrium equations proposed in the general orthogonal curvilinear coordinate system (α,β,z) are valid for both plates and shells with constant radii of curvature. For the sake of clarity, the plates and shells geometries within their local coordinate systems are depicted in Figure 1. The exact form of these equations is obtained using simply supported boundary conditions, while the differential equations in z are solved by means of the exponential matrix method (also known as transfer matrix method or the state-space approach).Figure 1
Geometries of plates and shells with constant radii of curvature and their local reference system.
(a) (b) (c) (d)Considering spherical shells with constant radii of curvatureRα and Rβ and embedding NL classical or FGM layers, the differential equations of equilibrium for the free vibration analysis written for a generic k layer are(1)Hβ∂σααk∂α+Hα∂σαβk∂β+HαHβ∂σαzk∂z+2HβRα+HαRβσαzk=ρkHαHβu¨k,Hβ∂σαβk∂α+Hα∂σββk∂β+HαHβ∂σβzk∂z+2HαRβ+HβRασβzk=ρkHαHβv¨k,Hβ∂σαzk∂α+Hα∂σβzk∂β+HαHβ∂σzzk∂z-HβRασααk-HαRβσββk+HβRα+HαRβσzzk=ρkHαHβw¨k.It is remarked that this set of equations is valid only for spherical shells whereinRα=Rβ=R, because for doubly curved shells with variable radii of curvature such equations take a different form. The authors decided to keep such system with the two radii of curvature distinct because it will be much easier to understand the cylindrical bending effect for plates and shells starting from a general point of view.In the case of shell geometries and/or FGM layers, each physical layerk must be divided in an opportune number of q mathematical layers to approximate the shell curvature and/or the FGM law through the thickness direction z. The index j indicates the total mathematical layers obtained from the subdivision of each k layer in q mathematical layers (j=k×q). The parametric coefficients for shells with constant radii of curvature can be written as(2)Hα=1+zRα,Hβ=1+zRβ,Hz=1.A closed form solution for simply supported shells and plates can be developed by writing harmonic displacements as(3)ujα,β,z,t=Ujzeiωtcos⁡α¯αsin⁡β¯β,vjα,β,z,t=Vjzeiωtsin⁡α¯αcos⁡β¯β,wjα,β,z,t=Wjzeiωtsin⁡α¯αsin⁡β¯β,where Uj(z), Vj(z), and Wj(z) indicate the displacement amplitudes in α, β, and z directions, respectively. i indicates the imaginary unit. ω=2πf is the circular frequency where f is the frequency value in Hz. t is the symbol for the time. In the coefficients α¯=mπ/a and β¯=nπ/b, m and n are the half-wave numbers and a and b are the shell dimensions in α and β directions, respectively (all these quantities are evaluated on the midsurface Ω0). ρj is the mass density. (σααj,σββj,σzzj,σβzj,σαzj,σαβj) are the six stress components and u¨j, v¨j, and w¨j are the second temporal derivative of the three displacement components. Rα and Rβ indicate the radii of curvature defined with respect the reference surface Ω0 of the whole multilayered plate or shell. Hα and Hβ proposed in (2) vary with continuity through the thickness direction of the multilayered structure. Therefore, they depend on the thickness coordinate z. The reference surface Ω0 of the structure is defined as the locus of points which are positioned midway between the external surfaces. Displacement components are indicated as u, v, and w in α, β, and z directions, respectively.The final system in closed form is developed substituting (3) and constitutive and geometrical equations (which are not here given for the sake of brevity but they can be found in detail in [32]) in (1):(4)-C55jHβjHαjRα2-C55jRαRβ-α-2C11jHβjHαj-β-2C66jHαjHβj+ρjHαjHβjω2Uj+-α-β-C12j-α-β-C66jVj+α-C11jHβjHαjRα+α-C12jRβ+α-C55jHβjHαjRα+α-C55jRβWj+C55jHβjRα+C55jHαjRβU,zj+α-C13jHβj+α-C55jHβjW,zj+C55jHαjHβjU,zzj=0,-α-β-C66j-α-β-C12jUj+-C44jHαjHβjRβ2-C44jRαRβ-α-2C66jHβjHαj-β-2C22jHαjHβj+ρjHαjHβjω2Vj+β-C44jHαjHβjRβ+β-C44jRα+β-C22jHαjHβjRβ+β-C12jRαWj+C44jHαjRβ+C44jHβjRαV,zj+β-C44jHαj+β-C23jHαjW,zj+C44jHαjHβjV,zzj=0,α¯C55jHβjHαjRα-α¯C13jRβ+α¯C11jHβjHαjRα+α¯C12jRβUj+β¯C44jHαjHβjRβ-β¯C23jRα+β¯C22jHαjHβjRβ+β¯C12jRαVj+C13jRαRβ+C23jRαRβ-C11jHβjHαjRα2-2C12jRαRβ-C22jHαjHβjRβ2-α¯2C55jHβjHαj-β¯2C44jHαjHβj+ρjHαjHβjω2Wj+-α¯C55jHβj-α¯C13jHβjU,zj+-β¯C44jHαj-β¯C23jHαjV,zj+C33jHβjRα+C33jHαjRβW,zj+C33jHαjHβjW,zzj=0. The cylindrical bending requirements for analytical methods are obtained choosing one of the two half-wave numbers which equals zero. This condition means α¯=mπ/a=0 or β¯=nπ/b=0 which allows the derivatives made in a particular direction (α or β) which equals zero too. This mathematical condition means that each section in the cylindrical bending direction has always the same behavior.It is obvious that the cylindrical bending effect is proper for cylindrical structures only. However, the authors started the present treatise from the governing equations (1) of spherical shell structures because one of the aims of the present contribution is to underline which are the extra terms (related to the curvature) that do not allow having the structure to be in cylindrical bending conditions.The compact form of (4) has been elaborated in detail in [32]. The final version is here proposed:(5)Dj∂Uj∂z=AjUj,where ∂Uj/∂z=Uj′ and Uj=UjVjWjUj′Vj′Wj′. Equation (5) is rewritten as(6)DjUj′=AjUj,(7)Uj′=Dj-1AjUj,(8)Uj′=Aj∗Uj,with Aj∗=Dj-1Aj.The exponential matrix method (also known as transfer matrix method or the state-space approach) is employed for the solution of (8) (details are given in [32]):(9)Ujzj=exp⁡Aj∗zjUj0withzj∈0,hj,where zj is the thickness coordinate for the j mathematical layer with values from 0 at the bottom to hj at the top. The exponential matrix is calculated using zj=hj for each j mathematical layer:(10)Aj∗∗=exp⁡Aj∗hj=I+Aj∗hj+Aj∗22!hj2+Aj∗33!hj3+⋯+Aj∗NN!hjN,where I indicates the 6×6 identity matrix. The expansion used for the exponential matrix in (10) has a fast convergence and it is not time consuming from the computational point of view.j = M mathematical layers are used to approximate the shell curvature and/or the FGM law through the thickness. M-1 transfer matrices must be defined by means of the interlaminar continuity conditions of displacements and transverse shear/normal stresses at each interface. The closed form solution for plates or shells is obtained using simply supported boundary conditions and free stress conditions at the top and at the bottom surfaces. The above conditions allow the following final system:(11)EU10=0,and the 6×6 matrix E does not change, independently of the number of layers M, even if a layer-wise approach is applied. U1(0) is the vector U defined at the bottom of the whole multilayered shell or plate (first layer 1 with z1=0). All the steps and details, here omitted for the sake of brevity, can be found in [32–40] where the proposed exact 3D shell model has been elaborated and tested for the first time for several applications. Several validations for the 3D exact shell model have been shown in [32–40]. For this reason, the model is now used with confidence for the results and comparisons proposed in Section 4.The nontrivial solution ofU1(0) in (11) allows conducting the free vibration analysis imposing the zero determinant of matrix E:(12)det⁡E=0,Equation (12) is used to find the roots of an higher order polynomial written in λ=ω2. The imposition of a couple of (m,n) allows a certain number of circular frequencies ω=2πf (from I to ∞) to be calculated. This number of frequencies depends on the order N used in each exponential matrix Aj∗∗ and the number M of mathematical layers employed for the shell curvature and/or the FGM law approximations. The results shown in [32–40] confirm that N=3 for the exponential matrix and M=100 for the mathematical layers are appropriate to always calculate the correct free frequencies for each structure, lamination sequence, number of layers, material, FGM law, and thickness value.


## 3. 2D Finite Element Analysis and 2D Generalized Differential Quadrature Models

The 3D exact natural frequencies are compared with two numerical 2D solutions given by a weak and a strong formulation. The former is a classic 2D finite element (FE) model carried out through the commercial software Strand7. The latter is a recent widespread technique known as Generalized Differential Quadrature (GDQ) method [41]. Details related to convergence and accuracy of these techniques can be found in the works [29–31]. It is noted that the FE solution employs a classic linear theory, known also as Reissner-Mindlin theory, with a linear rotation of the fibers through the shell thickness. The present GDQ solution exploits higher order theories with a general expansion through the thickness of the structure. This general approach falls within the framework of the unified theories because several structural models can be obtained via the repetition of a fundamental nucleus. Another difference between the GDQ applications and the FE ones is that the former includes the curvature effects and the rotary inertias in its formulation.The 2D kinematic hypothesis [42–50] considers the following expansion of the degrees of freedom:(13)U=∑τ=0Nc+1Fτuτ,where U represents the vector of the 3D displacement components, which is function of (α,β,z), and uτ is the vector of the displacement parameters, function of (α,β), which depends on the free index τ. Therefore uτ lays on the midsurface of the shell [42]. The matrix Fτ(ij)=δijFτ (for i,j=1,2,3) contains the so-called thickness functions where δ is the Kronecker delta function. The thickness functions describe how the displacement parameters characterize the solution through the thickness. They can have any form of any smooth function of z; however the classic form of these functions is the power function. Thus, Fτ=zτ can be 1 (τ=0), z(τ=1), z2(τ=2), z3(τ=3), and so on. Please note that according to the mathematical background of the present Unified Formulation the index Nc+1 is reserved to the (optional) zigzag effect. Since the aim of the present paper is to investigate the cylindrical bending effect for functionally graded plates and cylindrical panels, the following numerical applications will consider a fixed value of Nc=3, so that τ=0,1,2,3. This choice is justified by the fact that third order theories are sufficient, with a reasonably limited number of degrees of freedom, to obtain accurate solutions when compared to 3D shell models [42–50]. For multilayered structures, such as sandwich configurations, the zigzag function must be added into the formulation. Once again, as shown in (13), the zigzag function is embedded at the end of the expansion with the exponent Nc+1. When such effect is neglected the sum of (13) stops at Nc.The same theory used for the FE solution can be carried out in the GDQ method settingτ=0,1 and Nc=1; however it has been shown in other recent works that the strong form approaches led to more accurate results with respect to the FE ones [46, 48].The relation between generalized strainsε(τ) and displacements u(τ) can be set as(14)ετ=DΩuτforτ=0,1,2,…,Nc,Nc+1,where DΩ is a differential operator explicitly defined in [42]. The present constitutive material model is linear and elastic; therefore for a general kth lamina, Hooke’s law can be applied:(15)σk=C¯kεk,where σ(k) is the vector containing the stress components and C¯(k) represents the constitutive matrix for each kth ply [42] of the structure. Each component of the constitutive matrix is indicated by the symbol C¯ξη(k) and it is valid in the curvilinear reference system O′αβz of the reference geometry. These coefficients are evaluated with the application of the transformation matrix [42].The stress resultants are the terms which correspond to the punctual 3D stresses of the 3D doubly curved solid. Their definition can be achieved by developing the Hamilton principle for the present structural theory [42]:(16)Sτ=∑k=1NL∫zkzk+1ZτTστHαHβdzforτ=0,1,2,…,Nc,Nc+1,where Z(τ) is a 6×9 matrix containing the thickness functions Fτ and the curvature coefficients Hα and Hβ; the definition of such matrix with some other details can be found in [42]. S(τ) include the stress resultants as a function of the order τ.The definition of (16) can be rewritten as(17)Sτ=∑s=0Nc+1Aτsεsforτ=0,1,2,…,Nc,Nc+1,which is a function of generalized strains of order s  ε(s).Note that the so-called elastic coefficients are defined as(18)Aτs=∑k=1NL∫zkzk+1ZτTC¯kZsHαHβdz,and they can be carried out by the following integral expressions:(19)Aξηϱζτs=∑k=1NL∫zkzk+1C¯ξηkzszτHαHβHαϱHβζdzAξηϱζτ~s=∑k=1NL∫zkzk+1C¯ξηkzsτzτ-1HαHβHαϱHβζdzAξηϱζτs~=∑k=1NL∫zkzk+1C¯ξηkszs-1zτHαHβHαϱHβζdzAξηϱζτ~s~=∑k=1NL∫zkzk+1C¯ξηkszs-1τzτ-1HαHβHαϱHβζdzforτ,s=0,1,2,…,Nc,Nc+1,ξ,η=1,2,3,4,5,6,ϱ,ζ=0,1,2.It is recalled that the exponent Nc+1 is reserved for the Murakami function; otherwise this order is neglected. The subscripts τ~ and s~ indicate the derivatives of the thickness functions Fτ and Fs with respect to z. The subscripts ϱ and ζ are the exponents of the quantities Hα and Hβ, whereas ξ and η are the indices of the material constants C¯ξη(k) for the kth ply [42–49]. Note that C¯ξη(k), Hα, and Hβ are quantities that depend on the thickness z. Therefore, the integrals of (19) cannot be computed analytically. The present solution for solving such integrals is to employ a numerical integration algorithm. In particular, the so-called Generalized Integral Quadrature (GIQ) method is utilized. For an accurate evaluation of the integrals, 51 points per layer are considered in the integration.The Hamilton principle supplies for the dynamic equilibrium equations and the correspondent boundary conditions. In detail, the fundamental nucleus of three motion equations takes the form(20)DΩ⋆Sτ=∑s=0Nc+1Mτsu¨sforτ=0,1,2,…,Nc,Nc+1,where DΩ⋆ is the equilibrium operator and M(τs) is the inertial matrix. Their explicit expressions have been presented in a previous authors’ paper [42].Combining the kinematic equation (14), constitutive equation (17), and the motion equation (20), the fundamental governing system in terms of displacement parameters is carried out:(21)∑s=0Nc+1Lτsus=∑s=0Nc+1Mτsu¨sforτ=0,1,2,…,Nc,Nc+1,where L(τs)=DΩ⋆A(τs)DΩ is the fundamental operator [42].Boundary conditions must be included for the solution of the differential equations (21). The present paper deals only with free and simply supported boundary conditions. The former is used to reproduce the cylindrical bending conditions considering alternatively one free edge and one simply supported edge. When the four edges are all simply supported, the problem is the same as the ones investigated in the previous papers [29–31]. Simply supported conditions are(22)Nατ=0,uβτ=uzτ=0τ=0,1,2,…,Nc,Nc+1atα=α0orα=α1,β0≤β≤β1,uατ=0,Nβτ=0,uzτ=0τ=0,1,2,…,Nc,Nc+1atβ=β0orβ=β1,α0≤α≤α1,Free edge conditions are(23)Nατ=0,Nβτ=0,Tατ=0τ=0,1,2,…,Nc,Nc+1atα=α0orα=α1,β0≤β≤β1,Nατ=0,Nβτ=0,Tβτ=0τ=0,1,2,…,Nc,Nc+1atβ=β0orβ=β1,α0≤α≤α1.


## 4. Results

The present section includes some numerical applications for the comparison of the 3D exact shell vibrations with two numerical 2D models. In particular classical simply supported solutions and cylindrical bending (CB) behavior of plates and shells are presented. The following computations focus on plates, circular cylinders, and cylindrical panels due to the fact that the CB effect occurs only when at least two parallel straight edges are present in the structure.The first case is a square plate (a=b=1m) with simply supported edges. Two thickness ratios are analyzed: thin a/h=100 and moderately thick a/h=20 structure. The plate is made of a single functionally graded layer (h1=h) with volume fraction of ceramic (c) phase:(24)Vcz=0.5+zhpwithz∈-h2,h2,where p is the power exponent which varies the volume fraction of the ceramic component through the thickness. Two limit cases occur for (24): Vc=1 when p=0 and Vc→0 when p=∞. In the following computations it is set as p=0.5. Note that the volume fraction function in (24) gives a material which is ceramic (c) rich at the top and metal (m) rich at the bottom of the structure. The mechanical properties (Young modulus and mass density) are given by the theory of mixtures as(25)Ez=Em+Ec-EmVc,ρz=ρm+ρc-ρmVc,where Em and ρm are the elastic modulus and mass density of the metallic phase, respectively. Ec and ρc are the corresponding properties of the ceramic phase. Note that in all the following computations the Poisson ratio is considered constant through the thickness; thus νc=νm=ν=0.3. The two constituents have the following mechanical properties: (26)Em=73GPa,ρm=2800kg/m3,Ec=380GPa,ρc=3800kg/m3.Since the FE commercial software does not have an embedded tool for the computation of functionally graded properties through the thickness, the present FE applications are carried out building a fictitious laminated structure with isotropic plies. In detail, the thickness of the FGM lamina has been divided into 100 fictitious layers and for each layer the corresponding mechanical properties have been computed through (25). The numerical model has been made of a uniform element division of 40×40 Quad8 elements for a total number of 29280 degrees of freedom (dofs). The GDQ solution is computed with a not-uniform Chebyshev-Gauss-Lobatto grid [41] with 31×31 points for both thickness ratios. It is underlined that since a single ply is set, no zigzag function is needed in the computation. It is recalled that the present higher order GDQ model considers a third order expansion with τ=0,1,2,3, whereas the FE one has a first order model (also known as Reissner-Mindlin model). Table 1 reports the results in terms of natural frequency (Hz), wherein the first column reports the number of half-waves (m,n) of the 3D exact solution. Note that, for a fixed set of half-wave numbers, different mode shapes can occur (indicated by the second column “Mode”). In particular, I indicates the first mode and II indicates the second mode. In the case of in-plane modes, the caption (w=0) is added to underline that the vibration occurs in the plane and not out of plane. The third column shows all the results for the 3D exact solution, both classical and cylindrical bending (CB) results. It is noted that these frequencies are obtained with a single run of the 3D exact code keeping simply supported boundary conditions for all the four edges of the plate. The other four columns of Table 1 refer to the numerical solutions for simply supported structures and for cylindrical bending conditions (setting two free opposite edges). The cylindrical bending results are indicated by the acronym CB in the same table.Table 1
One-layered FGM plate withp=0.5. Analytical 3D versus 2D numerical models. Frequency f in Hz.
m , n Mode 3D exact 2D FE 2D FE (CB) 2D GDQ 2D GDQ (CB) a / h = 100 1,0 I 40.19 — 39.07 — 39.21 0,1 I 40.19 — 39.07 — 39.21 1,1 I 80.37 80.06 — 80.36 — 2,0 I 160.7 — 158.0 — 158.4 0,2 I 160.7 — 158.0 — 158.4 2,1 I 200.8 200.1 — 200.8 — 1,2 I 200.8 200.1 — 200.8 — 2,2 I 321.2 320.2 — 321.1 — 3,0 I 361.3 — 356.8 — 357.7 0,3 I 361.3 — 356.8 — 357.7 3,1 I 401.4 400.2 — 401.3 — 1,3 I 401.4 400.2 — 401.3 — 3,2 I 521.5 520.2 — 521.5 — 2,3 I 521.5 520.2 — 521.5 — 4,0 I 641.6 — 635.4 — 636.5 0,4 I 641.6 — 635.4 — 636.5 4,1 I 681.6 680.1 — 681.5 — 1,4 I 681.6 680.1 — 681.5 —  a / h = 20 1,0 I 200.2 — 195.1 — 195.1 0,1 I 200.2 — 195.1 — 195.1 1,1 I 398.9 399.5 — 398.8 — 2,0 I 791.6 — 786.4 — 779.7 0,2 I 791.6 — 786.4 — 779.7 2,1 I 985.8 995.4 — 985.7 — 1,2 I 985.8 995.4 — 985.7 — 2,2 I 1560 1588 — 1560 — 3,0 I 1749 — 1767 — 1729 0,3 I 1749 — 1767 — 1729 3,1 I 1936 1980 — 1936 — 1,3 I 1936 1980 — 1936 — 3,2 I 2490 2566 — 2490 — 2,3 I 2490 2566 — 2490 — 1,0 II (w=0) 2775 2773 — 2775 — 0,1 II (w=0) 2775 2773 — 2775 — 4,0 I 3034 — 3124 — 3007 0,4 I 3034 — 3124 — 3007There is a very good agreement among the three solutions when simply supported conditions are considered. On the contrary a small error is observed in cylindrical bending because the actual CB condition is not completely satisfied due to the fact that the plate has a limited dimension in the simply supported direction. This feature can be detected in Figures2 and 3 because the deformed shapes of the plate in the CB configuration are not properly correct. In fact, the fibers parallel to the supported edges are not parallel among them but there is a curvature effect due to the free boundaries and to the limited length of the plate. On the contrary, the in-plane mode shapes (w=0), when the plate is simply supported on the four edges, have all the fibers parallel among them also after the deformation. These modes are the results obtained in the 3D solution with m=0 or n=0 combined with w=0. They are not cylindrical bending modes; in fact they can be found with 2D numerical models without modifying the boundary conditions. 3D frequencies with m=0 or n=0 (withw different from 0) are CB modes which can be obtained considering two simply supported edges and two free edges in 2D numerical models.Figure 2
Differences between in-plane and cylindrical bending vibration modes for thea/h=20 one-layered FGM plate with p=0.5. 3D modes versus 2D modes.
(a)
Mode I: (0,1), 3D (b)
Mode I: (0,1), 2D (c)
Mode II: (w=0) (0,1), 3D (d)
Mode II: (w=0) (0,1), 2D (e)
Mode I: (1,0), 3D (f)
Mode I: (1,0), 2D  (g)
Mode II: (w=0) (1,0), 3D (h)
Mode II: (w=0) (1,0), 2DFigure 3
First four cylindrical bending vibration modes for thea/h=20 one-layered FGM plate with p=0.5. 3D modes versus 2D modes.
(a)
Mode I: (1,0), 3D (b)
Mode I: (1,0), 2D (c)
Mode I: (2,0), 3D (d)
Mode I: (2,0), 2D (e)
Mode I: (3,0), 3D (f)
Mode I: (3,0), 2D  (g)
Mode I: (4,0), 3D (h)
Mode I: (4,0), 2DAs it is well-known the “perfect” CB would occur when the simply supported direction has an “infinite” length. This cannot be proven numerically, but a trend can be observed modeling different plates in cylindrical bending configuration with different lengths. These analyses are made keeping one side of the plate constant and enlarging the length of the other edge. In the present computations,a is kept constant and b=L varies from 1 to 20. The results of the aforementioned analysis are listed in Table 2 and depicted in Figures 4 and 5. Table 2 includes the numerical values of the first mode shape in cylindrical bending for the 3D exact solution, the 2D FE, and the 2D GDQ models. Figure 4 shows that both GDQ and FE solutions tend to the 3D exact one increasing the length of the simply supported edge. In particular, the GDQ solution is slightly more accurate than the FE one. The 2D GDQ solution is carried out by enlarging the length of one of the sides of the plate and the same number of dofs is kept. On the contrary, the FE model is generated by 40 × 40 blocks of mesh. In other words, the plate with L=2a has a mesh of 40 × 80, the one with L=3a has a mesh 40 × 120, and so on. When L=20a, the error among the 3D exact and the 2D numerical solutions is very small and this is proven also by the mode shapes of Figure 5, wherein the boundary effect is negligible with respect the length of the simply supported edges. It must be underlined that it is impossible to completely avoid the boundary effects due to the free edges.Table 2
First frequencyf in Hz for the a/h=100 one-layered FGM plate with p=0.5: convergence of 2D GDQ and 2D FE solutions increasing the length of the simply supported sides.
L/a 1 2 3 4 5 6 8 10 15 20 3D exact 2D GDQ 39.21 39.64 39.80 39.88 39.93 39.97 40.01 40.04 40.07 40.08 40.19 2D FEM 39.07 39.49 39.66 39.74 39.79 39.83 39.87 39.89 39.93 39.94 40.19Figure 4
First frequencyf in Hz for the a/h=100 one-layered FGM plate with p=0.5: convergence of 2D GDQ and 2D FE solutions increasing the length of the simply supported sides.Figure 5
First vibration mode for thea/h=100 one-layered FGM plate with p=0.5: 2D GDQ solution and edge effects versus the length of the simply supported sides.
(a)
L / a = 2 (b)
L / a = 3 (c)
L / a = 4 (d)
L / a = 5 (e)
L / a = 6 (f)
L / a = 10The second example is related to a simply supported cylinder withRα=10m, Rβ=∞, a=2πRα, and b=20m and thickness ratios Rα/h=100 and Rα/h=10. A one-layered FGM scheme is considered with h1=h and the same volume fraction function of (24) for the previous example but setting p=2.0. The present material properties are the same as the flat plate. The FE model has a 20×80 Quad8 mesh with 27966 dofs and the GDQ solution has a 15×51 Chebyshev-Gauss-Lobatto grid. Numerically speaking, compatibility conditions must be enforced on the closing meridian in order to obtain a closed revolution shell [42]. To obtain the CB conditions numerically, the cylinder is modeled completely free and the first 6 rigid mode shapes are neglected and removed from the numerical solution. It is noted that the GDQ model uses a third-order kinematic expansion with no zigzag effect as the previous case, since a one-layer structure has been taken into account. Table 3 reports the results with the same meaning of the symbols seen in Table 1. A very good agreement is observed when the cylinder is simply supported, whereas, when the cylinder is completely free, some small errors can be observed but the numerical solutions follow the 3D exact one.Table 3
One-layered FGM cylinder withp=2.0. Analytical 3D versus 2D numerical models. Frequency f in Hz.
m , n Mode 3D exact 2D FE 2D FE (CB) 2D GDQ 2D GDQ (CB) R α / h = 100 4,0 I 0.896 — 0.897 — 0.883 6,0 I 2.532 — 2.538 — 2.511 8,0 I 4.854 — 4.867 — 4.819 10,0 I 7.848 — 7.871 — 7.801 12,0 I 11.51 — 11.55 — 11.46 10,1 I 13.56 13.60 — 13.57 — 12,1 I 14.41 14.49 — 14.45 — 14,0 I 15.84 — 15.89 — 15.77 8,1 I 16.38 16.40 — 16.36 — 14,1 I 17.56 17.67 — 17.63 — 16,0 I 20.83 — 20.90 — 20.70 16,1 I 22.07 22.22 — 22.16 — 6,1 I 24.41 24.42 — 24.38 — 18,0 I 26.49 — 26.58 — 26.40 14,2 I 27.42 27.52 — 27.46 — 18,1 I 27.51 27.72 — 27.63 — 16,2 I 28.75 28.90 — 28.83 — 12,2 I 29.21 29.26 — 29.22 —  Rα/h=10 4,0 I 8.774 — 9.085 — 8.679 6,0 I 24.56 — 25.71 — 24.30 6,1 I 39.78 40.89 — 39.79 — 4,1 I 43.88 43.57 — 43.88 — 8,0 I 46.54 — 49.20 — 46.02 8,1 I 55.75 59.05 — 55.76 — 2,1 I 72.02 71.01 — 72.02 — 2,0 I (w=0) 73.22 73.90 — 73.19 — 10,0 I 74.19 — 79.33 — 73.47 6,2 I 81.10 82.61 — 81.11 — 10,1 I 81.89 88.11 — 81.91 — 8,2 I 87.83 91.73 — 87.85 — 4,2 I 90.16 90.84 — 90.17 — 12,0 I 107.1 — 115.9 — 106.1 10,2 I 108.1 115.5 — 108.1 — 2,2 I 108.9 110.5 — 108.9 — 12,1 I 114.2 124.5 — 114.3 — 4,3 I 126.3 128.6 — 126.3 —Results withn=0 and w≠0 indicate the CB conditions. Such results are obtained by means of 2D numerical models modifying the boundary conditions (see 2D FE(CB) and 2D GDQ(CB) models). The cases n=0 combined with w=0 are not CB conditions and the frequencies can be found without modifying the boundary conditions in the 2D numerical models.The final test is related to a cylindrical panel withRα=10m, Rβ=∞, a=π/3Rα, and b=20m and thickness ratios Rα/h=100 and Rα/h=10. The present shell is a sandwich structure embedding a FGM core with h2=0.7h and two external skins with h1=h3=0.15h (h1+h2+h3=h). The bottom layer is made of isotropic metallic material and the top layer is made of ceramic material; the mechanical properties are(27)Em=73GPa,ρm=2800kg/m3,Ec=200GPa,ρc=5700kg/m3. The FGM core has a volume fraction function defined as(28)Vcz=0.5+z0.7hpwithz∈-0.35h,0.35h,where p=1.0 is employed in the present computations. The Poisson ratio is constant (νc=νm=ν=0.3). Equations (25) are still valid with Vc defined as in (28).The FE model is made of40×40 Quad8 elements for total 29280 dofs and the GDQ model has a 31×31 Chebyshev-Gauss-Lobatto grid. Since a sandwich structure is analyzed in the following, a third-order theory with zigzag function is considered in the computations for the GDQ method. The present natural frequencies are reported in Table 4 with the same meaning of the columns seen in the previous cases. In order to achieve the cylindrical bending conditions, the numerical models have the straight edges simply supported and the curved edges as free. Due to the CB definition, it is not possible to achieve CB conditions when m=0 and n=1,2,3,… because the cylindrical bending should occur on the curved edge and this is not physically possible because of the variation of displacement components due to the curvature presence. This aspect is indicated as NA (Not Applicable) solution in Table 4. Therefore, this kind of 3D exact solutions cannot be achieved by any numerical model. Once again, solution with n=0 or m=0 combined with w=0 (in-plane modes) can be found with 2D numerical models without modifying the simply supported boundary conditions.Table 4
Sandwich cylindrical shell panel embedding FGM core withp=1.0. Analytical 3D versus 2D numerical models. Frequency f in Hz.
m,n Mode 3D exact 2D FE 2D FE (CB) 2D GDQ 2D GDQ (CB) R α / h = 100 1,0 I 1.948 — 1.951 — 1.941 2,0 I 8.856 — 8.872 — 8.823 2,1 I 11.04 11.10 — 11.04 — 1,1 I 18.50 18.73 — 18.50 — 3,0 I 20.38 — 20.42 — 20.30 3,1 I 21.16 21.28 — 21.16 — 2,2 I 22.21 22.42 — 22.21 — 3,2 I 24.84 24.96 — 24.84 — 3,3 I 32.35 32.54 — 32.35 — 4,0 I 36.49 — 36.60 — 36.34 2,3 I 36.92 37.39 — 36.92 — 4,1 I 37.14 37.45 — 37.14 — 4,2 I 39.40 39.65 — 39.41 — 3,4 I 42.26 42.61 — 42.26 — 4,3 I 43.78 44.03 — 43.78 — 1,2 I 45.68 46.64 — 45.68 — 4,4 I 50.35 50.64 — 50.35 — 2,4 I 50.40 51.28 — 50.40 —  Rα/h=10 1,0 I 19.05 — 19.63 — 18.84 1,1 I 30.54 31.45 — 30.55 — 1,2 I 62.09 64.09 — 62.09 — 2,0 I 83.26 — 88.57 — 82.54 0,1 I 87.10 — NA — NA 0,1 II (w=0) 88.20 89.56 — 88.20 — 2,1 I 88.80 95.31 — 88.82 — 0,2 I 91.76 — NA — NA 1,3 I 95.76 99.83 — 95.67 — 0,3 I 103.6 — NA — NA 2,2 I 106.6 114.7 — 106.7 — 0,4 I 129.1 — NA — NA 1,4 I 133.5 141.3 — 133.5 — 2,3 I 136.4 147.5 — 136.4 — 1,0 II (w=0) 166.5 173.8 — 166.5 — 0,5 I 168.3 — NA — NA 2,4 I 175.5 191.8 — 175.5 — 0,2 II (w=0) 176.4 179.1 — 176.4 —


## 5. Conclusions

The present paper analyzes the cylindrical bending in the free frequency analysis of functionally graded material (FGM) plates and cylindrical shells. This condition is investigated by means of a 3D exact shell model and two different 2D numerical models, a weak formulation based on classical Finite Elements (FEs) and a strong formulation based on the GDQ method. The cylindrical bending (CB) behavior in the 3D exact shell model can be analyzed using completely simply supported structures and supposing one of the two half-wave numbers equals zero. The numerical models obtain the CB behavior using mixed simply supported and free boundary conditions. It has been demonstrated that the actual cylindrical bending (CB) is possible only with respect to straight edges. The CB natural frequencies of the numerical solutions converge to the 3D exact CB ones if the length of straight simply supported edges is increased up to 20 times the value of the original edges. The CB conditions are not possible in 2D numerical models when the half-wave number is zero in the curvilinear edge; in this case the modification of boundary conditions does not properly work. In-plane vibration modes with zero transverse displacement combined withm=0 or n=0 cannot be defined as CB frequencies and they can be calculated via 2D numerical models without the modification of the simply supported conditions for the four edges. All the considerations about the cylindrical bending are valid for each geometry (plate, cylinder, and cylindrical shell), lamination scheme, FGM law, and thickness ratio.



---

*Source: 2373862-2016-12-04.xml*