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

Authors: D. A. Maturi; A. J. M. Ferreira; A. M. Zenkour; D. S. Mashat
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

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.


---

## Body

## 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]).Figure 1
Scheme of the zig-zag assumptions for a three-layered laminate.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–16] and some relevant papers on the analysis of sandwich structures were presented in [17–19].The most common numerical procedure for the analysis of the shells is the finite element method [20–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–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, ϕ(∥x-xj∥), is a spline that depends on the Euclidian distance between distinct data centers xj,j=1,2,…,N∈ℝn, 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–46]. The authors have recently applied the RBF collocation to the static deformations of composite beams and plates [47–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.


## 2. Applying the Unified Formulation to MZZF

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


### 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, withz being the thickness coordinate of the whole plate while zk is the layer thickness coordinate; see Figure 1. a and h are length and thickness of the square laminated plate, respectively. The adimensional layer coordinate ζk=(2zk)/hk is further introduced (hk is the thickness of the kth layer). Murakami's zig-zag function Z(z) was defined according to the following formula [11]:
(1)Z(z)=(-1)kζz.Z
(
z
) has the following properties.  (1)
It is a piecewise linear function of layer coordinateszk.(2)
Z(z) has unit amplitude for the whole layers.(3)
The slopeZ′(z)=dZ/dz 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:
(2)u=u0+zu1+(-1)k2hk(z-12(zk+zk+1))uZ,v=v0+zv1+(-1)k2hk(z-12(zk+zk+1))vZ,w=w0.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:(3)u=u0+zu1+(-1)k2hk(z-12(zk+zk+1))uZ,(4)v=v0+zv1+(-1)k2hk(z-12(zk+zk+1))vZ,(5)w=w0+zw1+(-1)k2hk(z-12(zk+zk+1))wZ,wherezk and zk+1 are the bottom and top z-coordinates at each layer. The additional degrees of freedom uZ and vZ 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 inz 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
(6)dsk2=Hαk2dα2+Hβk2dβ2+Hzk2dz2,dΩk=HαkHβkdαdβ,dV=HαkHβkHzkdαdβdz,where the metric coefficients are(7)Hαk=Ak(1+zRαk),Hβk=Bk(1+zRβk),Hzk=1.kdenotes the k-layer of the multilayered shell; Rαk and Rβk are the principal radii of curvature along the coordinates α and β, respectively. Ak and Bk are the coefficients of the first fundamental form of Ωk (Γk is the Ωk boundary). In this work, the attention has been restricted to shells with constant radii of curvature (cylindrical, spherical, and toroidal geometries) for which Ak=Bk=1.Although one can use the UF for one-layer, isotropic shell, a multilayered shell withNl layers is considered. The Principle of Virtual Displacement (PVD) for the pure-mechanical case reads
(8)∑k=1Nl‍∫Ωk‍∫Ak‍{δϵpGkTσpCk+δϵnGkTσnCk}×dΩkdz=∑k=1Nl‍δLek,
where Ωk and Ak are the integration domains in plane (α,β) and z direction, respectively. Here, k indicates the layer and T the transpose of a vector, and δLek is the external work for the kth layer. G means geometrical relations and C constitutive equations.The steps to obtain the governing equations are(i)
substitution of the geometrical relations (subscriptG),(ii)
substitution of the appropriate constitutive equations (subscriptC),(iii)
introduction of the Unified Formulation.Stresses and strains are separated into in-plane and normal components, denoted, respectively, by the subscriptsp and n. The mechanical strains in the kth layer can be related to the displacement field uk={uαk,uβk,uzk} via the geometrical relations:
(9)ϵpGk=[ϵααk,ϵββk,ϵαβk]T=(Dpk+Apk)uk,ϵnGk=[ϵαzk,ϵβzk,ϵzzk]T=(DnΩk+Dnzk-Ank)uk.The explicit form of the introduced arrays is as follows(10)Dpk=[∂αHαk000∂βHβk0∂βHβk∂αHαk0],DnΩk=[00∂αHαk00∂βHβk000],Dnzk=[∂z000∂z000∂z],(11)Apk=[001HαkRαk001HβkRβk000],Ank=[1HαkRαk0001HβkRβk0000].The 3D constitutive equations are given as(12)σpCk=CppkϵpGk+CpnkϵnGk,σnCk=CnpkϵpGk+CnnkϵnGk
with
(13)Cppk=[C11kC12kC16kC12kC22kC26kC16kC26kC66k],Cpnk=[00C13k00C23k00C36k],Cnpk=[000000C13kC23kC36k],Cnnk=[C55kC45k0C45kC44k000C33k].According to the Unified Formulation by Carrera, the three displacement componentsuα, uβ, and uz and their relative variations can be modelled as
(14)(uα,uβ,uz)=Fτ(uατ,uβτ,uzτ),(δuα,δuβ,δuz)=Fs(δuαs,δuβs,δuzs)
with Taylor expansions from the first up to the 4th order:  F0=z0=1, F1=z1=z,…, FN=zN,…, and F4=z4 if an Equivalent Single Layer (ESL) approach is used.Resorting to the displacement field in (3), we choose vectors Ft=[1z(-1)k(2/hk)(z-(1/2)(zk+zk+1))] for displacement u, v, and w. 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 ofFt, in order to solve this type of problems. We designed a MATLAB code that just by changing Ft can analyse static deformations, free vibrations, and buckling loads for any type of C∘ 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 Figure2 the assembling procedures are shown on layer k for ESL approach.Figure 2
Assembling procedure for ESL approach.Geometry and notations for a multilayered shell (doubly curved).
(a)
(b)Substituting the geometrical relations, the constitutive equations, and the Unified Formulation into the variational statement PVD, for thekth layer, one has
(15)∑k=1Nl‍{∫Ωk‍∫Ak‍{((Dp+Ap)δuk)T11111111111×(Cppk(Dp+Ap)uk111111111111111+Cpnk(DnΩ+Dnz-An)uk)11111111111+((DnΩ+Dnz-An)δuk)T11111111111×(Cnpk(Dp+Ap)uk111111111111111+Cnnk(DnΩ+Dnz-An)uk){((Dp+Ap)δuk)T}dΩkdzk{∫Ωk‍∫Ak‍{((Dp+Ap)δuk)T}11111111111=∑k=1Nl‍δLek.At this point, the formula of integration by parts is applied(16)∫Ωk‍((DΩ)δak)TakdΩk=-∫Ωk‍δakT((DΩT)ak)dΩk∫Ωk‍((DΩ)δak)TakdΩk=+∫Γk‍δakT((IΩ)ak)dΓk,
where IΩ matrix is obtained applying the gradient theorem(17)∫Ω‍∂ψ∂xidυ=∮Γ‍niψds
with ni being the components of the normal n^ to the boundary along the direction i. After integration by parts and the substitution of CUF, the governing equations and boundary conditions for the shell in the mechanical case are obtained:
(18)∑k=1Nl‍{∫Ωk‍∫Ak‍{δuskT[(-Dp+Ap)T11111111111111111×Fs(Cppk(Dp+Ap)Fτuτk1111111111111111111111+Cpnk(DnΩ+Dnz-An)Fτuτk)[(-Dp+Ap)T]111111111111+δuskT[+Cnnk(DnΩ+Dnz-An)Fτuτk)](-DnΩ+Dnz-An)T11111111111111111111×Fs(Cnpk(Dp+Ap)Fτuτk1111111111111111111111111+Cnnk(DnΩ+Dnz1111111111111111111111111-An)Fτuτk)]{δuskT[(-Dp+Ap)T}dΩkdzk{∫Ωk‍∫Ak‍{δuskT[(-Dp+Ap)T}+∑k=1Nl‍{∫Γk‍∫Ak‍{×Fτuτk)]}δuskT[IpTFs(Cppk(Dp+Ap)Fτuτk1111111111111111111111111+Cpnk(DnΩ+Dnz11111111111111111111111111111111-An)Fτuτk)]111111111111111+δuskT[InpTFs(Cnpk(Dp-Ap)Fτuτk111111111111111111111111111+Cnnk(DnΩ+Dnz-An)111111111111111111111111111×Fτuτk)]}dΓkdzk{∫Γk‍∫Ak‍{δuskT[IpTFs(Cppk(Dp+Ap)Fτuτk}=∑k=1Nl‍{∫ΩkδuskTFspuk},
where Ipk and Inpk depend on the boundary geometry:
(19)Ip=[nαHα000nβHβ0nβHβnαHα0],Inp=[00nαHα00nβHβ000].The normal to the boundary of domainΩ is
(20)n^=[nαnβ]=[cos(φα)cos(φβ)],
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(21)δuskT:Kuukτsuτk=Puτk,
where the fundamental nucleus Kuukτs is obtained as
(22)Kuukτs=∫Ak‍[[-Dp+Ap]TCppk[Dp+Ap]1111111111+[-Dp+Ap]TCpnk[DnΩ+Dnz-An]1111111111+[-DnΩ+Dnz-An]TCnpk[Dp+Ap]1111111111+[-DnΩ+Dnz-An]T1111111111×Cnnk[DnΩ+Dnz-An][[-Dp+Ap]TCppk[Dp+Ap]]FτFsHαkHβkdz.
And the corresponding Neumann-type boundary conditions on Γk are
(23)Πdkτsuτk=Πdkτsu-τk,
where
(24)Πdkτs=∫Ak‍[×[DnΩ+Dnz-Anτ]]IpTCppk[Dp+Apτ]+IpTCpnk[DnΩ+Dnz-Anτ]1111111111+InpTCnpk[Dp+Apτ]+InpTCnnk1111111111×[DnΩ+Dnz-Anτ][×[DnΩ+Dnz-Anτ]]IpTCppk[Dp+Apτ]+IpTCpnk[DnΩ+Dnz-Anτ]]FτFsHαkHβkdz.
And Puτk are variationally consistent loads with applied pressure.


### 2.3. Fundamental Nuclei

The following integrals are introduced to perform the explicit form of fundamental nuclei:(25)(Jkτs,Jαkτs,Jβkτs,Jα/βkτs,Jβ/αkτs,Jαβkτs)=∫Ak‍FτFs(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτzs,Jαkτzs,Jβkτzs,Jα/βkτzs,Jβ/αkτzs,Jαβkτzs)=∫Ak‍∂Fτ∂zFs(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτsz,Jαkτsz,Jβkτsz,Jα/βkτsz,Jβ/αkτsz,Jαβkτsz)=∫Ak‍Fτ∂Fs∂z(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτzsz,Jαkτzsz,Jβkτzsz,Jα/βkτzsz,Jβ/αkτzsz,Jαβkτzsz)=∫Ak‍∂Fτ∂z∂Fs∂z(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz.The fundamental nucleusKuukτs is reported for doubly curved shells (radii of curvature in both α and β directions; see Figure 3):
(26)(Kuuτsk)11=-C11kJβ/αkτs∂αs∂ατ-C16kJkτs∂ατ∂βs(Kuuτsk)11=-C16kJkτs∂αs∂βτ-C66kJα/βkτs∂βs∂βτ(Kuuτsk)11=+C55k(Jαβkτzsz-1RαkJβkτzs1111111111111111-1RαkJβkτsz+1Rαk2Jβ/αkτs),(Kuuτsk)12=-C12kJkτs∂ατ∂βs-C16kJβ/αkτs∂αs∂ατ(Kuuτsk)12=-C26kJα/βkτs∂βs∂βτ-C66kJkτs∂αs∂βτ(Kuuτsk)12=+C45k(Jαβkτzsz-1RβkJαkτzs1111111111111111-1RαkJβkτsz+1Rαk1RβkJkτs),(Kuuτsk)13=-C11k1RαkJβ/αkτs∂ατ-C12k1RβkJkτs∂ατ(Kuuτsk)13=-C13kJβkτsz∂ατ-C16k1RαkJkτs∂βτ(Kuuτsk)13=-C26k1RβkJα/βkτs∂βτ-C36kJαkτsz∂βτ(Kuuτsk)13=+C45k(Jαkτzs∂βs-1RαkJkτs∂βs)(Kuuτsk)13=+C55k(Jβkτzs∂αs-1RαkJβ/αkτs∂αs),(Kuuτsk)21=-C12kJkτs∂αs∂βτ-C16kJβ/αkτs∂αs∂ατ(Kuuτsk)21=-C26kJα/βkτs∂βs∂βτ-C66kJkτs∂ατ∂βs(Kuuτsk)21=+C45k(Jαβkτzsz-1RβkJαkτsz1111111111111111-1RαkJβkτzs+1Rαk1RβkJkτs),(Kuuτsk)22=-C22kJα/βkτs∂βs∂βτ-C26kJkτs∂αs∂βτ(Kuuτsk)22=-C26kJkτs∂ατ∂βs-C66kJβ/αkτs∂αs∂ατ(Kuuτsk)22=+C44k(Jαβkτzsz-1RβkJαkτzs1111111111111111-1RβkJαkτsz+1Rβk2Jα/βkτs),(Kuuτsk)23=-C12k1RαkJkτs∂βτ-C22k1RβkJα/βkτs∂βτ(Kuuτsk)23=-C23kJαkτsz∂βτ-C16k1RαkJβ/αkτs∂ατ(Kuuτsk)23=-C26k1RβkJkτs∂ατ-C36kJβkτsz∂ατ(Kuuτsk)23=+C45k(Jβkτzs∂αs-1RβkJkτs∂αs)(Kuuτsk)23=+C44k(Jαkτzs∂βs-1RβkJα/βkτs∂βs),(Kuuτsk)31=C11k1RαkJβ/αkτs∂αs+C12k1RβkJkτs∂αs(Kuuτsk)31=+C13kJβkτzs∂αs+C16k1RαkJkτs∂βs(Kuuτsk)31=+C26k1RβkJα/βkτs∂βs+C36kJαkτzs∂βs(Kuuτsk)31=-C45k(Jαkτsz∂βτ-1RαkJkτs∂βτ)(Kuuτsk)31=-C55k(Jβkτsz∂ατ-1RαkJβ/αkτs∂ατ),(Kuuτsk)32=C12k1RαkJkτs∂βs+C22k1RβkJα/βkτs∂βs(Kuuτsk)32=+C23kJαkτzs∂βs+C16k1RαkJβ/αkτs∂αs(Kuuτsk)32=+C26k1RβkJkτs∂αs+C36kJβkτzs∂αs(Kuuτsk)32=-C45k(Jβkτsz∂ατ-1RβkJkτs∂ατ)(Kuuτsk)32=-C44k(Jαkτsz∂βτ-1RβkJα/βkτs∂βτ),(Kuuτsk)33=C11k1Rαk2Jβ/αkτs+C22k1Rβk2Jα/βkτs(Kuuτsk)33=+C33kJαβkτzsz+2C12k1Rαk1RβkJkτs(Kuuτsk)33=+C13k1Rαk(Jβkτzs+Jβkτsz)(Kuuτsk)33=+C23k1Rβk(Jαkτzs+Jαkτsz)(Kuuτsk)33=-C44kJα/βkτs∂βs∂βτ-C55kJβ/αkτs∂αs∂ατ(Kuuτsk)33=-C45kJkτs∂αs∂βτ-C45kJkτs∂ατ∂βs.The application of boundary conditions makes use of the fundamental nucleusΠd in the form:
(27)(Πuuτsk)11=nαC11kJβ/αkτs∂αs+nβC66kJα/βkτs∂βs(Πuuτsk)11=+nβC16kJkτs∂αs+nαC16kJkτs∂βs,(Πuuτsk)12=nαC16kJβ/αkτs∂αs+nβC26kJα/βkτs∂βs(Πuuτsk)12=+nαC12kJkτs∂βs+nβC66kJkτs∂αs,(Πuuτsk)13=nα1RαkC11kJβ/αkτs+nα1RβkC12kJkτs(Πuuτsk)13=+nαC13kJβkτsz+nβ1RαkC16kJkτs(Πuuτsk)13=+nβ1RβkC26kJαβkτs+nβC36kJαkτsz,(Πuuτsk)21=nαC16kJβ/αkτs∂αs+nβC26kJα/βkτs∂βs(Πuuτsk)21=+nβC12kJkτs∂αs+nαC66kJkτs∂βs,(Πuuτsk)22=nαC66kJβ/αkτs∂αs+nβC22kJα/βkτs∂βs(Πuuτsk)22=+nβC26kJkτs∂αs+nαC26kJkτs∂βs,(Πuuτsk)23=nα1RαkC16kJβ/αkτs+nα1RβkC26kJkτs(Πuuτsk)23=+nαC36kJβkτsz+nβ1RαkC12kJkτs(Πuuτsk)23=+nβ1RβkC22kJα/βkτs+nβC23kJαkτsz,(Πuuτsk)31=-nα1RαkC55kJβ/αkτs+nαC55kJβkτsz(Πuuτsk)31=-nβ1RαkC45kJkτs+nβC45kJαkτsz,(Πuuτsk)32=-nα1RβkC45kJkτs+nαC45kJβkτsz(Πuuτsk)32=-nβ1RβkC44kJα/βkτs+nβC44kJαkτsz,(Πuuτsk)33=nαC55kJβ/αkτs∂αs+nβC44kJα/βkτs∂βs(Πuuτsk)33=+nβC45kJkτs∂αs+nαC45kJkτs∂βs.One can note that all the equations written for the shell degenerate into those for the plate when1/Rαk=1/Rβk=0. In practice we set the radii of curvature to 109.


### 2.4. Dynamic Governing Equations

The PVD for the dynamic case is expressed as(28)∑k=1Nl‍∫Ωk‍∫Ak‍{δϵpGkTσpCk+δϵnGkTσnCk}dΩkdz=∑k=1Nl‍∫Ωk‍∫Ak‍ρkδukTu¨kdΩkdz+∑k=1Nl‍δLek,
where ρk is the mass density of the kth layer and double dots denote acceleration.By substituting the geometrical relations, the constitutive equations, and the Unified Formulation, we obtain the following governing equations:(29)δuskT:Kuukτsuτk=Mkτsu¨τk+Puτk.In the case of free vibrations one has(30)δuskT:Kuukτsuτk=Mkτsu¨τk,
where Mkτs is the fundamental nucleus for the inertial term. The explicit form of that is
(31)M11kτs=ρkJαβkτs,M12kτs=0,M13kτs=0,M21kτs=0,M22kτs=ρkJαβkτs,M23kτs=0,M31kτs=0,M32kτs=0,M33kτs=ρkJαβkτs,
where the meaning of the integral Jαβkτs 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.


## 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, withz being the thickness coordinate of the whole plate while zk is the layer thickness coordinate; see Figure 1. a and h are length and thickness of the square laminated plate, respectively. The adimensional layer coordinate ζk=(2zk)/hk is further introduced (hk is the thickness of the kth layer). Murakami's zig-zag function Z(z) was defined according to the following formula [11]:
(1)Z(z)=(-1)kζz.Z
(
z
) has the following properties.  (1)
It is a piecewise linear function of layer coordinateszk.(2)
Z(z) has unit amplitude for the whole layers.(3)
The slopeZ′(z)=dZ/dz 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:
(2)u=u0+zu1+(-1)k2hk(z-12(zk+zk+1))uZ,v=v0+zv1+(-1)k2hk(z-12(zk+zk+1))vZ,w=w0.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:(3)u=u0+zu1+(-1)k2hk(z-12(zk+zk+1))uZ,(4)v=v0+zv1+(-1)k2hk(z-12(zk+zk+1))vZ,(5)w=w0+zw1+(-1)k2hk(z-12(zk+zk+1))wZ,wherezk and zk+1 are the bottom and top z-coordinates at each layer. The additional degrees of freedom uZ and vZ 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 inz 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
(6)dsk2=Hαk2dα2+Hβk2dβ2+Hzk2dz2,dΩk=HαkHβkdαdβ,dV=HαkHβkHzkdαdβdz,where the metric coefficients are(7)Hαk=Ak(1+zRαk),Hβk=Bk(1+zRβk),Hzk=1.kdenotes the k-layer of the multilayered shell; Rαk and Rβk are the principal radii of curvature along the coordinates α and β, respectively. Ak and Bk are the coefficients of the first fundamental form of Ωk (Γk is the Ωk boundary). In this work, the attention has been restricted to shells with constant radii of curvature (cylindrical, spherical, and toroidal geometries) for which Ak=Bk=1.Although one can use the UF for one-layer, isotropic shell, a multilayered shell withNl layers is considered. The Principle of Virtual Displacement (PVD) for the pure-mechanical case reads
(8)∑k=1Nl‍∫Ωk‍∫Ak‍{δϵpGkTσpCk+δϵnGkTσnCk}×dΩkdz=∑k=1Nl‍δLek,
where Ωk and Ak are the integration domains in plane (α,β) and z direction, respectively. Here, k indicates the layer and T the transpose of a vector, and δLek is the external work for the kth layer. G means geometrical relations and C constitutive equations.The steps to obtain the governing equations are(i)
substitution of the geometrical relations (subscriptG),(ii)
substitution of the appropriate constitutive equations (subscriptC),(iii)
introduction of the Unified Formulation.Stresses and strains are separated into in-plane and normal components, denoted, respectively, by the subscriptsp and n. The mechanical strains in the kth layer can be related to the displacement field uk={uαk,uβk,uzk} via the geometrical relations:
(9)ϵpGk=[ϵααk,ϵββk,ϵαβk]T=(Dpk+Apk)uk,ϵnGk=[ϵαzk,ϵβzk,ϵzzk]T=(DnΩk+Dnzk-Ank)uk.The explicit form of the introduced arrays is as follows(10)Dpk=[∂αHαk000∂βHβk0∂βHβk∂αHαk0],DnΩk=[00∂αHαk00∂βHβk000],Dnzk=[∂z000∂z000∂z],(11)Apk=[001HαkRαk001HβkRβk000],Ank=[1HαkRαk0001HβkRβk0000].The 3D constitutive equations are given as(12)σpCk=CppkϵpGk+CpnkϵnGk,σnCk=CnpkϵpGk+CnnkϵnGk
with
(13)Cppk=[C11kC12kC16kC12kC22kC26kC16kC26kC66k],Cpnk=[00C13k00C23k00C36k],Cnpk=[000000C13kC23kC36k],Cnnk=[C55kC45k0C45kC44k000C33k].According to the Unified Formulation by Carrera, the three displacement componentsuα, uβ, and uz and their relative variations can be modelled as
(14)(uα,uβ,uz)=Fτ(uατ,uβτ,uzτ),(δuα,δuβ,δuz)=Fs(δuαs,δuβs,δuzs)
with Taylor expansions from the first up to the 4th order:  F0=z0=1, F1=z1=z,…, FN=zN,…, and F4=z4 if an Equivalent Single Layer (ESL) approach is used.Resorting to the displacement field in (3), we choose vectors Ft=[1z(-1)k(2/hk)(z-(1/2)(zk+zk+1))] for displacement u, v, and w. 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 ofFt, in order to solve this type of problems. We designed a MATLAB code that just by changing Ft can analyse static deformations, free vibrations, and buckling loads for any type of C∘ 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 Figure2 the assembling procedures are shown on layer k for ESL approach.Figure 2
Assembling procedure for ESL approach.Geometry and notations for a multilayered shell (doubly curved).
(a)
(b)Substituting the geometrical relations, the constitutive equations, and the Unified Formulation into the variational statement PVD, for thekth layer, one has
(15)∑k=1Nl‍{∫Ωk‍∫Ak‍{((Dp+Ap)δuk)T11111111111×(Cppk(Dp+Ap)uk111111111111111+Cpnk(DnΩ+Dnz-An)uk)11111111111+((DnΩ+Dnz-An)δuk)T11111111111×(Cnpk(Dp+Ap)uk111111111111111+Cnnk(DnΩ+Dnz-An)uk){((Dp+Ap)δuk)T}dΩkdzk{∫Ωk‍∫Ak‍{((Dp+Ap)δuk)T}11111111111=∑k=1Nl‍δLek.At this point, the formula of integration by parts is applied(16)∫Ωk‍((DΩ)δak)TakdΩk=-∫Ωk‍δakT((DΩT)ak)dΩk∫Ωk‍((DΩ)δak)TakdΩk=+∫Γk‍δakT((IΩ)ak)dΓk,
where IΩ matrix is obtained applying the gradient theorem(17)∫Ω‍∂ψ∂xidυ=∮Γ‍niψds
with ni being the components of the normal n^ to the boundary along the direction i. After integration by parts and the substitution of CUF, the governing equations and boundary conditions for the shell in the mechanical case are obtained:
(18)∑k=1Nl‍{∫Ωk‍∫Ak‍{δuskT[(-Dp+Ap)T11111111111111111×Fs(Cppk(Dp+Ap)Fτuτk1111111111111111111111+Cpnk(DnΩ+Dnz-An)Fτuτk)[(-Dp+Ap)T]111111111111+δuskT[+Cnnk(DnΩ+Dnz-An)Fτuτk)](-DnΩ+Dnz-An)T11111111111111111111×Fs(Cnpk(Dp+Ap)Fτuτk1111111111111111111111111+Cnnk(DnΩ+Dnz1111111111111111111111111-An)Fτuτk)]{δuskT[(-Dp+Ap)T}dΩkdzk{∫Ωk‍∫Ak‍{δuskT[(-Dp+Ap)T}+∑k=1Nl‍{∫Γk‍∫Ak‍{×Fτuτk)]}δuskT[IpTFs(Cppk(Dp+Ap)Fτuτk1111111111111111111111111+Cpnk(DnΩ+Dnz11111111111111111111111111111111-An)Fτuτk)]111111111111111+δuskT[InpTFs(Cnpk(Dp-Ap)Fτuτk111111111111111111111111111+Cnnk(DnΩ+Dnz-An)111111111111111111111111111×Fτuτk)]}dΓkdzk{∫Γk‍∫Ak‍{δuskT[IpTFs(Cppk(Dp+Ap)Fτuτk}=∑k=1Nl‍{∫ΩkδuskTFspuk},
where Ipk and Inpk depend on the boundary geometry:
(19)Ip=[nαHα000nβHβ0nβHβnαHα0],Inp=[00nαHα00nβHβ000].The normal to the boundary of domainΩ is
(20)n^=[nαnβ]=[cos(φα)cos(φβ)],
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(21)δuskT:Kuukτsuτk=Puτk,
where the fundamental nucleus Kuukτs is obtained as
(22)Kuukτs=∫Ak‍[[-Dp+Ap]TCppk[Dp+Ap]1111111111+[-Dp+Ap]TCpnk[DnΩ+Dnz-An]1111111111+[-DnΩ+Dnz-An]TCnpk[Dp+Ap]1111111111+[-DnΩ+Dnz-An]T1111111111×Cnnk[DnΩ+Dnz-An][[-Dp+Ap]TCppk[Dp+Ap]]FτFsHαkHβkdz.
And the corresponding Neumann-type boundary conditions on Γk are
(23)Πdkτsuτk=Πdkτsu-τk,
where
(24)Πdkτs=∫Ak‍[×[DnΩ+Dnz-Anτ]]IpTCppk[Dp+Apτ]+IpTCpnk[DnΩ+Dnz-Anτ]1111111111+InpTCnpk[Dp+Apτ]+InpTCnnk1111111111×[DnΩ+Dnz-Anτ][×[DnΩ+Dnz-Anτ]]IpTCppk[Dp+Apτ]+IpTCpnk[DnΩ+Dnz-Anτ]]FτFsHαkHβkdz.
And Puτk are variationally consistent loads with applied pressure.


## 2.3. Fundamental Nuclei

The following integrals are introduced to perform the explicit form of fundamental nuclei:(25)(Jkτs,Jαkτs,Jβkτs,Jα/βkτs,Jβ/αkτs,Jαβkτs)=∫Ak‍FτFs(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτzs,Jαkτzs,Jβkτzs,Jα/βkτzs,Jβ/αkτzs,Jαβkτzs)=∫Ak‍∂Fτ∂zFs(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτsz,Jαkτsz,Jβkτsz,Jα/βkτsz,Jβ/αkτsz,Jαβkτsz)=∫Ak‍Fτ∂Fs∂z(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz,(Jkτzsz,Jαkτzsz,Jβkτzsz,Jα/βkτzsz,Jβ/αkτzsz,Jαβkτzsz)=∫Ak‍∂Fτ∂z∂Fs∂z(1,Hα,Hβ,HαHβ,HβHα,HαHβ)dz.The fundamental nucleusKuukτs is reported for doubly curved shells (radii of curvature in both α and β directions; see Figure 3):
(26)(Kuuτsk)11=-C11kJβ/αkτs∂αs∂ατ-C16kJkτs∂ατ∂βs(Kuuτsk)11=-C16kJkτs∂αs∂βτ-C66kJα/βkτs∂βs∂βτ(Kuuτsk)11=+C55k(Jαβkτzsz-1RαkJβkτzs1111111111111111-1RαkJβkτsz+1Rαk2Jβ/αkτs),(Kuuτsk)12=-C12kJkτs∂ατ∂βs-C16kJβ/αkτs∂αs∂ατ(Kuuτsk)12=-C26kJα/βkτs∂βs∂βτ-C66kJkτs∂αs∂βτ(Kuuτsk)12=+C45k(Jαβkτzsz-1RβkJαkτzs1111111111111111-1RαkJβkτsz+1Rαk1RβkJkτs),(Kuuτsk)13=-C11k1RαkJβ/αkτs∂ατ-C12k1RβkJkτs∂ατ(Kuuτsk)13=-C13kJβkτsz∂ατ-C16k1RαkJkτs∂βτ(Kuuτsk)13=-C26k1RβkJα/βkτs∂βτ-C36kJαkτsz∂βτ(Kuuτsk)13=+C45k(Jαkτzs∂βs-1RαkJkτs∂βs)(Kuuτsk)13=+C55k(Jβkτzs∂αs-1RαkJβ/αkτs∂αs),(Kuuτsk)21=-C12kJkτs∂αs∂βτ-C16kJβ/αkτs∂αs∂ατ(Kuuτsk)21=-C26kJα/βkτs∂βs∂βτ-C66kJkτs∂ατ∂βs(Kuuτsk)21=+C45k(Jαβkτzsz-1RβkJαkτsz1111111111111111-1RαkJβkτzs+1Rαk1RβkJkτs),(Kuuτsk)22=-C22kJα/βkτs∂βs∂βτ-C26kJkτs∂αs∂βτ(Kuuτsk)22=-C26kJkτs∂ατ∂βs-C66kJβ/αkτs∂αs∂ατ(Kuuτsk)22=+C44k(Jαβkτzsz-1RβkJαkτzs1111111111111111-1RβkJαkτsz+1Rβk2Jα/βkτs),(Kuuτsk)23=-C12k1RαkJkτs∂βτ-C22k1RβkJα/βkτs∂βτ(Kuuτsk)23=-C23kJαkτsz∂βτ-C16k1RαkJβ/αkτs∂ατ(Kuuτsk)23=-C26k1RβkJkτs∂ατ-C36kJβkτsz∂ατ(Kuuτsk)23=+C45k(Jβkτzs∂αs-1RβkJkτs∂αs)(Kuuτsk)23=+C44k(Jαkτzs∂βs-1RβkJα/βkτs∂βs),(Kuuτsk)31=C11k1RαkJβ/αkτs∂αs+C12k1RβkJkτs∂αs(Kuuτsk)31=+C13kJβkτzs∂αs+C16k1RαkJkτs∂βs(Kuuτsk)31=+C26k1RβkJα/βkτs∂βs+C36kJαkτzs∂βs(Kuuτsk)31=-C45k(Jαkτsz∂βτ-1RαkJkτs∂βτ)(Kuuτsk)31=-C55k(Jβkτsz∂ατ-1RαkJβ/αkτs∂ατ),(Kuuτsk)32=C12k1RαkJkτs∂βs+C22k1RβkJα/βkτs∂βs(Kuuτsk)32=+C23kJαkτzs∂βs+C16k1RαkJβ/αkτs∂αs(Kuuτsk)32=+C26k1RβkJkτs∂αs+C36kJβkτzs∂αs(Kuuτsk)32=-C45k(Jβkτsz∂ατ-1RβkJkτs∂ατ)(Kuuτsk)32=-C44k(Jαkτsz∂βτ-1RβkJα/βkτs∂βτ),(Kuuτsk)33=C11k1Rαk2Jβ/αkτs+C22k1Rβk2Jα/βkτs(Kuuτsk)33=+C33kJαβkτzsz+2C12k1Rαk1RβkJkτs(Kuuτsk)33=+C13k1Rαk(Jβkτzs+Jβkτsz)(Kuuτsk)33=+C23k1Rβk(Jαkτzs+Jαkτsz)(Kuuτsk)33=-C44kJα/βkτs∂βs∂βτ-C55kJβ/αkτs∂αs∂ατ(Kuuτsk)33=-C45kJkτs∂αs∂βτ-C45kJkτs∂ατ∂βs.The application of boundary conditions makes use of the fundamental nucleusΠd in the form:
(27)(Πuuτsk)11=nαC11kJβ/αkτs∂αs+nβC66kJα/βkτs∂βs(Πuuτsk)11=+nβC16kJkτs∂αs+nαC16kJkτs∂βs,(Πuuτsk)12=nαC16kJβ/αkτs∂αs+nβC26kJα/βkτs∂βs(Πuuτsk)12=+nαC12kJkτs∂βs+nβC66kJkτs∂αs,(Πuuτsk)13=nα1RαkC11kJβ/αkτs+nα1RβkC12kJkτs(Πuuτsk)13=+nαC13kJβkτsz+nβ1RαkC16kJkτs(Πuuτsk)13=+nβ1RβkC26kJαβkτs+nβC36kJαkτsz,(Πuuτsk)21=nαC16kJβ/αkτs∂αs+nβC26kJα/βkτs∂βs(Πuuτsk)21=+nβC12kJkτs∂αs+nαC66kJkτs∂βs,(Πuuτsk)22=nαC66kJβ/αkτs∂αs+nβC22kJα/βkτs∂βs(Πuuτsk)22=+nβC26kJkτs∂αs+nαC26kJkτs∂βs,(Πuuτsk)23=nα1RαkC16kJβ/αkτs+nα1RβkC26kJkτs(Πuuτsk)23=+nαC36kJβkτsz+nβ1RαkC12kJkτs(Πuuτsk)23=+nβ1RβkC22kJα/βkτs+nβC23kJαkτsz,(Πuuτsk)31=-nα1RαkC55kJβ/αkτs+nαC55kJβkτsz(Πuuτsk)31=-nβ1RαkC45kJkτs+nβC45kJαkτsz,(Πuuτsk)32=-nα1RβkC45kJkτs+nαC45kJβkτsz(Πuuτsk)32=-nβ1RβkC44kJα/βkτs+nβC44kJαkτsz,(Πuuτsk)33=nαC55kJβ/αkτs∂αs+nβC44kJα/βkτs∂βs(Πuuτsk)33=+nβC45kJkτs∂αs+nαC45kJkτs∂βs.One can note that all the equations written for the shell degenerate into those for the plate when1/Rαk=1/Rβk=0. In practice we set the radii of curvature to 109.


## 2.4. Dynamic Governing Equations

The PVD for the dynamic case is expressed as(28)∑k=1Nl‍∫Ωk‍∫Ak‍{δϵpGkTσpCk+δϵnGkTσnCk}dΩkdz=∑k=1Nl‍∫Ωk‍∫Ak‍ρkδukTu¨kdΩkdz+∑k=1Nl‍δLek,
where ρk is the mass density of the kth layer and double dots denote acceleration.By substituting the geometrical relations, the constitutive equations, and the Unified Formulation, we obtain the following governing equations:(29)δuskT:Kuukτsuτk=Mkτsu¨τk+Puτk.In the case of free vibrations one has(30)δuskT:Kuukτsuτk=Mkτsu¨τk,
where Mkτs is the fundamental nucleus for the inertial term. The explicit form of that is
(31)M11kτs=ρkJαβkτs,M12kτs=0,M13kτs=0,M21kτs=0,M22kτs=ρkJαβkτs,M23kτs=0,M31kτs=0,M32kτs=0,M33kτs=ρkJαβkτs,
where the meaning of the integral Jαβkτs 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 operatorL and a bounded region Ω in ℝn with some boundary ∂Ω. In the static problems we seek the computation of displacement (u) from the global system of equations
(32)ℒu=finΩ,ℒBu=gon∂Ω,
where ℒ and ℒB 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 (u) that satisfy
(33)ℒu+λu=0inΩ,ℒBu=0on∂Ω.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 (u) is given by
(34)u~(x)=∑i=1N‍αiϕ(∥x-yi∥2),x∈ℝn,
where yi,i=1,…,N, is a finite set of distinct points (centers) in ℝn. The most common RBFs are
(35)Cubic:ϕ(r)=r3,Thinplatesplines:ϕ(r)=r2log(r),Wendlandfunctions:ϕ(r)=(1-r)+mp(r),Gaussian:ϕ(r)=e-(cr)2,Multiquadrics:ϕ(r)=c2+r2,InverseMultiquadrics:ϕ(r)=(c2+r2)-1/2,
where the Euclidian distance r is real and nonnegative and c 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 N distinct interpolations and knowing u(xj), j=1,2,…,N, we find αi by the solution of a N×N linear system
(36)Aα=u,
where A=[ϕ(∥x-yi∥2)]N×N, α=[α1,α2,…,αN]T, and u=[u(x1),u(x2),…,u(xN)]T.


### 3.4. Solution of the Static Problem

The solution of a static problem by radial basis functions considersNI nodes in the domain and NB nodes on the boundary, with a total number of nodes N=NI+NB. We denote the sampling points by xi∈Ω, i=1,…,NI, and xi∈∂Ω, i=NI+1,…,N. At the points in the domain we solve the following system of equations:
(37)∑i=1N‍αiℒϕ(∥x-yi∥2)=f(xj),j=1,2,…,NI,
or
(38)ℒIα=F,
where
(39)ℒI=[ℒϕ(∥x-yi∥2)]NI×N.At the points on the boundary, we impose boundary conditions as(40)∑i=1N‍αiℒBϕ(∥x-yi∥2)=g(xj),j=NI+1,⋯,N,
or
(41)Bα=G,
where
(42)B=ℒBϕ[(∥xNI+1-yj∥2)]NB×N.Therefore, we can write a finite-dimensional static problem as(43)[ℒIB]α=[FG]By inverting the system (43), we obtain the vector α. We then obtain the solution u using the interpolation (34).


### 3.5. Solution of the Eigenproblem

We considerNI nodes in the interior of the domain and NB nodes on the boundary, with N=NI+NB. We denote interpolation points by xi∈Ω, i=1,…,NI, and xi∈∂Ω, i=NI+1,…,N. At the points in the domain, we define the eigenproblem as
(44)∑i=1N‍αiℒϕ(∥x-yi∥2)=λu~(xj),j=1,2,…,NI,
or
(45)ℒIα=λu~I,
where
(46)ℒI=[ℒϕ(∥x-yi∥2)]NI×N.At the points on the boundary, we enforce the boundary conditions as(47)∑i=1N‍αiℒBϕ(∥x-yi∥2)=0,j=NI+1,…,N,
or
(48)Bα=0.Equations (45) and (48) can now be solved as a generalized eigenvalue problem
(49)[ℒIB]α=λ[AI0]α,
where
(50)AI=ϕ[(∥xNI-yj∥2)]NI×N.


### 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
(51)α=[LIB]-1[FG].Thisα vector is then used to obtain solution u~, by using (5). If derivatives of u~ are needed, such derivatives are computed as
(52)∂u~∂x=∑j=1N‍αj∂ϕj∂x,∂2u~∂x2=∑j=1N‍αj∂2ϕj∂x2,and so forth.In the present collocation approach, we need to impose essential and natural boundary conditions. Consider, for example, the conditionw=0, on a simply supported or clamped edge. We enforce the conditions by interpolating as
(53)w=0⟶∑j=1N‍αjWϕj=0.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,u0,u1,v0,v1,…, as
(54)u0=U0(w,y)eiωt;u1=U1(w,y)eiωt,uZ=UZ(w,y)eiωt,v0=V0(w,y)eiωt,v1=V1(w,y)eiωt,vZ=VZ(w,y)eiωt,w0=W0(w,y)eiωt,w1=W1(w,y)eiωt,w2=W2(w,y)eiωt,
where ω is the frequency of natural vibration. Substituting the harmonic expansion into (49) in terms of the amplitudes U0, U1, UZ, V0, V1, VZ, W0, W1, and W2, we may obtain the natural frequencies and vibration modes for the plate or shell problem, by solving the eigenproblem
(55)[ℒ-ω2𝒢]X=0,
where ℒ collects all stiffness terms and 𝒢 collects all terms related to the inertial terms. In (55) X are the modes of vibration associated with the natural frequencies defined as ω.


## 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 operatorL and a bounded region Ω in ℝn with some boundary ∂Ω. In the static problems we seek the computation of displacement (u) from the global system of equations
(32)ℒu=finΩ,ℒBu=gon∂Ω,
where ℒ and ℒB 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 (u) that satisfy
(33)ℒu+λu=0inΩ,ℒBu=0on∂Ω.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 (u) is given by
(34)u~(x)=∑i=1N‍αiϕ(∥x-yi∥2),x∈ℝn,
where yi,i=1,…,N, is a finite set of distinct points (centers) in ℝn. The most common RBFs are
(35)Cubic:ϕ(r)=r3,Thinplatesplines:ϕ(r)=r2log(r),Wendlandfunctions:ϕ(r)=(1-r)+mp(r),Gaussian:ϕ(r)=e-(cr)2,Multiquadrics:ϕ(r)=c2+r2,InverseMultiquadrics:ϕ(r)=(c2+r2)-1/2,
where the Euclidian distance r is real and nonnegative and c 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 N distinct interpolations and knowing u(xj), j=1,2,…,N, we find αi by the solution of a N×N linear system
(36)Aα=u,
where A=[ϕ(∥x-yi∥2)]N×N, α=[α1,α2,…,αN]T, and u=[u(x1),u(x2),…,u(xN)]T.


## 3.4. Solution of the Static Problem

The solution of a static problem by radial basis functions considersNI nodes in the domain and NB nodes on the boundary, with a total number of nodes N=NI+NB. We denote the sampling points by xi∈Ω, i=1,…,NI, and xi∈∂Ω, i=NI+1,…,N. At the points in the domain we solve the following system of equations:
(37)∑i=1N‍αiℒϕ(∥x-yi∥2)=f(xj),j=1,2,…,NI,
or
(38)ℒIα=F,
where
(39)ℒI=[ℒϕ(∥x-yi∥2)]NI×N.At the points on the boundary, we impose boundary conditions as(40)∑i=1N‍αiℒBϕ(∥x-yi∥2)=g(xj),j=NI+1,⋯,N,
or
(41)Bα=G,
where
(42)B=ℒBϕ[(∥xNI+1-yj∥2)]NB×N.Therefore, we can write a finite-dimensional static problem as(43)[ℒIB]α=[FG]By inverting the system (43), we obtain the vector α. We then obtain the solution u using the interpolation (34).


## 3.5. Solution of the Eigenproblem

We considerNI nodes in the interior of the domain and NB nodes on the boundary, with N=NI+NB. We denote interpolation points by xi∈Ω, i=1,…,NI, and xi∈∂Ω, i=NI+1,…,N. At the points in the domain, we define the eigenproblem as
(44)∑i=1N‍αiℒϕ(∥x-yi∥2)=λu~(xj),j=1,2,…,NI,
or
(45)ℒIα=λu~I,
where
(46)ℒI=[ℒϕ(∥x-yi∥2)]NI×N.At the points on the boundary, we enforce the boundary conditions as(47)∑i=1N‍αiℒBϕ(∥x-yi∥2)=0,j=NI+1,…,N,
or
(48)Bα=0.Equations (45) and (48) can now be solved as a generalized eigenvalue problem
(49)[ℒIB]α=λ[AI0]α,
where
(50)AI=ϕ[(∥xNI-yj∥2)]NI×N.


## 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
(51)α=[LIB]-1[FG].Thisα vector is then used to obtain solution u~, by using (5). If derivatives of u~ are needed, such derivatives are computed as
(52)∂u~∂x=∑j=1N‍αj∂ϕj∂x,∂2u~∂x2=∑j=1N‍αj∂2ϕj∂x2,and so forth.In the present collocation approach, we need to impose essential and natural boundary conditions. Consider, for example, the conditionw=0, on a simply supported or clamped edge. We enforce the conditions by interpolating as
(53)w=0⟶∑j=1N‍αjWϕj=0.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,u0,u1,v0,v1,…, as
(54)u0=U0(w,y)eiωt;u1=U1(w,y)eiωt,uZ=UZ(w,y)eiωt,v0=V0(w,y)eiωt,v1=V1(w,y)eiωt,vZ=VZ(w,y)eiωt,w0=W0(w,y)eiωt,w1=W1(w,y)eiωt,w2=W2(w,y)eiωt,
where ω is the frequency of natural vibration. Substituting the harmonic expansion into (49) in terms of the amplitudes U0, U1, UZ, V0, V1, VZ, W0, W1, and W2, we may obtain the natural frequencies and vibration modes for the plate or shell problem, by solving the eigenproblem
(55)[ℒ-ω2𝒢]X=0,
where ℒ collects all stiffness terms and 𝒢 collects all terms related to the inertial terms. In (55) X 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(56)ϕ(r)=(1-cr)8+(32(cr)3+25(cr)2+8cr+1),
where the shape parameter (c) 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 sidea and thickness h, to be composed of layers oriented at [0∘/90∘/0∘] and [0∘/90∘/90∘/0∘]. The shell is subjected to a sinusoidal vertical pressure of the form
(57)pz=Psin(πxa)sin(πya)
with the origin of the coordinate system located at the lower left corner on the midplane and P the maximum load (at center of shell).The orthotropic material properties for each layer are given by(58)E1=25.0E2,G12=G13=0.5E2,G23=0.2E2,ν12=0.25.The in-plane displacement, the transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(59)2w¯=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,τ¯xy=τxy(0,0,h/2)h2Pa2.The shell is simply supported on all edges.In Table1 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 R/a and two values of a/h (10 and 100). Results are in good agreement for various a/h ratios with the higher-order results of Reddy and Liu [56].Table 1
Nondimensional central deflection,w¯=w(102E2h3/P0a4), and variation with various numbers of grid points per unit length, N, for different R/a ratios, for R1=R2.
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
[
0
∘
/
9
0
∘
/
0
∘
]
10
Present (13×13)
6.8510
7.0818
7.1429
7.1609
7.1637
7.1649
10
Present (17×17)
6.8516
7.0822
7.1433
7.1612
7.1640
7.1653
10
Present (21×21)
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 (13×13)
1.0245
2.3658
3.5167
4.0714
4.1652
4.1975
100
Present (17×17)
1.0250
2.3667
3.5181
4.0728
4.1667
4.1990
100
Present (21×21)
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
[
0
∘
/
9
0
∘
/
9
0
∘
/
0
∘
]
10
Present (13×13)
6.7737
7.0012
7.0614
7.0791
7.0819
7.0831
10
Present (17×17)
6.7742
7.0015
7.0618
7.0795
7.0822
7.0834
10
Present (21×21)
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 (13×13)
1.0190
2.3583
3.5125
4.0702
4.1647
4.1972
100
Present (17×17)
1.0194
2.3593
3.5138
4.0717
4.1662
4.1987
100
Present (21×21)
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 with13×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 2
Nondimensionalized fundamental frequencies of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, and laminate ([0∘/90∘/90∘/0∘]).
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
10
Present (13×13)
12.0527
11.8889
11.8474
11.8357
11.8340
11.8335
Present (17×17)
12.0523
11.8886
11.8471
11.8355
11.8338
11.8332
Present (21×21)
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 (13×13)
31.2170
20.5742
16.8698
15.6744
15.4960
15.4361
Present (17×17)
31.2072
20.5679
16.8648
15.6698
15.4915
15.4316
Present (21×21)
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.170Table 3
Nondimensionalized fundamental frequencies of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, and laminate ([0∘/90∘/0∘]).
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
10
Present (13×13)
11.9831
11.8192
11.7777
11.7660
11.7643
11.7638
Present (17×17)
11.9827
11.8190
11.7774
11.7658
11.7641
11.7635
Present (21×21)
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 (13×13)
31.1343
20.5420
16.8595
15.6721
15.4950
15.4355
Present (17×17)
31.1244
20.5357
16.8545
15.6675
15.4905
15.4310
Present (21×21)
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.170Table4 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.Table 4
Nondimensionalized fundamental frequencies of cross-ply cylindrical shells,ω¯=ω(a2/h)ρ/E2.
R
/
a
Method
[0/90/0]
[0/90/90/0]
a
/
h
=
100
a
/
h
=
10
a
/
h
=
100
a
/
h
=
10
5
Present (13×13)
20.4956
11.7774
20.5298
11.8524
Present (17×17)
20.4839
11.7770
20.5201
11.8521
Present (21×21)
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 (13×13)
16.8475
11.7672
16.8583
11.8382
Present (17×17)
16.8409
11.7669
16.8522
11.8380
Present (21×21)
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 (13×13)
15.8006
11.7646
15.8039
11.8346
Present (17×17)
15.7955
11.7644
15.7990
11.8344
Present (21×21)
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 (13×13)
15.4945
11.7639
15.4955
11.8336
Present (17×17)
15.4899
11.7637
15.4909
11.8334
Present (21×21)
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 (13×13)
15.4503
11.7638
15.4510
11.8335
Present (17×17)
15.4457
11.7636
15.4464
11.8333
Present (21×21)
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 (13×13)
15.4355
11.7638
15.4361
11.8335
Present (17×17)
15.4310
11.7635
15.4316
11.8332
Present (21×21)
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.780In Figure4 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, ω¯=ω(a2/h)ρ/E2, for a laminate ([0∘/90∘/90∘/0∘]), using a grid of 21×21 points, for a/h=100, R/a=10. The modes of vibration are quite stable.The first 4 vibrational modes of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, for laminate ([0∘/90∘/90∘/0∘]) using a grid of 21×21 points, for a/h=100, R/a=10.
(a)
(b)
(c)
(d)In Figure5 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, ω¯=ω(a2/h)ρ/E2, for a laminate ([0∘/90∘/0∘]), using a grid of 21×21 points, for a/h=100, R/a=10e9. Again the modes of vibration are quite stable.The first 4 vibrational modes of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, for laminate ([0∘/90∘/0∘]) using a grid of 21×21 points, for a/h=100, R/a=10e9 (plate).
(a)
(b)
(c)
(d)


## 4.1. Spherical Shell in Bending

A laminated composite spherical shell is here considered, of sidea and thickness h, to be composed of layers oriented at [0∘/90∘/0∘] and [0∘/90∘/90∘/0∘]. The shell is subjected to a sinusoidal vertical pressure of the form
(57)pz=Psin(πxa)sin(πya)
with the origin of the coordinate system located at the lower left corner on the midplane and P the maximum load (at center of shell).The orthotropic material properties for each layer are given by(58)E1=25.0E2,G12=G13=0.5E2,G23=0.2E2,ν12=0.25.The in-plane displacement, the transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(59)2w¯=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,τ¯xy=τxy(0,0,h/2)h2Pa2.The shell is simply supported on all edges.In Table1 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 R/a and two values of a/h (10 and 100). Results are in good agreement for various a/h ratios with the higher-order results of Reddy and Liu [56].Table 1
Nondimensional central deflection,w¯=w(102E2h3/P0a4), and variation with various numbers of grid points per unit length, N, for different R/a ratios, for R1=R2.
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
[
0
∘
/
9
0
∘
/
0
∘
]
10
Present (13×13)
6.8510
7.0818
7.1429
7.1609
7.1637
7.1649
10
Present (17×17)
6.8516
7.0822
7.1433
7.1612
7.1640
7.1653
10
Present (21×21)
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 (13×13)
1.0245
2.3658
3.5167
4.0714
4.1652
4.1975
100
Present (17×17)
1.0250
2.3667
3.5181
4.0728
4.1667
4.1990
100
Present (21×21)
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
[
0
∘
/
9
0
∘
/
9
0
∘
/
0
∘
]
10
Present (13×13)
6.7737
7.0012
7.0614
7.0791
7.0819
7.0831
10
Present (17×17)
6.7742
7.0015
7.0618
7.0795
7.0822
7.0834
10
Present (21×21)
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 (13×13)
1.0190
2.3583
3.5125
4.0702
4.1647
4.1972
100
Present (17×17)
1.0194
2.3593
3.5138
4.0717
4.1662
4.1987
100
Present (21×21)
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 with13×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 2
Nondimensionalized fundamental frequencies of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, and laminate ([0∘/90∘/90∘/0∘]).
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
10
Present (13×13)
12.0527
11.8889
11.8474
11.8357
11.8340
11.8335
Present (17×17)
12.0523
11.8886
11.8471
11.8355
11.8338
11.8332
Present (21×21)
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 (13×13)
31.2170
20.5742
16.8698
15.6744
15.4960
15.4361
Present (17×17)
31.2072
20.5679
16.8648
15.6698
15.4915
15.4316
Present (21×21)
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.170Table 3
Nondimensionalized fundamental frequencies of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, and laminate ([0∘/90∘/0∘]).
a
/
h
Method
R
/
a
5
10
20
50
100
1
0
9
10
Present (13×13)
11.9831
11.8192
11.7777
11.7660
11.7643
11.7638
Present (17×17)
11.9827
11.8190
11.7774
11.7658
11.7641
11.7635
Present (21×21)
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 (13×13)
31.1343
20.5420
16.8595
15.6721
15.4950
15.4355
Present (17×17)
31.1244
20.5357
16.8545
15.6675
15.4905
15.4310
Present (21×21)
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.170Table4 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.Table 4
Nondimensionalized fundamental frequencies of cross-ply cylindrical shells,ω¯=ω(a2/h)ρ/E2.
R
/
a
Method
[0/90/0]
[0/90/90/0]
a
/
h
=
100
a
/
h
=
10
a
/
h
=
100
a
/
h
=
10
5
Present (13×13)
20.4956
11.7774
20.5298
11.8524
Present (17×17)
20.4839
11.7770
20.5201
11.8521
Present (21×21)
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 (13×13)
16.8475
11.7672
16.8583
11.8382
Present (17×17)
16.8409
11.7669
16.8522
11.8380
Present (21×21)
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 (13×13)
15.8006
11.7646
15.8039
11.8346
Present (17×17)
15.7955
11.7644
15.7990
11.8344
Present (21×21)
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 (13×13)
15.4945
11.7639
15.4955
11.8336
Present (17×17)
15.4899
11.7637
15.4909
11.8334
Present (21×21)
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 (13×13)
15.4503
11.7638
15.4510
11.8335
Present (17×17)
15.4457
11.7636
15.4464
11.8333
Present (21×21)
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 (13×13)
15.4355
11.7638
15.4361
11.8335
Present (17×17)
15.4310
11.7635
15.4316
11.8332
Present (21×21)
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.780In Figure4 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, ω¯=ω(a2/h)ρ/E2, for a laminate ([0∘/90∘/90∘/0∘]), using a grid of 21×21 points, for a/h=100, R/a=10. The modes of vibration are quite stable.The first 4 vibrational modes of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, for laminate ([0∘/90∘/90∘/0∘]) using a grid of 21×21 points, for a/h=100, R/a=10.
(a)
(b)
(c)
(d)In Figure5 we illustrate the first 4 vibrational modes of cross-ply laminated spherical shells, ω¯=ω(a2/h)ρ/E2, for a laminate ([0∘/90∘/0∘]), using a grid of 21×21 points, for a/h=100, R/a=10e9. Again the modes of vibration are quite stable.The first 4 vibrational modes of cross-ply laminated spherical shells,ω¯=ω(a2/h)ρ/E2, for laminate ([0∘/90∘/0∘]) using a grid of 21×21 points, for a/h=100, R/a=10e9 (plate).
(a)
(b)
(c)
(d)


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



---

*Source: 123465-2013-11-11.xml*