Journal of Applied Mathematics

Journal of Applied Mathematics / 2013 / Article

Research Article | Open Access

Volume 2013 |Article ID 123465 |

D. A. Maturi, A. J. M. Ferreira, A. M. Zenkour, D. S. Mashat, "Analysis of Laminated Shells by Murakami’s Zig-Zag Theory and Radial Basis Functions Collocation", Journal of Applied Mathematics, vol. 2013, Article ID 123465, 14 pages, 2013.

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

Academic Editor: Song Cen
Received10 Jul 2013
Accepted15 Sep 2013
Published11 Nov 2013


The static and free vibration analysis of laminated shells is performed by radial basis functions collocation, according to Murakami’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.

1. 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 [1216] and some relevant papers on the analysis of sandwich structures were presented in [1719].

The most common numerical procedure for the analysis of the shells is the finite element method [2024]. 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 [2530]. 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 , 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 [3246]. The authors have recently applied the RBF collocation to the static deformations of composite beams and plates [4749]. 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.

2. Applying the Unified Formulation to MZZF

The Unified Formulation (UF) proposed by Carrera [13, 5153], 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.

2.1. 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 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 are the bottom and top -coordinates at each layer. The additional degrees of freedom and have a meaning of displacement, and their amplitude is layer independent.

2.2. 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 in-plane 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 .

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 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:  , , , and if an Equivalent Single Layer (ESL) approach is used.

Resorting to the displacement field in (3), we choose vectors for displacement , , 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 matrix is obtained applying the gradient theorem with being the components of the normal 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 and depend on the boundary geometry:

The normal to the boundary of domain is where and are the angles between the normal and the directions and , respectively.

The governing equations for a multilayered shell subjected to mechanical loadings are where the fundamental nucleus is obtained as And the corresponding Neumann-type boundary conditions on are where And are variationally consistent loads with applied pressure.

2.3. Fundamental Nuclei

The following integrals are introduced to perform the explicit form of fundamental nuclei:

The fundamental nucleus is reported for doubly curved shells (radii of curvature in both and directions; see Figure 3):

The application of boundary conditions makes use of the fundamental nucleus in the form:

One can note that all the equations written for the shell degenerate into those for the plate when . In practice we set the radii of curvature to .

2.4. 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 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.

3. 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 RBF-based method to compute elliptic operators is presented.

Consider a linear elliptic partial differential operator and a bounded region in with some boundary . In the static problems we seek the computation of displacement from the global system of equations where and 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.

3.2. The Eigenproblem

The eigenproblem looks for eigenvalues () and eigenvectors that satisfy

As in the static problem, the eigenproblem defined in (33) is replaced by a finite-dimensional eigenvalue problem, based on RBF approximations.

3.3. Radial Basis Functions Approximations

The radial basis function approximation of a function is given by where , is a finite set of distinct points (centers) in . The most common RBFs are 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 , , we find by the solution of a linear system where , , and .

3.4. Solution of the Static Problem

The solution of a static problem by radial basis functions considers nodes in the domain and nodes on the boundary, with a total number of nodes . We denote the sampling points by , , and , . 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 using the interpolation (34).

3.5. Solution of the Eigenproblem

We consider nodes in the interior of the domain and nodes on the boundary, with . We denote interpolation points by , , and , . 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

3.6. 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 , on a simply supported or clamped edge. We enforce the conditions by interpolating as

Other boundary conditions are interpolated in a similar way.

3.7. Free Vibrations Problems

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

4. 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].

4.1. Spherical Shell in Bending

A laminated composite spherical shell is here considered, of side and thickness , to be composed of layers oriented at and . The shell is subjected to a sinusoidal vertical pressure of the form 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 , , and 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].

5 10 20 50 100

10 Present ( ) 6.8510 7.0818 7.1429 7.1609 7.1637 7.1649
10 Present ( ) 6.8516 7.0822 7.1433 7.1612 7.1640 7.1653
10 Present ( ) 6.8516 7.0822 7.1433 7.1612 7.1640 7.1653
10 HSDT [56] 6.7688 7.0325 7.1016 7.1212 7.1240 7.125
10 FSDT [56] 6.4253 6.6247 6.6756 6.6902 6.6923 6.6939
100 Present ( ) 1.0245 2.3658 3.5167 4.0714 4.1652 4.1975
100 Present ( ) 1.0250 2.3667 3.5181 4.0728 4.1667 4.1990
100 Present ( ) 1.0250 2.3669 3.5183 4.0731 4.1669 4.1992
100 HSDT [56] 1.0321 2.4099 3.617 4.2071 4.3074 4.3420
100 FSDT [56] 1.0337 2.4109 3.6150 4.2027 4.3026 4.3370

