Analysis of Laminated Shells by Murakami ’ s Zig-Zag Theory and Radial Basis Functions Collocation

The static and free vibration analysis of laminated shells is performed by radial basis functions collocation, according toMurakami’s zig-zag (ZZ) function (MZZF) theory . The MZZF theory accounts for through-the-thickness deformation, by considering a ZZ evolution of the transverse displacement with the thickness coordinate. The equations of motion and the boundary conditions are obtained by Carrera’s Unified Formulation and further interpolated by collocation with radial basis functions.


Introduction
The efficient load-carrying capabilities of shell structures make them very useful in a variety of engineering applications.The continuous development of new structural materials leads to the ever increasingly complex structural designs that require careful analysis.Although analytical techniques are very important, the use of numerical methods to solve shell mathematical models of complex structures has become an essential ingredient in the design process.
The most common mathematical models used to describe shell structures may be classified into two classes according to different physical assumptions: the Koiter model [1], based on the Kirchhoff hypothesis, and the Naghdi model [2], based on the Reissner-Mindlin assumptions that take into account the transverse shear deformation.But these theories are not adequate to describe the so-called zig-zag (ZZ) effect in sandwich structures or layered composites, due to the discontinuity of mechanical properties between faces and core at the interfaces; see Figure 1 (to trace accurate responses of sandwich structures, see the books by Zenkert [3] and Vinson [4]).
The ZZ effect can be captured by the layerwise theories which typically assume independent degrees of freedom per layer.Unfortunately the computation can be prohibitive.The layerwise theories are reviewed in Burton and Noor [5], Noor et al. [6], Altenbach [7], Librescu and Hause [8], Vinson [9], and Demasi [10].In order to overcome the computational cost of the layerwise theories, Murakami [11] proposed a zig-zag function (ZZF) that is able to reproduce the slope discontinuity.Equivalent Single Layer models with only displacement unknowns can be developed on the basis of ZZF.A review of the application of ZZF in plates and shells was presented by Carrera [12][13][14][15][16] and some relevant papers on the analysis of sandwich structures were presented in [17][18][19].
The most common numerical procedure for the analysis of the shells is the finite element method [20][21][22][23][24].It is known that the phenomenon of numerical locking may arise from hidden constrains that are not well represented in the finite element approximation and, in the scientific literature, it is possible to find many methods to overcome this problem [25][26][27][28][29][30].The present paper, that performs the bending and free vibration analysis of laminated shells by collocation with radial basis functions, avoids the locking phenomenon.A radial basis function, (‖ −   ‖), is a spline that depends on the Euclidian distance between distinct data centers  , = 1, 2, . . .,  ∈ R  , also called nodal or collocation points.We use the so-called unsymmetrical Kansa method that was introduced by Kansa [31].The use of radial basis function for the analysis of structures and materials has been previously studied by numerous authors [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46].The authors have recently applied the RBF collocation to the static deformations of composite beams and plates [47][48][49].One of the authors has already combined Reddy's theory with radial basis functions in [50].In fact, Reddy's theory is quite efficient for laminated (monolithic) composite plates or shells but not as efficient (or adequate) for sandwich structures because of the very high difference on material properties from the skins and the core.Reddy's theory does not allow the thickness stretching, but our formulation is more general and allows any expansion in the thickness direction.This is where the present paper shows more interest.
In this paper for the first time how the Unified Formulation can be combined with radial basis functions to the analysis of thin and thick laminated shells, using Murakami's zig-zag function, allowing for through-the-thickness deformations, is investigated.The quality of the present method in predicting static deformations and free vibrations of thin and thick laminated shells is compared and discussed with other methods in some numerical examples.