10 Present ( ) 6.7737 7.0012 7.0614 7.0791 7.0819 7.0831
10 Present ( ) 6.7742 7.0015 7.0618 7.0795 7.0822 7.0834
10 Present ( ) 6.7742 7.0015 7.0618 7.0795 7.0822 7.0834
10 HSDT [56] 6.7865 7.0536 7.1237 7.1436 7.1464 7.1474
10 FSDT [56] 6.3623 6.5595 6.6099 6.6244 6.6264 6.6280
100 Present ( ) 1.0190 2.3583 3.5125 4.0702 4.1647 4.1972
100 Present ( ) 1.0194 2.3593 3.5138 4.0717 4.1662 4.1987
100 Present ( ) 1.0195 2.3594 3.5140 4.0719 4.1664 4.1989
100 HSDT [56] 1.0264 2.4024 3.6133 4.2071 4.3082 4.3430
100 FSDT [56] 1.0279 2.4030 3.6104 4.2015 4.3021 4.3368

4.2. Free Vibration of Spherical and Cylindrical Laminated Shells

We consider nodal grids with , , and 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.

5 10 20 50 100

10 Present ( ) 12.0527 11.8889 11.8474 11.8357 11.8340 11.8335
Present ( ) 12.0523 11.8886 11.8471 11.8355 11.8338 11.8332
Present ( ) 12.0522 11.8885 11.8470 11.8355 11.8338 11.8332
HSDT [56] 12.040 11.840 11.790 11.780 11.780 11.780

100 Present ( ) 31.2170 20.5742 16.8698 15.6744 15.4960 15.4361
Present ( ) 31.2072 20.5679 16.8648 15.6698 15.4915 15.4316
Present ( ) 31.2059 20.5672 16.8642 15.6693 15.4910 15.4311
HSDT [56] 31.100 20.380 16.630 15.420 15.230 15.170

5 10 20 50 100

10 Present ( ) 11.9831 11.8192 11.7777 11.7660 11.7643 11.7638
Present ( ) 11.9827 11.8190 11.7774 11.7658 11.7641 11.7635
Present ( ) 11.9827 11.8190 11.7774 11.7658 11.7641 11.7635
HSDT [56] 12.060 11.860 11.810 11.790 11.790 11.790

100 Present ( ) 31.1343 20.5420 16.8595 15.6721 15.4950 15.4355
Present ( ) 31.1244 20.5357 16.8545 15.6675 15.4905 15.4310
Present ( ) 31.1231 20.5350 16.8540 15.6671 15.4901 15.4306
HSDT [56] 31.020 20.350 16.620 15.420 15.240 15.170

Table 4 contains nondimensionalized natural frequencies obtained using the the present Murakami theory for cross-ply 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.

Method [0/90/0] [0/90/90/0]

5 Present ( ) 20.4956 11.7774 20.5298 11.8524
Present ( ) 20.4839 11.7770 20.5201 11.8521
Present ( ) 20.4829 11.7770 20.5189 11.8521
FSDT [56] 20.332 12.207 20.361 12.267
HSDT [56] 20.330 11.850 20.360 11.830

10 Present ( ) 16.8475 11.7672 16.8583 11.8382
Present ( ) 16.8409 11.7669 16.8522 11.8380
Present ( ) 16.8401 11.7669 16.8516 11.8380
FSDT [56] 16.625 12.173 16.634 12.236
HSDT [56] 16.620 11.800 16.630 11.790

20 Present ( ) 15.8006 11.7646 15.8039 11.8346
Present ( ) 15.7955 11.7644 15.7990 11.8344
Present ( ) 15.7950 11.7644 15.7985 11.8344
FSDT [56] 15.556 12.166 15.559 12.230
HSDT [56] 15.55 11.79 15.55 11.78

50 Present ( ) 15.4945 11.7639 15.4955 11.8336
Present ( ) 15.4899 11.7637 15.4909 11.8334
Present ( ) 15.4895 11.7637 15.4905 11.8334
FSDT [56] 15.244 12.163 15.245 12.228
HSDT [56] 15.24 11.79 15.23 11.78

100 Present ( ) 15.4503 11.7638 15.4510 11.8335
Present ( ) 15.4457 11.7636 15.4464 11.8333
Present ( ) 15.4453 11.7636 15.4460 11.8333
FSDT [56] 15.198 12.163 15.199 12.227
HSDT [56] 15.19 11.79 15.19 11.78

Plate Present ( ) 15.4355 11.7638 15.4361 11.8335
Present ( ) 15.4310 11.7635 15.4316 11.8332
Present ( ) 15.4306 11.7635 15.4311 11.8332
FSDT [56] 15.183 12.162 15.184 12.226
HSDT [56] 15.170 11.790 15.170 11.780

In Figure 4 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, , for a laminate (), using a grid of points, for , . The modes of vibration are quite stable.

In Figure 5 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, , for a laminate (), using a grid of points, for , . Again the modes of vibration are quite stable.

5. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This paper was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under Grant no. HiCi/1433/363-4. The authors, therefore, acknowledge with thanks DSR technical and financial support.


  1. W. T. Koiter, “On the foundations of the linear theory of thin elastic shells,” Koninklijke Nederlandse Akademie van Wetenschappen, vol. 73, pp. 169–195, 1970. View at: Google Scholar | Zentralblatt MATH | MathSciNet
  2. P. M. Naghdi, “A survey of recent progress in the theory of elastic shells,” Applied Mechanics Reviews, vol. 9, pp. 365–368, 1956. View at: Google Scholar
  3. D. Zenkert, An Introduction to Sandwich Structures, Chamelon Press, Oxford,