Applying the Unified Formulation to MZZF
The Unified Formulation (UF) proposed by Carrera [13,[51][52][53], also known as CUF, is a powerful framework for the analysis of beams, plates, and shells.This formulation has been applied in several finite element analyses, either by using the Principle of Virtual Displacement, or by using Reissner's mixed variational theorem.The stiffness matrix components, the external force terms, or the inertia terms can be obtained directly with this UF, irrespective of the shear deformation theory being considered.
In this section Carrera's Unified Formulation is briefly reviewed.How to obtain the fundamental nuclei, which allows the derivation of the equations of motion and boundary conditions, in weak form for the finite element analysis and in strong form for the present RBF collocation, is shown.

The MZZF Theory.
Let us consider a sandwich plate (translation to shells becomes evident later in the paper) composed of three perfectly bonded layers, with  being the thickness coordinate of the whole plate while   is the layer thickness coordinate; see Figure 1. and ℎ are length and thickness of the square laminated plate, respectively.The adimensional layer coordinate   = (2  )/ℎ  is further introduced (ℎ  is the thickness of the kth layer).Murakami's zig-zag function () was defined according to the following formula [11]: () has the following properties.
(1) It is a piecewise linear function of layer coordinates   .
(2) () has unit amplitude for the whole layers.
(3) The slope   () = / assumes opposite sign between two adjacent layers.Its amplitude is layer thickness independent.
A possible FSDT theory has been investigated by Carrera [14] and Demasi [15], ignoring the through-the-thickness deformations: A refinement of FSDT by inclusion of ZZ effects and transverse normal strains was introduced in Murakami's original ZZF, defined by the following displacement field: where   and  +1 are the bottom and top -coordinates at each layer.The additional degrees of freedom   and V  have a meaning of displacement, and their amplitude is layer independent.

Governing Equations and Boundary Conditions in the
Framework of Unified Formulation.Shells are bidimensional structures in which one dimension (in general the thickness in  direction) is negligible with respect to the other two inplane dimensions.Geometry and the reference system are indicated in Figure 3.The square of an infinitesimal linear segment in the layer and the associated infinitesimal area and volume are given by where the metric coefficients are denotes the -layer of the multilayered shell;    and    are the principal radii of curvature along the coordinates  and , respectively.  and   are the coefficients of the first fundamental form of Ω  (Γ  is the Ω  boundary).In this work, the attention has been restricted to shells with constant radii of curvature (cylindrical, spherical, and toroidal geometries) for which   =   = 1.
Although one can use the UF for one-layer, isotropic shell, a multilayered shell with   layers is considered.The Principle of Virtual Displacement (PVD) for the pure-mechanical case reads where Ω  and   are the integration domains in plane (, ) and  direction, respectively.Here,  indicates the layer and  the transpose of a vector, and    is the external work for the th layer. means geometrical relations and  constitutive equations.
The steps to obtain the governing equations are (i) substitution of the geometrical relations (subscript ), (ii) substitution of the appropriate constitutive equations (subscript ), (iii) introduction of the Unified Formulation.
Stresses and strains are separated into in-plane and normal components, denoted, respectively, by the subscripts  and .The mechanical strains in the th layer can be related to the displacement field u  = {   ,    ,    } via the geometrical relations: The explicit form of the introduced arrays is as follows The 3D constitutive equations are given as with According to the Unified Formulation by Carrera, the three displacement components   ,   , and   and their relative variations can be modelled as with Taylor expansions from the first up to the 4th order:  0 =  0 = 1,  1 =  1 = , . ..,   =   , . .., and  4 =  4 if an Equivalent Single Layer (ESL) approach is used.
Resorting to the displacement field in (3), we choose vectors for displacement , V, and .We then obtain all terms of the equations of motion by integrating through-the-thickness direction.
It is interesting to note that, under this combination of the Unified Formulation and RBF collocation, the collocation code depends only on the choice of   , in order to solve this type of problems.We designed a MATLAB code that just by changing   can analyse static deformations, free vibrations, and buckling loads for any type of  ∘ shear deformation theory.An obvious advantage of the present methodology is that the tedious derivation of the equations of motion and boundary conditions for a particular shear deformation theory is no longer an issue, as this MATLAB code does all that work for us.
In Figure 2 the assembling procedures are shown on layer  for ESL approach.
Substituting the geometrical relations, the constitutive equations, and the Unified Formulation into the variational statement PVD, for the th layer, one has At this point, the formula of integration by parts is applied where I Ω matrix is obtained applying the gradient theorem with   being the components of the normal n to the boundary along the direction .After integration by parts and the substitution of CUF, the governing equations and boundary conditions for the shell in the mechanical case are obtained: where I   and I   depend on the boundary geometry: The normal to the boundary of domain where   and   are the angles between the normal n and the directions  and , respectively.The governing equations for a multilayered shell subjected to mechanical loadings are Figure 3: Geometry and notations for a multilayered shell (doubly curved).
where the fundamental nucleus K   is obtained as And the corresponding Neumann-type boundary conditions on Γ  are where And P   are variationally consistent loads with applied pressure.

Fundamental
The fundamental nucleus K   is reported for doubly curved shells (radii of curvature in both  and  directions; see Figure 3): One can note that all the equations written for the shell degenerate into those for the plate when 1/  = 1/  = 0.
In practice we set the radii of curvature to 10 9 .

Dynamic Governing Equations.
The PVD for the dynamic case is expressed as where   is the mass density of the th layer and double dots denote acceleration.By substituting the geometrical relations, the constitutive equations, and the Unified Formulation, we obtain the following governing equations: In the case of free vibrations one has where M  is the fundamental nucleus for the inertial term.
The explicit form of that is where the meaning of the integral    has been illustrated in (25).The geometrical and mechanical boundary conditions are the same of the static case.Because we consider the static case only, the mass terms will be neglected.

The Radial Basis Function Method
3.1.The Static Problem.Radial basis functions (RBFs) approximations are mesh-free numerical schemes that can exploit accurate representations of the boundary, are easy to implement, and can be spectrally accurate.In this section the formulation of a global unsymmetrical collocation RBFbased method to compute elliptic operators is presented.
Consider a linear elliptic partial differential operator  and a bounded region Ω in R  with some boundary Ω.In the static problems we seek the computation of displacement (u) from the global system of equations where L and L  are linear operators in the domain and on the boundary, respectively.The right-hand sides of (32) represent the external forces applied on the plate or shell and the boundary conditions applied along the perimeter of the plate or shell, respectively.The PDE problem defined in (32) will be replaced by a finite problem, defined by an algebraic system of equations, after the radial basis expansions.

The Eigenproblem.
The eigenproblem looks for eigenvalues () and eigenvectors (u) that satisfy As in the static problem, the eigenproblem defined in ( 33) is replaced by a finite-dimensional eigenvalue problem, based on RBF approximations.

Radial Basis Functions Approximations. The radial basis function (𝜙) approximation of a function (u) is given by
where   ,  = 1, . . ., , is a finite set of distinct points (centers) in R  .The most common RBFs are Cubic:  () =  3 , Thin plate splines:  () =  2 log () , Wendland functions:  () = (1 − )  +  () , where the Euclidian distance  is real and nonnegative and  is a positive shape parameter.Hardy [54] introduced multiquadrics in the analysis of scattered geographical data.In the 1990s Kansa [31] used multiquadrics for the solution of partial differential equations.Considering  distinct interpolations and knowing (  ),  = 1, 2, . . ., , we find   by the solution of a  ×  linear system where At the points in the domain we solve the following system of equations: or where At the points on the boundary, we impose boundary conditions as or where Therefore, we can write a finite-dimensional static problem as By inverting the system (43), we obtain the vector .We then obtain the solution u using the interpolation (34).

Solution of the Eigenproblem.
We consider   nodes in the interior of the domain and   nodes on the boundary, with  =   +   .We denote interpolation points by   ∈ Ω,  = 1, . . .,   , and   ∈ Ω,  =   + 1, . . ., .At the points in the domain, we define the eigenproblem as or where At the points on the boundary, we enforce the boundary conditions as or Equations ( 45) and ( 48) can now be solved as a generalized eigenvalue problem where

Discretization of the Equations of Motion and Boundary
Conditions.The radial basis collocation method follows a simple implementation procedure.Taking (11), we compute This  vector is then used to obtain solution ũ, by using (5).
If derivatives of ũ are needed, such derivatives are computed as and so forth.
In the present collocation approach, we need to impose essential and natural boundary conditions.Consider, for example, the condition  = 0, on a simply supported or clamped edge.We enforce the conditions by interpolating as Other boundary conditions are interpolated in a similar way.

Free Vibrations Problems.
For free vibration problems we set the external force to zero and assume harmonic solution in terms of displacement,  0 ,  1 , V 0 , V 1 , . .., as where  is the frequency of natural vibration.Substituting the harmonic expansion into (49) in terms of the amplitudes  0 ,  1 ,   ,  0 ,  1 ,   ,  0 ,  1 , and  2 , we may obtain the natural frequencies and vibration modes for the plate or shell problem, by solving the eigenproblem where L collects all stiffness terms and G collects all terms related to the inertial terms.In (55) X are the modes of vibration associated with the natural frequencies defined as .

Numerical Examples
All numerical examples consider a Chebyshev grid and a Wendland function, defined as where the shape parameter () was obtained by an optimization procedure, as detailed in Ferreira and Fasshauer [55].

Spherical Shell in
with the origin of the coordinate system located at the lower left corner on the midplane and  the maximum load (at center of shell).The orthotropic material properties for each layer are given by The in-plane displacement, the transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as The shell is simply supported on all edges.
In Table 1 we compare the static deflections for the present shell model with results of Reddy shell formulation using first-order and third-order shear deformation theories [56].We consider nodal grids with 13 × 13, 17 × 17, and 21 × 21 points.We consider various values of / and two values of /ℎ (10 and 100).Results are in good agreement for various /ℎ ratios with the higher-order results of Reddy and Liu [56].

Free Vibration of Spherical and Cylindrical Laminated
Shells.We consider nodal grids with 13 × 13, 17 × 17, and 21 × 21 points.In Tables 2 and 3 we compare the nondimensionalized natural frequencies from the present Murakami theory for various cross-ply spherical shells, with analytical solutions done by Reddy and Liu [56] who considered both the first-order (FSDT) and the third-order (HSDT) theories.The first-order theory overpredicts the fundamental natural frequencies of symmetric thick shells and symmetric shallow thin shells.The present radial basis function method is  compared with analytical results by Reddy and Liu [56] and shows excellent agreement.
Table 4 contains nondimensionalized natural frequencies obtained using the the present Murakami theory for crossply cylindrical shells with lamination schemes [0/90/0], [0/90/90/0].Present results are compared with analytical solutions done by Reddy and Liu [56] who considered both the first-order (FSDT) and the third-order (HSDT) theories.The present radial basis function method is compared with analytical results by Reddy and Liu [56] and shows excellent agreement.

Concluding Remarks
In this paper Murakami's theory was implemented for the first time for laminated orthotropic elastic shells through a multiquadric discretization of equations of motion and boundary conditions.The multiquadric radial basis function method for the solution of shell bending and free vibration problems was presented.Results for static deformations and natural frequencies were obtained and compared with other sources.This meshless approach demonstrated that is very successful in the static deformations and free vibration analysis of laminated composite shells.Advantages of radial basis functions are absence of mesh, ease of discretization of boundary conditions and equations of equilibrium or motion, and very easy coding.We show that the static displacement and stresses and the natural frequencies obtained from present method are in excellent agreement with analytical solutions.

Figure 1 :
Figure 1: Scheme of the zig-zag assumptions for a three-layered laminate.