Foam core composite sandwich plates and shells with variable stiffness: Effect of the curvilinear fiber path on the modal response

Abstract

This paper presents the free vibration analysis of composite sandwich plates and doubly curved shells with variable stiffness. The reinforcing fibers are located in the external skins of the sandwich structures according to curved paths. These curvilinear paths are described by a general expression that combines power-law, sinusoidal, exponential, Gaussian and ellipse-shaped functions. As a consequence, the reinforcing fibers are placed in these orthotropic layers in an arbitrary manner, in order to achieve the desired mechanical properties. The effect of this variable fiber orientation on the natural frequencies is investigated by means of several parametric studies. As far as the structural theory is concerned, an equivalent single layer approach based on the well-known Carrera Unified Formulation is employed. The Murakami’s function is added to the kinematic model to capture the zig-zag effect, when the soft-core effect is significant. Thus, several higher order shear deformation theories are taken into account in a unified manner. The differential geometry is employed to describe the reference surface of doubly curved shells and panels, which are characterized by variable radii of curvature. The numerical solution is obtained using the generalized differential quadrature method, due to its accuracy and stability features. The present solution is compared with the results available in the literature or obtained by finite element commercial codes.

Keywords

Composite sandwich plates and shells free vibration analysis curvilinear reinforcing fibers higher order equivalent single layer theories generalized differential quadrature

Introduction

In the last decades, the need of lightweight structures with high levels of strength has inevitably encouraged the development and the use of different kinds of composite materials in various engineering fields, such as in aerospace, automotive sector and boat hull manufacturing [1,2]. It is clear that the advancements of composite material production arose in order to obtain stiffer and lighter mediums with remarkable improvements in terms of static and dynamic behavior, if compared to the conventional materials, such as isotropic bodies. Among them, the class of laminated composites has achieved an impressive appeal in both manufacturing and academic research, as clearly stated in the works by Groh and Weaver [3 –6]. In these papers, the limitless potentialities of laminated composite structures have been highlighted. By superimposing layers with different mechanical properties, it is possible to obtain the most appropriate stacking sequence that adapts mainly to the structural problem under consideration. Thus, peculiar lamination schemes are designed to solve some issues related to the structural buckling, to achieve specific stress profiles and to modify the first natural frequencies. It is worth spending few words also on the class of functionally graded materials (FGMs). Due to their continuous gradual variation of the mechanical properties along a particular direction, which is commonly given by the structure thickness, these materials have been developed to solve the issues that affect laminated composites, such as stress concentrations or delamination [7 –10]. Nevertheless, reliable theories must be developed to analyze in a correct way the mechanical behavior of multilayered structures, especially if they are doubly curved shells. The development of the so-called Carrera Unified Formulation (CUF) has laid the foundations for the growth of enormous number of structural models for the study of laminated composite plates and shells [11]. As a consequence, several papers are available in the literature to prove that higher order structural theories based on CUF represent an excellent tool to analyze the mechanical behavior of layered structures [12 –26]. The inadequacy of first-order or classical shell theories has been underlined mainly in the study of sandwich structures [27 –30]. In general, these structural elements are made of three layers: two external stiff skins and a soft-core. This particular configuration generates some effects, such as the stretching and squeezing ones, that make unusable first-order structural theories as a matter of fact. However, accurate results have been obtained in these cases introducing peculiar functions to the kinematic model, such as the Murakami’s function, which are able to capture the so-called zig-zag effect. This aspect is investigated in depth in the works by Maturi et al. [31] and Carrera [32,33]. Recently, the well-known and hugely exploited fiber-reinforced layers, commonly used in the stacking sequence of laminated composite structures, have been subjected to such an innovative development that has led to a new class of composite materials. This particular kind of composite material is characterized by the reinforcing phase placed along curvilinear paths [34 –36]. The first works about this topic have been presented in the 1970s, as proven by the paper presented by Sendeckyj [37]. Then, important contributions for the development of the curvilinear reinforcing fibers have been given by Hyer and Charette [38] in the 1990s. Currently, the interest of many researchers is focused on the optimum design and optimization processes of composite structures with curved reinforcing fibers, as it can be noticed in the works by Parnas et al. [39] and by Blom et al. [40 –42]. In these works, as well as in the more recent papers by Rouhi et al. [43] and by Coburn and Weaver [44], the problems related to the manufacture of these composites are also investigated. As illustrated in depth in the works by Nik et al. [45] and by Akbarzadeh et al. [46], such composites present necessarily some defects. The most typical issues that affect variable stiffness composites are gaps and overlaps. For this purpose, the work by Kim et al. [47] must be mentioned, since an innovative production technology has been developed to overcome these problems. For further comments and details concerning the manufacturing of laminates reinforced by curvilinear fibers, the readers can refer to the papers just mentioned. The manufacture techniques are not taken into account in the present work, assuming that the presence of such defects is limited.

Nevertheless, it has been proven that a proper choice of the curvilinear paths of the fibers can produce significant improvements to face some problems related to the buckling [48 –50] and post-buckling behavior [51 –54] of this kind of laminated composite structures. Analogously, many papers available in the literature have proven the advantages of such composites in the free vibrations analysis, both in linear [55 –58] and non-linear [59 –61] hypotheses. In particular, as far as the linear field is considered, it should be mentioned that Akhavan and Ribeiro [55] have performed the linear modal analysis in order to evaluate the natural frequencies of laminated variable stiffness plates by using a third-order shear deformation theory. The same analysis is carried out by Yazdani and Ribeiro [56] employing a layer-wise formulation. In addition, Tornabene et al. [57] has extended the variable stiffness concept also to doubly curved shells reinforced by curvilinear fibers. Finally, it is worth noting that Honda and Narita [58] have studied the dynamic behavior of variable stiffness plates reinforced by arbitrarily shaped fibers. In particular, they have employed spline functions to describe the curved path of the fibers. The static behavior of variable stiffness structure has been investigated too. As proven by the papers [62 –66], the curvilinear placement of the fibers affects both the displacement and the stress and strain profiles of these variable stiffness structures. According to the authors’ knowledge, the paper by Vescovini and Dozio [67] represents the first example of combining together the effect of variable stiffness layers and the advantages of sandwich structures. In their work, they have performed both the free vibration and buckling analysis of sandwich plates whose external face-sheets are reinforced by curvilinear fibers. For completeness purposes, it should be mentioned that their solutions have been obtained by using thin plate theory and the Ritz method. The aim of the present paper is to illustrate a general mathematical formulation able to describe the curvilinear paths of the reinforcing fibers of the external orthotropic face-sheets. In particular, the effect of the curvilinear placement of the reinforcing fibers on the modal response of sandwich plates and shells with variable stiffness is investigated by means of several parametric analyses. Starting from the comparison with the results available in the literature [67], the free vibrations of variable stiffness composite sandwich plates and shells are obtained by using the well-known generalized differential quadrature (GDQ) method [68]. The main features of this numerical approach have been explained in the review paper by Tornabene et al. [69].

Shell geometry

The geometry of a doubly curved shell structure is defined by using the differential geometry principles. In the current section, only the fundamental aspects are shown. Nevertheless, the interested readers can find a complete treatise concerning the differential geometry in the book by Tornabene and Fantuzzi [2]. If a two-dimensional structural model is taken into account, any shell structure is represented by its middle surface, which is taken as the reference domain where the governing equations are written. The doubly curved surface in hand is identified by the position vector $r (α_{1}, α_{2})$ that allows to locate each point $P'$ of the reference surface in the global reference system $O x_{1} x_{2} x_{3}$ , as shown in Figure 1. It should be stated that $α_{1}, α_{2}$ denote the parametric lines of curvature of the middle surface, which are taken as orthogonal and principal coordinates in the present work. It is important to specify that coordinates $α_{1}, α_{2}$ are completely general and can have different meanings in dependence on the structural typology. For example, it is convenient to assign the following spherical coordinates $α_{1} = ϕ, α_{2} = ϑ$ to analyze a generic doubly curved shell of revolution. On the other hand, they become $α_{1} = x, α_{2} = y$ in the simpler case of a flat panel. It is worth noting that a shell structure is a three-dimensional solid in any case. Therefore, it occupies a well-defined region in the three-dimensional space. In connection with Figure 1, any point P within the shell is identified by the position vector $R (α_{1}, α_{2}, ζ)$ , defined as follows

R (α_{1}, α_{2}, ζ) = r (α_{1}, α_{2}) + \frac{h}{2} z n (α_{1}, α_{2})

(1)

where

z = 2 ζ / h (α_{1}, α_{2}) \in [- 1, 1]

is a non-dimensional parameter that stands for the distance of the point P from its projection

P'

on the shell middle surface, whereas

n (α_{1}, α_{2})

is the outward unit normal. Consequently, ζ represents the direction which is orthogonal to the middle surface along which the shell thickness h is measured. The three orthogonal coordinates

α_{1}, α_{2}, ζ

allow to define the local reference system

O' α_{1} α_{2} ζ

of origin

O'

(Figure 1). By setting the proper limits along these three principal directions, the three-dimensional solid is univocally defined. Such limitations can be written as

α_{1} \in [α_{1}^{0}, α_{1}^{1}], α_{2} \in [α_{2}^{0}, α_{2}^{1}], ζ \in [- h / 2, h / 2]

(2)

At this point, it is important to specify that the overall thickness of a generic laminated composite shell structure takes into account the superimposition of l laminae of thickness h_k, as it is shown in Figure 2. Thus, the following relation allows to define the total thickness of the structure

h = \sum_{k = 1}^{l} h_{k}

(3)

If a sandwich structures is considered, the shell is made of three layers only: two external face-sheets of thickness h_s and a central soft-core of thickness h_c. As a consequence, the overall thickness is given by the following expression

h = 2 h_{s} + h_{c}

(4)

assuming that the face-sheets have the same thickness. Once the position vector

r (α_{1}, α_{2})

of the middle surface is specified, the fundamental forms, and the corresponding fundamental coefficients, can be evaluated. They are needed to perform a complete characterization of the doubly curved surface at issue and play an important role in the governing equations. Following the same approach illustrated by the authors in their previous work [25], the differential geometry parameters can be presented by using a matrix notation. The matrix

I_{P}

is useful to collect the fundamental coefficients

E (α_{1}, α_{2}), F (α_{1}, α_{2}), G (α_{1}, α_{2})

of the first fundamental form

I_{P} = [E F F G] = [\frac{\partial r}{\partial α_{1}} \cdot \frac{\partial r}{\partial α_{1}} \frac{\partial r}{\partial α_{1}} \cdot \frac{\partial r}{\partial α_{2}} \frac{\partial r}{\partial α_{1}} \cdot \frac{\partial r}{\partial α_{2}} \frac{\partial r}{\partial α_{2}} \cdot \frac{\partial r}{\partial α_{2}}]

(5)

where the symbol “ċ” denotes the scalar product. On the other hand, the outward unit normal vector

n (α_{1}, α_{2})

is expressed as follows

n = \frac{\frac{\partial r}{\partial α_{1}} \times \frac{\partial r}{\partial α_{2}}}{| \frac{\partial r}{\partial α_{1}} \times \frac{\partial r}{\partial α_{2}} |}

(6)

where the symbol “×” stands for the vector product. The outward unit normal vector allows to define the second fundamental coefficients

L (α_{1}, α_{2}), M (α_{1}, α_{2}), N (α_{1}, α_{2})

. They are collected in the matrix

{II}_{P}

written below

{II}_{P} = [L M M N] = [- \frac{\partial^{2} r}{\partial α_{1}^{2}} \cdot n - \frac{\partial^{2} r}{\partial α_{1} \partial α_{2}} \cdot n - \frac{\partial^{2} r}{\partial α_{1} \partial α_{2}} \cdot n - \frac{\partial^{2} r}{\partial α_{2}^{2}} \cdot n]

(7)

By combining properly the two operators $I_{P}, {II}_{P}$ , a third matrix $S_{P}$ can be introduced. Recalling the linear algebra properties, it is defined as follows

S_{P} = I_{P}^{- 1} {II}_{P} = \frac{1}{EG - F^{2}} [LG - FM GM - FN EM - FL EN - FM]

(8)

Since the determinant of $I_{P}$ is never equal to zero, it can always be inverted. By means of the matrix $S_{P}$ , the Gaussian curvature $Γ_{G}$ and the mean curvature $Γ_{M}$ of the surface can be introduced as

Γ_{G} = det (S_{P}) = \frac{LN - M^{2}}{EG - F^{2}} = \frac{1}{R_{1} R_{2}}

(9)

Γ_{M} = \frac{1}{2} tr (S_{P}) = \frac{1}{2} \frac{LG - 2 FM + EN}{EG - F^{2}} = \frac{1}{2} (\frac{1}{R_{1}} + \frac{1}{R_{2}})

(10)

where

R_{1} = R_{1} (α_{1}, α_{2})

and

R_{2} = R_{2} (α_{1}, α_{2})

are the principal radii of curvature of the surface that can be written as

R_{1} = \frac{1}{Γ_{M} - \sqrt{Γ_{M}^{2} - Γ_{G}}} R_{2} = \frac{1}{Γ_{M} + \sqrt{Γ_{M}^{2} - Γ_{G}}}

(11)

According to the hypothesis of orthogonal and principal coordinates, the expressions of the principal radii of curvature can be simplified as

R_{1} = \frac{E}{L}, R_{2} = \frac{G}{N}

(12)

since the relation

F = M = 0

must be taken into account. At this point, the Lamé parameters

A_{1} = A_{1} (α_{1}, α_{2})

and

A_{2} = A_{2} (α_{1}, α_{2})

can be introduced using the first fundamental coefficients

A_{1} = \sqrt{E} A_{2} = \sqrt{G}

(13)

Finally, it is important to underline the fact that the current approach is suitable to analyze the mechanical behavior of thick and moderately thick shell structures, for which the following limitation is valid

0.01 \leq max (\frac{h}{R_{min}}, \frac{h}{L_{min}}) \leq 0.2

(14)

Shell fundamental equations

Even though shell structures are three-dimensional solids, the well-known three-dimensional theory of elasticity is too burdensome to be solved due to the huge number of degrees of freedom (dof) that are involved. Thus, a simplified two-dimensional model is introduced to analyze these structural problems. In this paper, an equivalent single layer (ESL) approach is presented. According to this theory, the original three-dimensional problem is converted into a two-dimensional problem due to proper hypotheses [12,14]. In the following, the three aspects of the elastic problem (equilibrium, kinematic and the constitutive laws) are presented.

The well-known Carrera Unified Formulation (CUF) [11] represents the framework which this work is based on. According to this formulation, several higher-order shear deformation theories (HSDTs) can be developed and employed to analyze the structural behavior of plates and shells in a unified manner. It is important to notice that a HSDT should be used to deal with particular mechanical configurations, such as laminated composite or sandwich structures. In this circumstance, the classic Reissner-Mindlin Theory or first-order shear deformation theory (FSDT) is inadequate and the solution does not represent correctly the real behavior of the considered structure. In general, the three-dimensional displacement component vector $U = U (α_{1}, α_{2}, ζ, t)$ can be written as

U = \sum_{τ = 0}^{N + 1} F_{τ} u^{(τ)}

(15)

in which

u^{(τ)} = u^{(τ)} (α_{1}, α_{2}, t)

is τ-th order generalized displacement component vector of the points located on the shell middle surface, whereas

F_{τ} = F_{τ} (ζ)

is the thickness function matrix. The index τ represents the order of kinematic expansion and it can be assumed in an arbitrary way. According to the value of τ, it is possible to define a different kinematic model. For the sake of clarity, the displacement field (15) becomes

U_{1} = F_{0} u_{1}^{(0)} + F_{1} u_{1}^{(1)} + F_{2} u_{1}^{(2)} + F_{3} u_{1}^{(3)} + \dots + F_{N} u_{1}^{(N)} + F_{N + 1} u_{1}^{(N + 1)} U_{2} = F_{0} u_{2}^{(0)} + F_{1} u_{2}^{(1)} + F_{2} u_{2}^{(2)} + F_{3} u_{2}^{(3)} + \dots + F_{N} u_{2}^{(N)} + F_{N + 1} u_{2}^{(N + 1)} U_{3} = F_{0} u_{3}^{(0)} + F_{1} u_{3}^{(1)} + F_{2} u_{3}^{(2)} + F_{3} u_{3}^{(3)} + \dots + F_{N} u_{3}^{(N)} + F_{N + 1} u_{3}^{(N + 1)}

(16)

It is important to underline that the choice of the thickness functions $F_{τ} = F_{τ} (ζ)$ is unrestricted, as it was proven by the huge number of functions of practical interest listed in the work [12]. The $(N + 1)$ -th order of kinematic expansion is linked to the Murakami’s function $Z = Z (ζ)$ , which can be defined as follows

Z = (- 1)^{k} (\frac{2}{ζ_{k + 1} - ζ_{k}} ζ - \frac{ζ_{k + 1} + ζ_{k}}{ζ_{k + 1} - ζ_{k}})

(17)

where

ζ_{k}

is the coordinate of the k-th layer along the normal direction ζ, as it can be observed in Figure 2. The Murakami’s function in hand allows to consider the zig-zag effect along the shell thickness. In other words, this function permits to take into account a continuous displacement field along the thickness that is characterized by a different slope between two adjacent layers [31 –33]. In general, the zig-zag effect is caused by a different stiffness of the constituent materials along the transverse direction. This is the case of sandwich structures with a soft-core highly deformable in the thickness direction. By selecting the well-known power function

ζ^{τ}

(with

τ = 0, 1, 2, \dots, N

) as thickness functions, the displacement field (16) takes the following aspect

U_{1} = u_{1}^{(0)} + ζ u_{1}^{(1)} + ζ^{2} u_{1}^{(2)} + ζ^{3} u_{1}^{(3)} + \dots + ζ^{N} u_{1}^{(N)} + Z (ζ) u_{1}^{(N + 1)} U_{2} = u_{2}^{(0)} + ζ u_{2}^{(1)} + ζ^{2} u_{2}^{(2)} + ζ^{3} u_{2}^{(3)} + \dots + ζ^{N} u_{2}^{(N)} + Z (ζ) u_{2}^{(N + 1)} U_{3} = u_{3}^{(0)} + ζ u_{3}^{(1)} + ζ^{2} u_{3}^{(2)} + ζ^{3} u_{3}^{(3)} + \dots + ζ^{N} u_{3}^{(N)} + Z (ζ) u_{3}^{(N + 1)}

(18)

According to the maximum order of kinematic expansion $τ = N$ , several theories can be defined. The following higher order structural models are obtained varying the order $τ = N$ up to four

where the letter “E” specifies that an ESL approach is followed, “D” indicates that the governing system of equations is written in terms of generalized displacements, whereas “Z” stands for the insertion of the Murakami’s function. It should be noticed that the well-known FSDT can be included in the present general approach choosing appropriately the degrees of freedom, even with the addition of the Murakami’s function (FSDTZ). The strain component vector

ε = ε (α_{1}, α_{2}, ζ, t)

defined as

ε = [ɛ_{1} ɛ_{2} γ_{12} γ_{13} γ_{23} ɛ_{3}]^{T}

(20)

can be written as a function of the τ-th order generalized displacement component vector

u^{(τ)}

, in the hypothesis of a principal curvilinear coordinate reference system

ε = \sum_{τ = 0}^{N + 1} Z^{(τ)} D_{Ω} u^{(τ)}

(21)

in which

D_{Ω}

is the differential operator related the curvatures terms defined as

D_{Ω} = [\frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} - \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} - \frac{1}{R_{1}} 0100 \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} \frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} - \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} 0 - \frac{1}{R_{2}} 010 \frac{1}{R_{1}} \frac{1}{R_{2}} 00 \frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} 001]^{T}

(22)

where the quantities

H_{1} = 1 + ζ / R_{1}

and

H_{2} = 1 + ζ / R_{2}

have been introduced. On the other hand,

Z^{(τ)}

is the matrix that comprises the terms linked to the shell thickness defined as follows

Z^{(τ)} = [\frac{F_{τ}}{H_{1}} 000000000 \frac{F_{τ}}{H_{2}} 000000000 \frac{F_{τ}}{H_{1}} \frac{F_{τ}}{H_{2}} 000000000 \frac{F_{τ}}{H_{1}} 0 \frac{\partial F_{τ}}{\partial ζ} 0000000 \frac{F_{τ}}{H_{2}} 0 \frac{\partial F_{τ}}{\partial ζ} 000000000 \frac{\partial F_{τ}}{\partial ζ}]

(23)

The τ-th order vector $ε^{(τ)} = ε^{(τ)} (α_{1}, α_{2}, t)$ that collects the generalized strain components defined on the shell reference surface, for $τ = 0, 1, 2, \dots, N, N + 1$ , assumes the aspect shown below

ε^{(τ)} = [ɛ_{1}^{(τ)} ɛ_{2}^{(τ)} γ_{1}^{(τ)} γ_{2}^{(τ)} γ_{13}^{(τ)} γ_{23}^{(τ)} ω_{13}^{(τ)} ω_{23}^{(τ)} ɛ_{3}^{(τ)}]^{T}

(24)

Having in mind expression (21), τ-th order generalized strain component vector is defined eventually as follows, for $τ = 0, 1, 2, \dots, N, N + 1$

ε^{(τ)} = D_{Ω} u^{(τ)}

(25)

The constitutive relations correlate the strain components and the corresponding stress components, introducing also the mechanical properties of the constituent materials. If a laminated composite structure is considered, some hypotheses must be introduced. In particular, discontinuities and voids are not allowed, since each lamina is a continuous body. In addition, a linear-elastic material must be chosen for each layer. Finally, composite materials can be considered reasonably as homogeneous mediums from the macroscopic point of view.

In the present work, the face-sheets that compose sandwich structures are assumed to be fiber-reinforced layers. For the sake of completeness, it should be recalled that a fiber-reinforced lamina is made of many fibers embedded in a matrix material that can be a polymer for example. Some further hypotheses concerning fiber-reinforced layers must be added. It should be recalled that the matrix and the fibers have to be perfectly bonded, and the fibers are taken as continuous and parallel. In addition, both the fibers and the matrix must be considered as isotropic materials. Finally, it should be assumed that voids or micro cracks, as well as residual stress, are not allowed in the matrix. As a result, a fiber-reinforced layer can be considered as an orthotropic medium, oriented by a generic angle $θ^{(k)}$ with respect to the local reference system $O' α_{1} α_{2} ζ$ shown in Figure 1. In general, the fiber orientation of the k-th layer $θ^{(k)}$ is kept constant in the whole domain. Nevertheless, the advancements in manufacturing processes have allowed a curvilinear placement of the reinforcing fibers. As a consequence, the orientation $θ^{(k)}$ can vary continuously within the reference domain. In the present paper, a general mathematical function is introduced by the authors to define the fiber path orientation in the k-th fiber-reinforced layer

θ^{(k)} (α_{1}, α_{2}) = θ_{1}^{(k)} ψ_{1}^{(k)} (α_{1}) + θ_{2}^{(k)} ψ_{2}^{(k)} (α_{2}) + θ_{3}^{(k)} ψ_{3}^{(k)} (α_{1}) + + θ_{4}^{(k)} ψ_{4}^{(k)} (α_{2}) + θ_{5}^{(k)} ψ_{5}^{(k)} (α_{1}, α_{2}) + θ_{6}^{(k)} ψ_{6}^{(k)} (α_{1}, α_{2}) + {\bar{θ}}^{(k)}

(26)

where

θ_{i}^{(k)}

are the angular parameters that control the variation of the fiber orientation,

{\bar{θ}}^{(k)}

is a constant value, whereas

ψ_{1}^{(k)} (α_{1})

ψ_{2}^{(k)} (α_{2})

ψ_{3}^{(k)} (α_{1})

ψ_{4}^{(k)} (α_{2})

ψ_{5}^{(k)} (α_{1}, α_{2})

and

ψ_{6}^{(k)} (α_{1}, α_{2})

are the functions that describe the curvilinear paths of the fibers within the reference domain. In particular, the functions

ψ_{1}^{(k)} (α_{1})

and

ψ_{3}^{(k)} (α_{1})

denote a variation along the first curvilinear coordinate

α_{1}

. They are defined as follows

in which

n_{1}, n_{3} \in ℕ_{0}^{+}

p_{1}, p_{3} \in ℝ

and

α_{1 m}, α_{3 m} \in [0, 1]

. Analogously, the functions

ψ_{2}^{(k)} (α_{2})

and

ψ_{4}^{(k)} (α_{2})

specify a variation along the second curvilinear coordinate

α_{2}

and assume the following representation

where

n_{2}, n_{4} \in ℕ_{0}^{+}

p_{2}, p_{4} \in ℝ

and

α_{2 m}, α_{4 m} \in [0, 1]

. These functions are all written as a function of the dimensionless coordinates

{\bar{α}}_{i}, {\tilde{α}}_{i}

defined in the following manner, having in mind the limitations in equation (2), for

i = 1, 2

{\bar{α}}_{i} = \frac{α_{i} - α_{i}^{0}}{α_{i}^{1} - α_{i}^{0}}, {\tilde{α}}_{i} = \frac{α_{i}^{1} - α_{i}}{α_{i}^{1} - α_{i}^{0}}

(29)

where

α_{i}^{0}, α_{i}^{1}

, for

i = 1, 2

, represent the boundary coordinates of the reference domain shown in equation (2). It is worth noting that the first relation in equations (27) and (28) is a power-law function, the second one is a sine-wave variation, whereas the third one is an exponential distribution. The last two functions

ψ_{5}^{(k)} (α_{1}, α_{2})

ψ_{6}^{(k)} (α_{1}, α_{2})

represent two-dimensional functions. The following definition can be employed to define both of them, assuming

i = 5, 6

ψ_{i}^{(k)} (α_{1}, α_{2}) = e^{- Ψ_{i}}

(30)

According to the value of $Ψ_{i}$ , the Gaussian distribution ( $i = 5$ ) centered in the point $(α_{5 m}, α_{6 m})$ , with $α_{5 m} \in [α_{1}^{0}, α_{1}^{1}]$ and $α_{6 m} \in [α_{2}^{0}, α_{2}^{1}]$ , and an elliptic variation ( $i = 6$ ), centered in the point $(α_{7 m}, α_{8 m})$ , with $α_{7 m} \in [α_{1}^{0}, α_{1}^{1}]$ and $α_{8 m} \in [α_{2}^{0}, α_{2}^{1}]$ can be obtained. The following value is given to the exponent in equation (30) to define a Gaussian function

Ψ_{5} = \frac{1}{2 (1 - ρ_{12}^{2})} ((\frac{α_{1} - α_{5 m}}{Δ_{1}})^{2} + (\frac{α_{2} - α_{6 m}}{Δ_{2}})^{2} - 2 ρ_{12} \frac{α_{1} - α_{5 m}}{Δ_{1}} \frac{α_{2} - α_{6 m}}{Δ_{2}})

(31)

where

ρ_{12} \in [- 1, 1]

and

Δ_{1}, Δ_{2} > 0

. On the other hand, the following value is used to define an elliptic variation

Ψ_{6} = \frac{((α_{1} - α_{7 m}) \cos β + (α_{2} - α_{8 m}) \sin β)^{2}}{a^{2}} + \frac{(- (α_{1} - α_{7 m}) \sin β + (α_{2} - α_{8 m}) \cos β)^{2}}{b^{2}}

(32)

In particular, it should be specified that the ellipse is identified by the semi-axes a and b and is affected by a rotation β, measured counter-clockwise from the first axis $α_{1}$ . The variation functions just introduced can be combined together to assign the desired orientation to the reinforcing fibers. As specified in the works by Nik et al. [45] and by Akbarzadeh et al. [46], it is not possible to leave the formation of defects, such as gaps and overlaps, out of consideration. Thus, the hypothesis that such defects are limited during the production process must be added. Once these fundamental hypotheses are illustrated, the constitutive relations can be introduced. The stress component vector $σ^{(k)} = σ^{(k)} (α_{1}, α_{2}, ζ, t)$ defined as follows for the k-th layer

σ^{(k)} = [σ_{1}^{(k)} σ_{2}^{(k)} τ_{12}^{(k)} τ_{13}^{(k)} τ_{23}^{(k)} σ_{3}^{(k)}]^{T}

(33)

can be related to the corresponding strain component vector in (20). The three-dimensional constitutive equation for an orthotropic lamina is defined by

σ^{(k)} = {\bar{C}}^{(k)} ε^{(k)}

(34)

in which

{\bar{C}}^{(k)}

is the constitutive operator written in geometric reference system

O' α_{1} α_{2} ζ

. For an orthotropic lamina, it becomes

{\bar{C}}^{(k)} = [{\bar{E}}_{11}^{(k)} {\bar{E}}_{12}^{(k)} {\bar{E}}_{16}^{(k)} 00 {\bar{E}}_{13}^{(k)} {\bar{E}}_{12}^{(k)} {\bar{E}}_{22}^{(k)} {\bar{E}}_{26}^{(k)} 00 {\bar{E}}_{23}^{(k)} {\bar{E}}_{16}^{(k)} {\bar{E}}_{26}^{(k)} {\bar{E}}_{66}^{(k)} 00 {\bar{E}}_{36}^{(k)} 000 {\bar{E}}_{44}^{(k)} {\bar{E}}_{45}^{(k)} 0000 {\bar{E}}_{45}^{(k)} {\bar{E}}_{55}^{(k)} 0 {\bar{E}}_{13}^{(k)} {\bar{E}}_{23}^{(k)} {\bar{E}}_{36}^{(k)} 00 {\bar{E}}_{33}^{(k)}]

(35)

where

{\bar{E}}_{nm}^{(k)}

represents the elastic constants. As specified in the works [25,26], they can stand for the plane stress-reduced elastic coefficients

{\bar{E}}_{nm}^{(k)} = {\bar{Q}}_{nm}^{(k)}

or for the non-reduced ones

{\bar{E}}_{nm}^{(k)} = {\bar{C}}_{nm}^{(k)}

. This choice depends on the fact that the structural theory can include (

{\bar{E}}_{nm}^{(k)} = {\bar{C}}_{nm}^{(k)}

) or not (

{\bar{E}}_{nm}^{(k)} = {\bar{Q}}_{nm}^{(k)}

) the stretching effects. In this case, the subscript

RS

will be added to the nomenclature used to indicate a structural theory that neglect this effect, such as the FSDT. The internal actions can be now expressed in terms of stress resultants

S^{(τ)} {= S}^{(τ)} (α_{1}, α_{2}, t)

. For a generic order τ of the kinematic expansion, the vector in hand assumes the following aspect

S^{(τ)} = [N_{1}^{(τ)} N_{2}^{(τ)} N_{12}^{(τ)} N_{21}^{(τ)} T_{1}^{(τ)} T_{2}^{(τ)} P_{1}^{(τ)} P_{2}^{(τ)} S_{3}^{(τ)}]^{T}

(36)

In general, the τ-th order generalized stress resultant vector, for $τ = 0, 1, 2, \dots, N, N + 1$ , is defined as follows

S^{(τ)} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} (Z^{(τ)})^{T} σ^{(k)} H_{1} H_{2} d ζ = \sum_{η = 0}^{N + 1} A^{(τ η)} D_{Ω} u^{(η)}

(37)

in which

A^{(τ η)}

, for

τ, η = 0, 1, 2, \dots, N, N + 1

, is the stiffness matrix. For an orthotropic medium, one gets

A^{(τ η)} = [A_{11 (20)}^{(τ η)} A_{12 (11)}^{(τ η)} A_{16 (20)}^{(τ η)} A_{16 (11)}^{(τ η)} 0000 A_{13 (10)}^{(τ \tilde{η})} A_{12 (11)}^{(τ η)} A_{22 (02)}^{(τ η)} A_{26 (11)}^{(τ η)} A_{26 (02)}^{(τ η)} 0000 A_{23 (01)}^{(τ \tilde{η})} A_{16 (20)}^{(τ η)} A_{26 (11)}^{(τ η)} A_{66 (20)}^{(τ η)} A_{66 (11)}^{(τ η)} 0000 A_{36 (10)}^{(τ \tilde{η})} A_{16 (11)}^{(τ η)} A_{26 (02)}^{(τ η)} A_{66 (11)}^{(τ η)} A_{66 (02)}^{(τ η)} 0000 A_{36 (01)}^{(τ \tilde{η})} 0000 A_{44 (20)}^{(τ η)} A_{45 (11)}^{(τ η)} A_{44 (10)}^{(τ \tilde{η})} A_{45 (10)}^{(τ \tilde{η})} 00000 A_{45 (11)}^{(τ η)} A_{55 (02)}^{(τ η)} A_{45 (01)}^{(τ \tilde{η})} A_{55 (01)}^{(τ \tilde{η})} 00000 A_{44 (10)}^{(\tilde{τ} η)} A_{45 (01)}^{(\tilde{τ} η)} A_{44 (00)}^{(\tilde{τ} \tilde{η})} A_{45 (00)}^{(\tilde{τ} \tilde{η})} 00000 A_{45 (10)}^{(\tilde{τ} η)} A_{55 (01)}^{(\tilde{τ} η)} A_{45 (00)}^{(\tilde{τ} \tilde{η})} A_{55 (00)}^{(\tilde{τ} \tilde{η})} 0 A_{13 (10)}^{(\tilde{τ} η)} A_{23 (01)}^{(\tilde{τ} η)} A_{36 (10)}^{(\tilde{τ} η)} A_{36 (01)}^{(\tilde{τ} η)} 0000 A_{33 (00)}^{(\tilde{τ} \tilde{η})}]

(38)

in which each element is defined as follows

A_{nm (pq)}^{(τ η)} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} {\bar{B}}_{nm}^{(k)} F_{η} F_{τ} \frac{H_{1} H_{2}}{H_{1}^{p} H_{2}^{q}} d ζ A_{nm (pq)}^{(\tilde{τ} η)} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} {\bar{B}}_{nm}^{(k)} F_{η} \frac{\partial F_{τ}}{\partial ζ} \frac{H_{1} H_{2}}{H_{1}^{p} H_{2}^{q}} d ζ A_{nm (pq)}^{(τ \tilde{η})} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} {\bar{B}}_{nm}^{(k)} \frac{\partial F_{η}}{\partial ζ} F_{τ} \frac{H_{1} H_{2}}{H_{1}^{p} H_{2}^{q}} d ζ A_{nm (pq)}^{(\tilde{τ} \tilde{η})} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} {\bar{B}}_{nm}^{(k)} \frac{\partial F_{η}}{\partial ζ} \frac{\partial F_{τ}}{\partial ζ} \frac{H_{1} H_{2}}{H_{1}^{p} H_{2}^{q}} d ζ

(39)

for

τ, η = 0, 1, 2, \dots, N, N + 1

n, m = 1, 2, 3, 4, 5, 6

and

p, q = 0, 1, 2

. It should be specified that the indices

τ, η

denote the corresponding thickness functions

F_{τ}, F_{η}

, whereas

\tilde{τ}, \tilde{η}

indicate the derivatives of the same thickness functions with respect to the coordinate ζ. The quantities

{\bar{B}}_{nm}^{(k)}

have the following meaning

{\bar{B}}_{nm}^{(k)} = {\bar{E}}_{nm}^{(k)} for n, m = 1, 2, 3, 6 {\bar{B}}_{nm}^{(k)} = κ {\bar{E}}_{nm}^{(k)} for n, m = 4, 5

(40)

in which

κ = 1 / χ

is the shear correction factor that cannot be neglected in those structural models that allow a non-parabolic shear stress profile along the thickness. In general, it is assumed equal to the constant value of

5 / 6

. For the sake of clarity, its value will be specified as superscript in the nomenclature used to identify those theories that could need it. As a consequence, the apex

χ = 1.2

specifies that the constant value of

5 / 6

is taken into account, whereas

χ = 1

denotes that the shear correction factor is neglected. According to the authors’ knowledge, the use of the shear correction factor is not always clear when a HSDT is considered. This is the case of the

{FSDTZ}_{RS}

and

EDZ 2

models. As explained in their previous work [25], the authors suggest that this factor should be added only if the laminae that compose the structure are all defined by the same mechanical properties. As a consequence, when a sandwich structure with a soft-core is analyzed, the shear correction factor will be set equal to

χ = 1

for the structural models in hand. In addition, as illustrated in the work [25], a numerical technique has to be employed to compute numerically the integrals (39).

Finally, the well-known Hamilton’s principle allows to deduce the equations of motion and the corresponding boundary conditions. For each order $τ = 0, 1, 2, \dots, N, N + 1$ of kinematic expansion, the three equations of motion can be written in matrix form as follows

D_{Ω}^{*} S^{(τ)} = \sum_{η = 0}^{N + 1} M^{(τ η)} {\overset{··}{u}}^{(η)}

(41)

in which

D_{Ω}^{*}

represents the equilibrium differential operator that assumes the following aspect

D_{Ω}^{*} = [\frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} - \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} - \frac{1}{R_{1}} - \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} - \frac{1}{R_{2}} \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} \frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} 0 \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} 0 \frac{1}{R_{1}} 0 \frac{1}{A_{1}} \frac{\partial}{\partial α_{1}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{2}}{\partial α_{1}} 0 \frac{1}{R_{2}} \frac{1}{A_{2}} \frac{\partial}{\partial α_{2}} + \frac{1}{A_{1} A_{2}} \frac{\partial A_{1}}{\partial α_{2}} - 1 000 - 1 000 - 1]^{T}

(42)

On the other hand, $M^{(τ η)}$ is the mass matrix defined as follows, for $τ, η = 0, 1, 2, \dots, N, N + 1$

M^{(τ η)} = [I^{(τ η)} 000 I^{(τ η)} 000 I^{(τ η)}]

(43)

where the inertia terms

I^{(τ η)}

are evaluated for

τ, η = 0, 1, 2, \dots, N, N + 1

once the mass density

ρ^{(k)}

of the k-th layer is introduced

I^{(τ η)} = \sum_{k = 1}^{l} \int_{ζ_{k}}^{ζ_{k + 1}} ρ^{(k)} F_{τ} F_{η} H_{1} H_{2} d ζ

(44)

Finally, ${\overset{··}{u}}^{(η)}$ is the vector that collects the generalized acceleration components that correspond to the generalized displacements. All the results just illustrated can be summarized in the so-called fundamental system of equations, containing $3 \times (N + 2)$ equilibrium equations. The compact form of this system of equations in written below for each order $τ = 0, 1, 2, \dots, N, N + 1$ of kinematic expansion

\sum_{η = 0}^{N + 1} L^{(τ η)} u^{(η)} = \sum_{η = 0}^{N + 1} M^{(τ η)} {\overset{··}{u}}^{(η)}

(45)

in which the fundamental operator

L^{(τ η)} = D_{Ω}^{*} A^{(τ η)} D_{Ω}

is a

3 \times 3

matrix defined as follows

L^{(τ η)} = [L_{11}^{(τ η)} L_{12}^{(τ η)} L_{13}^{(τ η)} L_{21}^{(τ η)} L_{22}^{(τ η)} L_{23}^{(τ η)} L_{31}^{(τ η)} L_{32}^{(τ η)} L_{33}^{(τ η)}]

(46)

The meaning of each term $L_{fg}^{(τ η)}$ , for $f, g = 1, 2, 3$ is illustrated in the work [14]. In order to complete the elastic problem in hand, the proper boundary conditions must be enforced for each order of kinematic expansion $τ = 0, 1, 2, \dots, N, N + 1$ . In the following, three different kinds of restraints are considered and mixed together. In particular, a generic edge can be completely clamped (C), simply supported (S) or free (F). For an edge identified by $α_{1} = α_{1}^{0}$ or $α_{1} = α_{1}^{1}$ , for $α_{2}^{0} \leq α_{2} \leq α_{2}^{1}$ , one gets

C \to u_{1}^{(τ)} = u_{2}^{(τ)} = u_{3}^{(τ)} = 0 S \to N_{1}^{(τ)} = 0, u_{2}^{(τ)} = u_{3}^{(τ)} = 0 F \to N_{1}^{(τ)} = N_{12}^{(τ)} = T_{1}^{(τ)} = 0

(47)

On the other hand, for an edge characterized by $α_{2} = α_{2}^{0}$ or $α_{2} = α_{2}^{1}$ , for $α_{1}^{0} \leq α_{1} \leq α_{1}^{1}$ , it is possible to write

C \to u_{1}^{(τ)} = u_{2}^{(τ)} = u_{3}^{(τ)} = 0 S \to N_{2}^{(τ)} = 0, u_{1}^{(τ)} = u_{3}^{(τ)} = 0 F \to N_{21}^{(τ)} = N_{2}^{(τ)} = T_{2}^{(τ)} = 0

(48)

In the particular case of structures with a closing meridian, such as complete shells of revolution or toroids, the structural compatibility must be enforced along the two coincident sides. If a shell is closed along $α_{1} = α_{1}^{0}$ and $α_{1} = α_{1}^{1}$ , the following relations have to be imposed for $τ = 0, 1, 2, \dots, N, N + 1$ and $α_{2}^{0} \leq α_{2} \leq α_{2}^{1}$

N_{1}^{(τ)} (α_{1}^{0}, α_{2}, t) = N_{1}^{(τ)} (α_{1}^{1}, α_{2}, t), u_{1}^{(τ)} (α_{1}^{0}, α_{2}, t) = u_{1}^{(τ)} (α_{1}^{1}, α_{2}, t) N_{12}^{(τ)} (α_{1}^{0}, α_{2}, t) = N_{12}^{(τ)} (α_{1}^{1}, α_{2}, t), u_{2}^{(τ)} (α_{1}^{0}, α_{2}, t) = u_{2}^{(τ)} (α_{1}^{1}, α_{2}, t) T_{1}^{(τ)} (α_{1}^{0}, α_{2}, t) = T_{1}^{(τ)} (α_{1}^{1}, α_{2}, t), u_{3}^{(τ)} (α_{1}^{0}, α_{2}, t) = u_{3}^{(τ)} (α_{1}^{1}, α_{2}, t)

(49)

Analogously, the conditions needed to impose the closure along $α_{2} = α_{2}^{0}$ and $α_{2} = α_{2}^{1}$ are written below

N_{21}^{(τ)} (α_{1}, α_{2}^{0}, t) = N_{21}^{(τ)} (α_{1}, α_{2}^{1}, t), u_{1}^{(τ)} (α_{1}, α_{2}^{0}, t) = u_{1}^{(τ)} (α_{1}, α_{2}^{1}, t) N_{2}^{(τ)} (α_{1}, α_{2}^{0}, t) = N_{2}^{(τ)} (α_{1}, α_{2}^{1}, t), u_{2}^{(τ)} (α_{1}, α_{2}^{0}, t) = u_{2}^{(τ)} (α_{1}, α_{2}^{1}, t) T_{2}^{(τ)} (α_{1}, α_{2}^{0}, t) = T_{2}^{(τ)} (α_{1}, α_{2}^{1}, t), u_{3}^{(τ)} (α_{1}, α_{2}^{0}, t) = u_{3}^{(τ)} (α_{1}, α_{2}^{1}, t)

(50)

With the aim to facilitate the identification of a specific edge, in the following, the sequence WSEN is employed. In other word, the boundary conditions follow the order just mentioned. With reference to the doubly curved panel of Figure 1, each edge is denoted by the coordinates written below

Westedge (W) Southedge (S) Eastedge (E) Northedge (N) \to \to \to \to α_{1}^{0} \leq α_{1} \leq α_{1}^{1}, α_{2} = α_{2}^{0} α_{1} = α_{1}^{1}, α_{2}^{0} \leq α_{2} \leq α_{2}^{1} α_{1}^{0} \leq α_{1} \leq α_{1}^{1}, α_{2} = α_{2}^{1} α_{1} = α_{1}^{0}, α_{2}^{0} \leq α_{2} \leq α_{2}^{1}

(51)

For the sake of clarity, some examples are now mentioned. In particular, the sequence CCCF means that the panel is completely clamped along the west (W), south (S), and east (E) edges, whereas it is free on the North (N) side. For a closed structure involved by the conditions (49) or (50), the edges on which the structural compatibility is enforced do not appear in the complete sequence WSEN.

Numerical procedure

The well-known GDQ method is used to solve numerically the fundamental system of equations (45) in terms of generalized displacements. For conciseness purposes, in the following only the main ideas of this approach are presented. Thus, the reader can consult the review paper by Tornabene et al. [69] for a more general and complete view on this technique. With reference to a one-dimensional domain, the GDQ method approximates the n-th derivative at a generic point x_i of a sufficiently smooth function $f (x)$ as a weighted linear sum of the function values at some chosen nodes

\frac{d^{n} f (x)}{{dx}^{n}} |_{x = x_{i}} ≅ \sum_{j = 1}^{T} ς_{ij}^{(n)} f (x_{j})

(52)

for

i = 1, 2, \dots, T

, in which T is the total number of grid points, whereas

ς_{ij}^{(n)}

are the weighting coefficients of the sum that can be evaluated by using the recursive formulae given by Shu [68], also illustrated in Tornabene et al. [69]. Taking into account a two-dimensional domain, such as the one used as a reference for a shell structure, the total number of grid points is set equal to

T = I_{N}

along

α_{1}

and equal to

T = I_{M}

along

α_{2}

. Due to the proven features of stability and accuracy, the discrete points are located within the reference domain according to the Chebyshev-Gauss-Lobatto grid distribution. Mathematically speaking, it is defined as follows

α_{ki} = (1 - \cos (\frac{i - 1}{T - 1} π)) \frac{(α_{k}^{1} - α_{k}^{0})}{2} + α_{k}^{0}

(53)

for

i = 1, 2, \dots, T

and

α_{k} \in [α_{k}^{0}, α_{k}^{1}]

. The index

k = 1, 2

is used to specify if the grid is applied along

α_{1}

direction (

k = 1

) or along the second one (

k = 2

). It should be specified that the same ideas of the GDQ method are employed to solve numerically the integrals of the elastic coefficients

A_{nm (pq)}^{(τ η)}

A_{nm (pq)}^{(\tilde{τ} η)}

A_{nm (pq)}^{(τ \tilde{η})}

A_{nm (pq)}^{(\tilde{τ} \tilde{η})}

, according to the rules of the so-called generalized integral quadrature (GIQ) technique. The reader can refer to the review paper by Tornabene et al. [69] for further details about this approach.

By means of the GDQ method, the fundamental system of equations (45) can be written in its discrete form, which takes the following aspect

K δ = ω^{2} M δ

(54)

in which K, M and

δ

denote, respectively, the global stiffness matrix, the mass matrix and the mode shape vector. On the other hand, ω is the circular frequency that allows to evaluate the corresponding natural frequency as

f = ω / 2 π

. It is well known that this last equation is a generalized linear eigenvalue. Thus, the free vibration analysis aims to solve this relation.

Results and discussion

The numerical tests presented in the current section aim to investigate the dependency of the natural frequencies of several sandwich structures with variable stiffness on the curvilinear path of the reinforcing fibers. Before showing some results related to doubly curved shell structures, the free vibration analysis of some variable stiffness plates is performed in order to prove the validity of the present approach. In particular, a preliminary comparison study is accomplished taking into account the results available in the literature. Then, the convergence and stability features of the numerical technique in hand are proven by the convergence analysis. Afterwards, some new results are shown. For convenience purposes, the mechanical properties of each material used in the following computations are listed in Table 1. All the solutions are obtained by using a MATLAB code [70] that implements the formulation illustrated above.

Table 1.

Mechanical properties of the materials.

Material	Mechanical properties
Ceramic foam	$E = 0.232 GPa$ , $ν = 0.2$ , $ρ = 320 kg / m^{3}$
Aluminum honey-comb	$E_{1} = E_{2} = 4.4 \cdot 10^{- 4} GPa$ , $E_{3} = 1.536 GPa$ , $G_{12} = 2.7 \cdot 10^{- 4} GPa$ , $G_{13} = 0.369 GPa$ , $G_{23} = 0.222 GPa$ , $ν_{12} = 0.99$ , $ν_{13} = ν_{23} = 8.64 \cdot 10^{- 5}$ , $ρ = 72 kg / m^{3}$
CFRP [67]	$E_{1} = 173 GPa$ , $E_{2} = E_{3} = 7.2 GPa$ , $G_{12} = G_{13} = G_{23} = 3.76 GPa$ , $ν_{12} = ν_{13} = ν_{23} = 0.29$ , $ρ = 1540 kg / m^{3}$
Glass-epoxy	$E_{1} = 53.78 GPa$ , $E_{2} = E_{3} = 17.93 GPa$ , $G_{12} = G_{13} = 8.96 GPa$ , $G_{23} = 3.45 GPa$ , $ν_{12} = ν_{13} = 0.25$ , $ν_{23} = 0.34$ , $ρ = 1900 kg / m^{3}$
Graphite-epoxy	$E_{1} = 137.9 GPa$ , $E_{2} = E_{3} = 8.96 GPa$ , $G_{12} = G_{13} = 7.1 GPa$ , $G_{23} = 6.21 GPa$ , $ν_{12} = ν_{13} = 0.3$ , $ν_{23} = 0.49$ , $ρ = 1450 kg / m^{3}$

Comparison with literature

The flat sandwich panel analyzed by Vescovini and Dozio [67] in their recent work is taken into account in this paragraph with the aim to compare the natural frequencies of a variable stiffness sandwich plate obtained through a variable-kinematic formulation based on the well-known Ritz method and the GDQ method. The FFFC plate under consideration is made of seven plies. The aluminum honey-comb is considered as core material of thickness h_c, whereas each face-sheet of thickness h_s is made of CFRP (Table 1). These orthotropic layers are affected by a variable distribution of the mechanical properties due to the curvilinear placement of the reinforcing fibers. The curvilinear paths are defined through the classic two-parameter law which is controlled by the fiber orientation in both the center of the domain and along the edges. The same law is employed by the authors in the works [57,66]. In particular, in the present case, the lay-up is set symmetrically with respect to the core. The stacking sequence of the three external laminae is given by the expression

(á 0 | 45 ñ_{x} / á - 45 | - 60 ñ_{x} / á 0 | 45 ñ_{x})

, in which the notation

á ϑ_{0} | ϑ_{1} ñ_{x}

denotes a variable stiffness ply. In particular,

ϑ_{0}

is the orthotropic angle in the center of the plate, whereas

ϑ_{1}

represents the fiber orientation along the external edges. The subscript x means that such variation is applied along x direction. The reader can find further details and comments about this approach in the works [57,66]. The size of the plate under consideration is specified by the parameters

L_{x}, L_{y}

, whereas

h = 6 h_{s} + h_{c}

is its total thickness. The ratio between the total thickness h and the thickness of each face-sheet h_s is assumed to be equal to

30

. The frequencies are expressed in their dimensionless form as

\bar{ω} = \frac{2 π f L_{x}^{2}}{h} \sqrt{\frac{ρ_{s}}{E_{2}}}

(55)

where

ρ_{s}

and E₂ are, respectively, the mass density and the transverse Young’s modulus of the CFRP. The results are shown in Table 2 for three different ratios of

L_{y} / h

and

L_{x} / L_{y}

. In general, all the GDQ solutions are in good agreement with the reference ones for all these cases. In particular, very similar values are achieved when the corresponding structural models are compared. In particular,

ED 2^{χ = 1.2}

EDZ 2^{χ = 1}

ED 4

and

EDZ 4

theories are used. It should be observed that the Murakami’s function does not introduce remarkable differences since the squeezing effect of the core is not so significant, as it can be noticed from the mechanical properties of the involved materials. A final observation should be addressed now. Since the results in Table 2 are shown for different ratios

L_{x} / L_{y}

, it is necessary to increase the number of discrete node to obtain an accurate solution by means of the GDQ method. In particular, when the plate is a square (

L_{x} / L_{y} = 1

), the same number of points is chosen along each principal direction (

I_{N} = I_{M} = 25

). On the other hand, for

L_{x} / L_{y} = 2

and

L_{x} / L_{y} = 3

, the GDQ solution is obtained by setting, respectively,

I_{N} = 41

and

I_{N} = 51

, whereas

I_{M} = 25

is unvaried.

Table 2.

Comparison of the non-dimensional frequencies $\bar{ω} = 2 π f L_{x}^{2} h^{- 1} \sqrt{ρ_{s} / E_{2}}$ for a FFFC sandwich rectangular plate with face-sheets of CFRP and core made of aluminum honey-comb, with $h / h_{s} = 30$ .

$\frac{L_{x}}{L_{y}}$	Structural model	Mode
		1	2	3	4	1	2	3	4	1	2	3	4
		$L_{y} / h = 10$				$L_{y} / h = 50$				$L_{y} / h = 100$
1	ED2 [67]	3.6821	7.9407	11.9758	18.1479	4.1002	9.5788	23.7127	31.8724	4.1429	9.7631	24.2744	32.5761
	ED2 GDQ	3.6189	7.7754	11.9616	17.6568	4.0937	9.5646	23.6054	31.6804	4.1387	9.7606	24.1524	32.4834
	EDZ2 GDQ	3.6539	7.8847	11.9592	17.9554	4.1102	9.5938	23.6748	31.7544	4.1519	9.7790	24.1962	32.4966
	ED4 [67]	3.3610	6.8627	11.9512	14.9936	4.0417	9.3462	22.9283	30.6815	4.1142	9.6468	23.9209	32.1586
	ED4 GDQ	3.3517	6.8732	11.9578	15.0061	4.0780	9.4006	22.9840	30.5763	4.1107	9.6641	24.0294	32.0689
	EDZ4 GDQ	3.3528	6.8731	11.9574	15.0031	4.0828	9.4088	23.0161	30.5721	4.1231	9.6996	24.1102	32.0673
	FEM [67]	3.3477	6.9041	11.5623	15.0779	4.0120	9.3128	22.8421	30.4028	4.0924	9.6315	23.9126	32.0361
2	ED2 [67]	3.5643	14.1935	15.4039	26.0942	3.8265	16.1481	29.8060	52.9823	3.8657	16.3793	30.0566	53.6760
	ED2 GDQ	3.5294	14.0157	15.3280	25.5615	3.8611	16.1149	29.6427	52.8270	3.6947	16.3974	30.2013	53.6024
	EDZ2 GDQ	3.5358	14.0883	15.3185	25.8592	3.8801	16.1471	29.6619	52.8990	3.7046	16.4262	30.2199	53.5738
	ED4 [67]	3.3970	12.8589	15.3777	23.0854	3.7824	15.8836	29.4262	52.0042	3.8409	16.2393	29.9095	53.2596
	ED4 GDQ	3.3749	12.8451	15.2255	23.0498	3.8991	16.0148	29.3414	52.2439	3.8886	16.4255	29.9390	53.8902
	EDZ4 GDQ	3.3763	12.8458	15.2252	23.0426	3.9778	16.1115	29.6305	52.7201	3.9532	16.5052	29.9015	53.8744
	FEM [67]	3.3728	12.9171	14.8626	23.1550	3.7612	15.8460	29.3271	51.9736	3.8250	16.2171	29.8703	53.3410
3	ED2 [67]	3.3612	16.2699	19.9127	27.9935	3.6067	22.1080	30.3972	62.2630	3.6465	22.3713	30.5915	62.8547
	ED2 GDQ	3.3247	16.2220	19.6639	27.6554	3.5341	22.0175	30.4769	62.3013	3.6422	22.3545	30.4898	62.1591
	EDZ2 GDQ	3.3435	16.2191	19.8102	27.8504	3.5454	22.0484	30.4945	62.3751	3.6547	22.3823	30.4708	62.0768
	ED4 [67]	3.2297	16.2496	18.3947	25.9012	3.5645	21.8159	30.1325	61.5151	3.6229	22.2130	30.4748	62.5038
	ED4 GDQ	3.2112	16.1498	18.3990	25.8712	3.5444	21.8239	30.2171	61.7529	3.5984	22.5301	30.8204	63.5057
	EDZ4 GDQ	3.2124	16.1495	18.4017	25.8678	3.5522	21.8790	30.2532	61.9104	3.6892	23.1492	31.1542	64.7257
	FEM [67]	3.1980	15.6472	18.4272	25.9416	3.5373	21.7200	30.1303	61.4477	3.6000	22.1270	30.5229	62.5005

Note: The plate is made of seven plies, and the stacking sequence is given by $(á 0 | 45 ñ_{x} / á - 45 | - 60 ñ_{x} / á 0 | 45 ñ_{x} / Core / á 0 | 45 ñ_{x} / á - 45 | - 60 ñ_{x} / á 0 | 45 ñ_{x})$ . All the face-sheets have the same thickness h_s.

Convergence analysis

In this paragraph, the CCCC sandwich square plate of side $L_{x} = L_{y} = 1 m$ depicted in Figure 3 is taken into account in order to show the convergence and stability characteristics of the GDQ method. The structure is made of three layers and the stacking sequence is given by $(θ_{1}^{(1)} ψ_{1}^{(1)} (x) / Soft - Core / 90 + θ_{1}^{(3)} ψ_{1}^{(3)} (x))$ . The soft-core is made of ceramic foam ( $h_{c} = 0.03 m$ ) and the graphite-epoxy is used as face-sheet material ( $h_{s} = 0.01 m$ ). Thus, the total thickness of the plate is $h = 0.05 m$ . Having in mind the general expression (26) for the curvilinear fiber path definition, a sine-wave variation is considered in each face-sheet. In particular, the functions $ψ_{1}^{(1)} (x) = ψ_{1}^{(3)} (x) = (\sin (π (n_{1} \bar{x} + x_{m})))^{p_{1}}$ , with $p_{1} = 1$ , $n_{1} = 1$ , $x_{m} = 0.5$ and $θ_{1}^{(1)} = θ_{1}^{(3)} = π / 3$ , are employed to describe the variable mechanical properties of the external layers. For the sake of clarity, the graphical representation of the curvilinear placement of the reinforcing phase is depicted in the lower part of Figure 3. Finally, the first 10 natural frequencies of the square plate at issue are shown in Table 3 for different structural theories, increasing the number of grid points $I_{N} = I_{M}$ along the two principal directions. A good convergence of the results can be observed as few as 21 points for all the considered structural models.

Table 3.

Convergence characteristics of the first 10 natural frequencies for a CCCC sandwich square plate increasing the number of points $I_{N} = I_{M}$ for the Chebyshev-Gauss-Lobatto grid distribution.

f (Hz)	$I_{N} = I_{M} = 11$	$I_{N} = I_{M} = 15$	$I_{N} = I_{M} = 21$	$I_{N} = I_{M} = 25$	$I_{N} = I_{M} = 31$	$I_{N} = I_{M} = 41$
	${FSDTZ}_{RS}^{χ = 1}$
1	230.539	243.727	243.790	243.793	243.794	243.794
2	401.756	401.154	400.958	400.941	400.935	400.934
3	407.872	414.855	414.894	414.896	414.897	414.898
4	539.959	539.200	539.086	539.079	539.078	539.077
5	634.560	625.817	626.134	626.150	626.155	626.156
6	672.100	636.312	636.331	636.331	636.332	636.332
7	672.100	734.740	734.665	734.660	734.659	734.659
8	735.539	748.182	748.521	748.528	748.530	748.530
9	889.280	875.318	874.978	874.969	874.965	874.965
10	896.570	891.445	891.489	891.489	891.489	891.490
	${EDZ 2}^{χ = 1}$
1	238.601	244.183	244.283	244.284	244.286	244.288
2	403.298	402.330	402.098	402.077	402.072	402.073
3	407.156	415.826	415.870	415.871	415.874	415.875
4	541.858	540.781	540.648	540.639	540.638	540.640
5	636.223	628.074	628.134	628.152	628.160	628.162
6	677.736	638.063	638.070	638.070	638.072	638.073
7	677.736	736.946	736.862	736.856	736.857	736.858
8	738.076	750.268	750.448	750.456	750.462	750.463
9	889.973	877.863	877.769	877.760	877.760	877.761
10	899.157	893.902	893.933	893.934	893.937	893.938
	$EDZ 3$
1	238.507	244.131	244.226	244.227	244.228	244.229
2	403.215	402.247	402.011	401.989	401.984	401.984
3	407.076	415.731	415.768	415.768	415.770	415.771
4	541.743	540.664	540.524	540.513	540.512	540.513
5	636.080	627.928	627.986	628.002	628.010	628.012
6	677.501	637.911	637.906	637.904	637.905	637.906
7	677.501	736.775	736.680	736.672	736.672	736.673
8	737.911	750.077	750.248	750.254	750.258	750.260
9	889.730	877.659	877.551	877.540	877.539	877.539
10	898.937	893.678	893.692	893.690	893.691	893.691
	$EDZ 4$
1	237.672	243.584	243.621	243.609	243.604	243.605
2	401.952	400.870	400.522	400.476	400.457	400.458
3	406.076	414.285	414.197	414.171	414.162	414.163
4	539.677	538.330	538.015	537.970	537.951	537.952
5	633.228	624.769	624.646	624.626	624.616	624.618
6	674.464	634.681	634.486	634.448	634.433	634.435
7	674.464	732.504	732.172	732.117	732.095	732.096
8	734.087	745.806	745.726	745.686	745.669	745.671
9	883.689	871.777	871.224	871.154	871.126	871.127
10	892.779	887.576	887.296	887.236	887.213	887.215

Note: The face-sheets are made of graphite-epoxy, whereas the soft-core is made of ceramic foam.

Free vibration analysis of variable stiffness singly and doubly curved shells

The free vibration analysis of four sandwich shell structures with variable stiffness is the main topic of the current section. Particular attention is given to the effect of the orthotropic angle variation on the dynamic response. Thus, the following results are organized as parametric studies in which the natural frequencies are obtained as a function of the angle parameter φ, whose meaning will be specified time by time. The structures in hand are depicted in Figure 4, as well as the global and local coordinate reference systems. It should be noticed that the local reference system $O' α_{1} α_{2} ζ$ has a great importance for the curvilinear path description according to the general law (26). In fact, the direction affected by the variable orientation of the fibers follows the principal parametric lines $α_{1}, α_{2}$ . Once the outward normal to the surface is identified, the other two directions can be immediately recognized since this reference system is always right-handed. The considered shells are the following: a truncated cone (Figure 4(a)), a doubly curved revolution panel with catenary meridian (Figure 4(b)), a doubly curved translational panel (Figure 4(c)) and an ellipsoidal panel (Figure 4(d)). Following the framework given by the differential geometry, the position vectors $r (α_{1}, α_{2})$ and all the geometric parameters required to describe their reference surfaces are listed in Table 4. In addition, also the lamination schemes of each sandwich structure are shown in the same table, together with the functions and the variables needed to characterize the curvilinear path of the reinforcing fibers. In all the following numerical computations, the ceramic foam of Table 1 is assumed as core material. Finally, it should be highlighted that the case of straight fibers is always contemplated in the parametric studies. Thus, the comparison with the corresponding three-dimensional model obtained by a finite element model (3D-FEM) is performed in this circumstance. For this purpose, the commercial finite element code Abaqus and Straus are employed. As far as the GDQ solution is concerned, four different ESL models which include the Murakami’s function are taken into account in each analysis: ${FSDTZ}_{RS}^{χ = 1}$ , ${EDZ 2}^{χ = 1}$ , $EDZ 3$ and $EDZ 4$ . Finally, the number of grid points employed to discretize the domain is specified in each case.

Table 4.

Definition of the position vectors $r (α_{1}, α_{2})$ and lamination schemes for the sandwich shell structures depicted in Figure 3.

(a) Conical shell (singly curved panel of revolution)

Reference surface

r (α_{1}, y) = (a \cos α_{1} + y \sin α \sin ϕ) e_{1} - y \cos α e_{2} + (b \sin α_{1} - y \sin α \cos ϕ) e_{3} ϕ (α_{1}) = \arctan (- \frac{b \cos α_{1}}{a \sqrt{1 - \cos^{2} α_{1}}}), a = b = 0.5 m, α = π / 12, α_{1} \in [α_{1}^{0}, α_{1}^{1}] = [0, 2 π], y \in [y_{0}, y_{1}] = [0, 1 m]

Lamination scheme

(90 + θ_{2}^{(1)} ψ_{2}^{(1)} (y) / Soft - Core / 90 + θ_{2}^{(3)} ψ_{2}^{(3)} (y)), θ_{2}^{(1)} = - θ_{2}^{(3)} = φ, h_{s} = 0.02 m, h_{c} = 0.06 m ψ_{2}^{(1)} (y) = {\bar{α}}_{2}^{p_{2}}, with p_{2} = 1 ψ_{2}^{(3)} (y) = {\bar{α}}_{2}^{p_{2}}, with p_{2} = 1

(b) Doubly curved panel of revolution with catenary meridian

Reference surface

r (ϕ, ϑ) = R_{0} (ϕ) \cos ϑ e_{1} - R_{0} (ϕ) \sin ϑ e_{2} + a arcsinh (\frac{1}{\tan ϕ}) e_{3} R_{0} (ϕ) = a \cosh (arcsinh (\frac{1}{\tan ϕ})), a = 2 m, c = 2 m, d = 1 m α_{1} \in [α_{1}^{0}, α_{1}^{1}] = [- 2.05117, 2.05117], α_{2} \in [α_{2}^{0}, α_{2}^{1}] = [- π / 6, π / 6]

Lamination scheme

(30 + θ_{1}^{(1)} ψ_{1}^{(1)} (ϕ) / Soft - Core / 45 + θ_{5}^{(3)} ψ_{5}^{(3)} (ϕ, ϑ)), θ_{1}^{(1)} = - 2 θ_{5}^{(3)} = φ, h_{s} = 0.05 m, h_{c} = 0.1 m ψ_{1}^{(1)} (ϕ) = e^{(- \frac{\bar{ϕ}}{ϕ_{m}})^{p_{1}}}, with p_{1} = 1, ϕ_{m} = 0.9 ψ_{5}^{(3)} (ϕ, ϑ) = e^{(- \frac{1}{2 (1 - ρ_{12}^{2})} ((\frac{ϕ - ϕ_{m}}{Δ_{1}})^{2} + (\frac{ϑ - ϑ_{m}}{Δ_{2}})^{2} - 2 ρ_{12} \frac{ϕ - ϕ_{m}}{Δ_{1}} \frac{ϑ - ϑ_{m}}{Δ_{2}}))}, with ϕ_{m} = 0, ϑ_{m} = \frac{ϑ_{1} - ϑ_{0}}{2}, Δ_{1} = 0.3, Δ_{2} = 100, ρ_{12} = 0

Reference surface

r (α_{1}, α_{2}) = (R_{0}^{α_{1}} (α_{1}) - x_{3}^{α_{2}} (α_{2}) \sin α_{1}) e_{1} + - R_{0}^{α_{2}} (α_{2}) e_{2} + (x_{3}^{α_{1}} (α_{1}) + x_{3}^{α_{2}} (α_{2}) \cos α_{1}) e_{3} R_{0}^{α_{1}} (α_{1}) = \frac{k^{α_{1}} \tan α_{1}}{2}, x_{3}^{α_{1}} (α_{1}) = \frac{(R_{0}^{α_{1}} (α_{1}))^{2}}{k^{α_{1}}}, k^{α_{1}} = \frac{a^{2} - d^{2}}{b}, a = 2 m, b = 0.5 m, c = - 2 m, d = 0 m R_{0}^{α_{2}} (α_{2}) = R \sin α_{2}, x_{3}^{α_{2}} (α_{2}) = R (1 - \cos α_{2}), R = 3 m, α_{1} \in [α_{1}^{0}, α_{1}^{1}] = [- \frac{π}{6}, \frac{π}{6}], α_{2} \in [α_{2}^{0}, α_{2}^{1}] = [- π / 4, π / 4]

Lamination scheme

(30 + θ_{3}^{(1)} ψ_{3}^{(1)} (α_{1}) / Soft - Core / 60 + θ_{1}^{(3)} ψ_{1}^{(3)} (α_{1})), θ_{3}^{(1)} = - θ_{1}^{(3)} = φ, h_{s} = 0.02 m, h_{c} = 0.06 m ψ_{3}^{(1)} (α_{1}) = {\tilde{α}}_{1}^{p_{3}}, with p_{3} = 1 ψ_{1}^{(3)} (α_{1}) = {\bar{α}}_{1}^{p_{1}}, with p_{1} = 1

(d) Ellipsoid (doubly curved panel)

Reference surface

r (α_{1}, α_{2}) = (\frac{a \sqrt{a^{2} - b^{2} \sin^{2} α_{1} - c^{2} \cos^{2} α_{1}}}{\sqrt{a^{2} - c^{2}}} \cos α_{2}) e_{1} + (b \cos α_{1} \sin α_{2}) e_{2} + (- \frac{c \sqrt{a^{2} \sin^{2} α_{2} + b^{2} \cos^{2} α_{2} - c^{2}}}{\sqrt{a^{2} - c^{2}}} \sin α_{1}) e_{3} a = 2 m, b = 1.5 m, c = 1 m α_{1} \in [α_{1}^{0}, α_{1}^{1}] = [π / 6, 5 π / 6], α_{2} \in [α_{2}^{0}, α_{2}^{1}] = [π / 4, 3 π / 4]

Lamination scheme

(θ_{2}^{(1)} ψ_{2}^{(1)} (α_{2}) + θ_{3}^{(1)} ψ_{3}^{(1)} (α_{1}) / Soft - Core / 90 + θ_{1}^{(3)} ψ_{1}^{(3)} (α_{1})), θ_{3}^{(1)} = - θ_{2}^{(1)} = θ_{1}^{(3)} = φ, h_{c} = (0.2 - 2 h_{s}) m ψ_{2}^{(1)} (α_{2}) = {\bar{α}}_{2}^{p_{2}}, with p_{2} = 1, ψ_{3}^{(1)} (α_{1}) = {\tilde{α}}_{1}^{p_{3}}, with p_{3} = 0 ψ_{1}^{(3)} (α_{1}) = (\sin (π (n_{1} {\bar{α}}_{1} + α_{1 m})))^{p_{1}}, with p_{1} = 1, n_{1} = 0.5, α_{1 m} = 0

The first application is related to a CC sandwich conical shell (Figure 4(a)). The external face-sheets are made of glass-epoxy and the curvilinear fiber paths are described by a linear variation in each ply. The first 10 natural frequencies are shown in Table 5 as a function of the angle parameter

θ_{2}^{(1)} = - θ_{2}^{(3)} = φ

included in the interval

[0^{\circ}, 60^{\circ}]

. When

φ = 0^{\circ}

is considered, the frequencies are in good agreement with the 3D-FEM solution carried out by Abaqus (48640 C3D20 elements), especially if the EDZ4 theory is considered. It can be noticed that all the structural models give comparable values. The variation of the first eight natural frequencies can be easily observed from the graphs depicted in Figure 5. It should be noticed that the peculiar placement of the curvilinear fibers implicates the mode shape symmetry; therefore, the natural frequencies are equal two by two.

Table 5.

Frequency variations of a CC sandwich conical shell for different structural models as a function of the parameter φ of the general variation.

f (Hz)	3D-FEM Abaqus C3D20 (650256dof) $(φ = 0^{\circ})$	$φ = 0^{\circ}$	$φ = 15^{\circ}$	$φ = 30^{\circ}$	$φ = 45^{\circ}$	$φ = 60^{\circ}$
		${FSDTZ}_{RS}^{χ = 1}$
1	390.024	393.752	398.487	409.460	419.528	422.350
2	390.636	393.752	398.487	409.460	419.528	422.350
3	418.566	424.422	427.617	434.848	441.035	442.400
4	418.935	424.422	427.617	434.848	441.035	442.400
5	426.413	428.327	435.163	451.484	468.274	476.703
6	427.215	428.327	435.163	451.484	468.274	476.703
7	484.628	492.385	494.892	500.622	505.909	508.615
8	484.805	492.385	494.892	500.622	505.909	508.615
9	552.187	553.216	560.854	579.652	595.018	600.772
10	553.022	553.216	560.854	579.652	595.018	600.772
		${EDZ 2}^{χ = 1}$
1	390.024	392.579	397.278	408.116	417.867	420.254
2	390.636	392.579	397.278	408.116	417.867	420.254
3	418.566	422.594	425.726	432.778	438.655	439.600
4	418.935	422.594	425.726	432.778	438.655	439.600
5	426.413	428.248	435.056	451.231	467.580	475.281
6	427.215	428.248	435.056	451.231	467.580	475.281
7	484.628	491.056	493.455	498.907	503.781	505.919
8	484.805	491.056	493.455	498.907	503.781	505.919
9	552.187	553.779	561.404	580.078	594.194	599.146
10	553.022	553.779	561.404	580.078	594.194	599.146
		$EDZ 3$
1	390.024	392.502	397.204	408.047	417.804	420.190
2	390.636	392.502	397.204	408.047	417.804	420.190
3	418.566	422.441	425.577	432.637	438.522	439.469
4	418.935	422.441	425.577	432.637	438.522	439.469
5	426.413	428.220	435.028	451.205	467.555	475.253
6	427.215	428.220	435.028	451.205	467.555	475.253
7	484.628	490.824	493.228	498.691	503.576	505.717
8	484.805	490.824	493.228	498.691	503.576	505.717
9	552.187	553.764	561.389	580.062	593.926	598.885
10	553.022	553.764	561.389	580.062	593.926	598.885
		$EDZ 4$
1	390.024	391.149	395.857	406.704	416.429	418.742
2	390.636	391.149	395.857	406.704	416.429	418.742
3	418.566	419.757	422.899	429.962	435.816	436.694
4	418.935	419.757	422.899	429.962	435.816	436.694
5	426.413	427.543	434.357	450.539	466.874	474.521
6	427.215	427.543	434.357	450.539	466.874	474.521
7	484.628	486.043	488.455	493.927	498.797	500.905
8	484.805	486.043	488.455	493.927	498.797	500.905
9	552.187	553.363	560.989	579.669	586.357	591.351
10	553.022	553.363	560.989	579.669	586.357	591.351

Note: The Chebyshev-Gauss-Lobatto grid distribution is used with $I_{N} = 31$ and $I_{M} = 25$ . The face-sheets are made of glass-epoxy, whereas the soft-core is made of ceramic foam.

The second example aims to evaluate the frequency variation of the CCCF doubly curved panel of revolution with catenary meridian depicted in Figure 4(b) as a function of the angle parameter

θ_{1}^{(1)} = - 2 θ_{5}^{(3)} = φ

, taken in the interval

[0^{\circ}, 60^{\circ}]

. The work by Viola et al. [12] can be taken as reference to understand the meaning of the geometric parameters listed in Table 4 for the description of this particular structure. In this circumstance, an exponential distribution along the first principal coordinate ϕ of the shell middle surface is employed to describe the curvilinear fiber in the lower face-sheet, whereas a Gaussian variation which involves both the principal coordinates

ϕ, ϑ

is added to the fiber orientation of the upper layer. Both face-sheets are made of graphite-epoxy. The numerical values are shown in Table 6, in which a good agreement can be seen by means of the comparison with the 3D-FEM solution of Straus (10800 Hexa20 elements) for

φ = 0^{\circ}

, whereas the corresponding graphical representation of the first four frequencies is depicted in Figure 6. Even in this case, comparable results are achieved with regard to each structural model.

Table 6.

Frequency variations of a CCCF doubly curved panel of revolution with a catenary meridian for different structural models as a function of the parameter φ of the general variation.

f (Hz)	3D-FEM Straus Hexa20 (146547dof) $(φ = 0^{\circ})$	$φ = 0^{\circ}$	$φ = 15^{\circ}$	$φ = 30^{\circ}$	$φ = 45^{\circ}$	$φ = 60^{\circ}$
		${FSDTZ}_{RS}^{χ = 1}$
1	203.641	203.933	214.971	220.665	224.639	234.173
2	227.702	228.922	231.931	249.367	264.541	272.253
3	235.320	236.232	239.867	252.437	269.537	279.488
4	272.263	272.855	268.703	268.631	289.605	309.519
5	286.970	287.399	279.761	283.683	309.726	325.450
6	300.445	302.111	302.873	316.101	326.743	357.024
7	325.185	327.306	320.641	318.974	348.583	364.206
8	338.968	341.850	342.314	345.318	353.096	382.252
9	352.388	355.480	357.885	363.867	375.112	383.320
10	363.738	366.505	371.649	386.041	400.892	417.377
		${EDZ 2}^{χ = 1}$
1	203.641	203.646	215.070	221.092	225.174	234.753
2	227.702	229.102	232.183	249.951	265.476	273.572
3	235.320	236.393	240.334	252.585	270.865	280.151
4	272.263	273.061	269.293	269.802	289.878	309.966
5	286.970	287.951	280.515	284.304	310.087	327.475
6	300.445	302.350	303.479	316.780	327.839	357.437
7	325.185	328.128	322.029	320.676	348.709	363.492
8	338.968	343.483	344.829	348.360	355.828	380.814
9	352.388	357.798	360.983	366.917	378.171	386.224
10	363.738	368.434	374.476	388.288	402.485	420.480
		$EDZ 3$
1	203.641	203.614	215.039	221.063	225.144	234.720
2	227.702	229.034	232.148	249.891	265.405	273.518
3	235.320	236.335	240.271	252.545	270.820	280.075
4	272.263	273.018	269.251	269.754	289.825	309.907
5	286.970	287.906	280.480	284.269	310.041	327.425
6	300.445	302.280	303.418	316.725	327.781	357.355
7	325.185	328.060	321.970	320.626	348.665	363.420
8	338.968	343.396	344.759	348.293	355.762	380.759
9	352.388	357.655	360.897	366.839	378.086	386.128
10	363.738	368.328	374.362	388.191	402.395	420.386
		$EDZ 4$
1	203.641	203.308	214.683	220.456	224.441	233.903
2	227.702	227.616	231.427	248.797	264.091	271.888
3	235.320	235.306	239.052	252.010	269.196	278.636
4	272.263	271.967	268.066	268.145	288.655	308.273
5	286.970	286.386	279.177	283.245	308.715	324.768
6	300.445	300.717	301.854	315.244	325.494	354.805
7	325.185	325.404	319.384	317.978	346.923	360.778
8	338.968	339.393	340.642	344.119	351.998	379.062
9	352.388	352.497	356.321	362.609	373.968	382.026
10	363.738	364.115	369.752	384.591	398.751	414.812

Note: The Chebyshev-Gauss-Lobatto grid distribution is used with $I_{N} = 31$ and $I_{M} = 25$ . The face-sheets are made of graphite-epoxy, whereas the soft-core is made of ceramic foam.

The next application is focused on the free vibration analysis of the CCCC doubly curved sandwich panel of translation [12] depicted in Figure 4(c). This panel is obtained by sliding a circumference over a parabola. The reader can find the meaning of each geometric parameter shown in Table 4, which are needed to describe the parabolic profile, in the previous paper by the authors [57]. The variable mechanical properties are now obtained by applying a linear variation along the first principal coordinate

α_{1}

to the fiber orientation in both the two face-sheets made of graphite-epoxy, as it can be noticed from the lamination scheme specified in Table 4. The results in term of natural frequencies are shown in Table 7 as a function of the parameter

θ_{3}^{(1)} = - θ_{1}^{(3)} = φ

, with

φ \in [0^{\circ}, 60^{\circ}]

, for several higher order structural models. When

φ = 0^{\circ}

is chosen, an excellent agreement is observed with the corresponding 3D-FEM solution given by Abaqus (59280 C3D20 elements). For the sake of clarity, the same results of Table 7 are shown in their graphic form in Figure 7, where the variation of the first four frequencies is depicted for the various structural models. It is clear that each of them gives comparable results, since each curve does not detach considerably from the others.

Table 7.

Frequency variations of a CCCC doubly curved sandwich panel of translation (a circumference slides over a parabola) for different structural models as a function of the parameter φ of the general variation.

f (Hz)	3D-FEM Abaqus C3D20 (779559dof) $(φ = 0^{\circ})$	$φ = 0^{\circ}$	$φ = 15^{\circ}$	$φ = 30^{\circ}$	$φ = 45^{\circ}$	$φ = 60^{\circ}$
		${FSDTZ}_{RS}^{χ = 1}$
1	118.467	118.492	122.108	126.803	131.857	135.158
2	128.039	128.053	131.693	135.116	141.075	146.128
3	153.573	153.679	156.002	157.082	155.967	154.080
4	158.498	158.566	166.776	163.997	158.738	155.907
5	163.200	163.241	171.937	173.895	175.008	176.712
6	170.531	170.605	175.784	183.556	185.957	186.455
7	182.425	182.557	183.206	186.791	191.806	189.504
8	189.186	189.322	188.034	193.180	201.024	200.876
9	197.504	197.605	204.742	210.109	206.533	202.548
10	197.709	197.806	208.997	212.044	213.874	216.603
		${EDZ 2}^{χ = 1}$
1	118.467	118.534	122.176	126.869	131.888	135.177
2	128.039	128.139	131.758	135.169	141.099	146.143
3	153.573	153.688	156.003	157.061	155.942	154.064
4	158.498	158.601	166.773	163.977	158.726	155.907
5	163.200	163.385	171.997	174.005	175.117	176.816
6	170.531	170.756	175.855	183.550	185.984	186.463
7	182.425	182.690	183.321	186.869	191.725	189.461
8	189.186	189.480	188.154	193.234	201.054	201.058
9	197.504	197.809	204.867	210.267	206.634	202.395
10	197.709	198.005	209.137	212.159	214.078	216.800
		$EDZ 3$
1	118.467	118.523	122.165	126.858	131.877	135.166
2	128.039	128.125	131.745	135.157	141.088	146.131
3	153.573	153.675	155.991	157.052	155.933	154.057
4	158.498	158.591	166.762	163.964	158.713	155.896
5	163.200	163.368	171.983	173.984	175.098	176.798
6	170.531	170.737	175.839	183.537	185.967	186.449
7	182.425	182.671	183.301	186.850	191.713	189.448
8	189.186	189.454	188.128	193.216	201.035	201.032
9	197.504	197.785	204.846	210.241	206.611	202.380
10	197.709	197.983	209.115	212.139	214.052	216.773
		$EDZ 4$
1	118.467	118.463	122.091	126.775	131.791	135.082
2	128.039	128.037	131.650	135.064	140.993	146.021
3	153.573	153.583	155.901	156.976	155.849	153.984
4	158.498	158.505	166.630	163.816	158.580	155.784
5	163.200	163.202	171.848	173.747	174.885	176.601
6	170.531	170.547	175.670	183.429	185.815	186.314
7	182.425	182.447	183.081	186.668	191.610	189.322
8	189.186	189.205	187.893	193.020	200.795	200.687
9	197.504	197.528	204.617	209.887	206.297	202.223
10	197.709	197.726	208.843	211.887	213.697	216.380

Note: The Chebyshev-Gauss-Lobatto grid distribution is used with $I_{N} = I_{M} = 31$ . The face-sheets are made of graphite-epoxy, whereas the soft-core is made of ceramic foam.

For completeness purposes, the results presented up to this moment are summarized in Figure 8, in which the relative variation $f_{i} / f_{0} - 1$ of the first four natural frequencies of the structures just analyzed is depicted for the EDZ4 theory. It is possible to notice that a relative variation of $0.15 \div 0.20$ can be reached specially for the first frequency by varying the fiber orientation of the face-sheets, according to what has been explained before. In addition, it is evident that a greater variation is obtained for higher values of φ.

The CFCC ellipsoidal sandwich panel of Figure 4(d) is the shell structure analyzed in the last numerical application. Even in this case, the face-sheets are made of graphite-epoxy, whereas the soft-core is made of ceramic foam. Further details about this peculiar geometry can be found in the work by Tornabene [23]. In contrast with the previous examples, only one structural theory is taken into account (EDZ4). Thus, the parametric study is carried out by varying the parameter

θ_{3}^{(1)} = - θ_{2}^{(1)} = θ_{1}^{(3)} = φ

in the interval

[0^{\circ}, 60^{\circ}]

, for different values of the face-sheet thickness h_s, keeping constant the total thickness

h = h_{tot} = 0.2 m

. Two fiber variations affect simultaneously the lower layer. In particular, a linear variation is applied along both the two principal curvilinear coordinates

α_{1}, α_{2}

. On the other hand, a sine-wave function is employed to describe the curvilinear path of the fiber along

α_{1}

in the upper ply. Table 8 shows the frequency variations of this structure for several values of face-sheet thicknesses as a function of the parameter φ of the general variation. A good agreement can be observed by comparing the 3D-FEM solution given by Straus (13872 Hexa20 elements) with the GDQ results for

φ = 0^{\circ}

. The variation of the first four frequencies is depicted in Figure 9 for various ratios

h_{s} / h_{tot}

. For the sake of completeness, the relative variation

f_{i} / f_{0} - 1

of the first four natural frequencies of the ellipsoid is shown in Figure 10. It is clear that the relative variation becomes more relevant increasing the face-sheet thickness. In this case, the maximum relative variation that can be reached is about

0.30 \div 0.35

for the first frequency.

Table 8.

Frequency variations of a CFCC ellipsoid for several values of face-sheet thicknesses as a function of the parameter φ of the general variation.

f (Hz)	3D-FEM Straus Hexa20 (184695dof) $(φ = 0^{\circ})$	$φ = 0^{\circ}$	$φ = 15^{\circ}$	$φ = 30^{\circ}$	$φ = 45^{\circ}$	$φ = 60^{\circ}$
		$h_{s} = 0.01 m$
1	129.118	129.169	122.702	112.390	105.097	101.020
2	171.667	170.707	171.239	170.997	168.772	165.818
3	213.090	213.568	209.373	206.758	201.809	193.261
4	258.428	258.789	257.351	256.587	254.763	249.451
5	266.460	265.572	264.053	261.289	260.379	260.183
6	278.315	277.186	277.306	275.762	272.947	272.333
7	308.100	307.653	306.221	304.194	302.743	301.723
8	309.273	308.959	309.653	314.098	319.722	323.713
9	359.023	356.475	354.974	353.169	350.709	346.919
10	365.419	364.700	363.558	361.118	358.848	355.597
		$h_{s} = 0.02 m$
1	151.437	151.657	142.268	127.722	117.909	112.685
2	180.688	180.843	181.788	180.670	177.139	173.436
3	244.885	245.275	240.231	234.869	226.262	215.238
4	287.073	287.253	282.141	276.033	269.253	261.618
5	289.527	289.894	285.630	282.786	280.975	277.996
6	295.040	295.607	299.096	299.370	299.667	301.894
7	330.011	330.391	329.331	326.039	324.243	322.411
8	333.847	334.546	333.356	338.356	344.005	348.148
9	367.914	368.469	366.515	363.140	357.583	350.851
10	382.787	383.365	381.777	378.948	376.695	374.221
		$h_{s} = 0.04 m$
1	172.570	172.950	160.103	140.955	128.264	121.751
2	196.090	196.695	197.205	192.981	186.407	181.089
3	275.161	275.622	270.401	263.763	252.412	238.884
4	321.543	321.988	312.358	298.367	286.286	276.371
5	326.646	327.249	321.246	317.480	314.149	308.602
6	330.033	331.170	335.606	336.880	338.723	340.973
7	370.847	371.612	372.065	367.662	364.273	361.067
8	380.217	381.181	377.153	381.111	379.918	372.068
9	425.627	427.531	418.410	401.376	392.075	392.474
10	434.959	436.314	434.153	433.840	432.732	428.088
		$h_{s} = 0.06 m$
1	188.324	188.727	174.349	151.853	136.721	128.991
2	223.004	223.968	221.896	213.586	203.737	196.737
3	301.310	301.928	295.520	289.795	278.382	262.965
4	361.558	362.353	359.384	337.929	320.835	308.451
5	368.395	369.587	360.081	357.673	354.926	348.870
6	385.231	387.377	384.492	382.858	382.106	381.265
7	427.374	428.955	429.732	425.261	419.610	414.971
8	439.923	441.928	436.902	438.836	434.936	424.149
9	506.183	508.586	503.269	469.869	452.032	450.300
10	525.370	529.564	513.047	516.682	516.975	505.843

Note: The Chebyshev-Gauss-Lobatto grid distribution is used with $I_{N} = I_{M} = 25$ . The face-sheets are made of graphite-epoxy, whereas the soft-core is made of ceramic foam. The GDQ solutions are obtained by using the EDZ4 theory.

Finally, the first three mode shapes for all the sandwich structures considered in the present study are depicted in Figure 11. These representations are shown taking into account the EDZ4 theory and the maximum value of the fiber orientation parameter φ. For the ellipsoidal panel, the ratio $h_{s} / h_{tot} = 0.3$ is employed. For the sake of completeness, the curvilinear paths of the fibers and the corresponding orientations for the face-sheets are depicted in Figure 12 in the dimensionless domain, taking into account the maximum value of φ.

Conclusions

A new general expression which is able to mix together different functions, such as power-law, sine-wave and exponential distributions, has been introduced in the present paper to define the curvilinear paths of the reinforcing fibers of the corresponding orthotropic layers. The effect of this variable orientation on the modal response of several composite sandwich plates and shells has been investigated by means of some parametric studies. It has been proven that a considerable variation in terms of natural frequencies can be obtained by varying the fiber orientation and increasing the thickness of the external face-sheets. The present results are helpful to design and manufacture doubly curved shells and plates with variable stiffness. Several higher order structural models have been employed combined with the Murakami’s function in order to catch the zig-zag effect which characterizes most of the common lamination schemes of sandwich structures. The solution has been obtained by the GDQ method, whose accuracy has been proven comparing the numerical values with the results available in the literature or given by some FE commercial codes. As a direct future development, the same general approach could be employed to investigate the effect of the curvilinear reinforcing fibers on the static behavior of shells and plates with variable stiffness. In particular, a parametric study could be performed to show how the stress and strain profiles along the thickness are affected by the fiber orientation.

Footnotes

Acknowledgements

The research topic is one of the subjects of the Centre of Study and Research for the Identification of Materials and Structures (CIMEST)-“M. Capurso” of the University of Bologna (Italy).

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Reddy

. Mechanics of laminated composite plates and shells, Boca Raton, FL, USA: CRC Press, 2004.

Tornabene

Fantuzzi

. Mechanics of laminated composite doubly-curved shell structures. The generalized differential quadrature method and the strong formulation finite Element method, Bologna, Italy: Esculapio, 2014.

Groh

RMJ

Weaver

. A computationally efficient 2D model for inherently equilibrated 3D stress predictions in heterogeneous laminated plates. Part I: model formulation. Compos Struct 2016; 156: 171–185.

Groh

RMJ

Weaver

. A computationally efficient 2D model for inherently equilibrated 3D stress predictions in heterogeneous laminated plates. Part II: model validation. Compos Struct 2015; 156: 186–217.

Groh

RMJ

Weaver

. Static inconsistencies in certain axiomatic higher-order shear deformation theories for beams, plates and shells. Compos Struct 2015; 120: 231–245.

Groh

RMJ

Weaver

. On displacement-based and mixed-variational equivalent single layer theories for modelling highly heterogeneous laminated beams. Int J Solids Struct 2015; 59: 147–170.

Camanho

Davila

De Moura

. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003; 37: 1415–1438.

Borrelli

Riccio

Sellitto

et al.

On the use of global-local kinematic coupling approaches for delamination growth simulation in stiffened composite panels. Compos Sci Technol 2015; 115: 43–51.

Borrelli

Caputo

Riccio

et al.

Kinematic approach for a global-local coupling: compressive behaviour of a delaminated panel. Adv Compos Lett 2011; 20: 150–156.

10.

Riccio

Damiano

Raimondo

et al.

A fast numerical procedure for the simulation of inter-laminar damage growth in stiffened composite panels. Compos Struct 2016; 145: 203–216.

11.

Carrera

. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Arch Comput Methods Eng 2002; 9: 87–140.

12.

Viola

Tornabene

Fantuzzi

. General higher-order shear deformation theories for the free vibration analysis of completely doubly-curved laminated shells and panels. Compos Struct 2013; 95: 639–666.

13.

Tornabene

Reddy

. FGM and laminated doubly-curved and degenerate shells resting on nonlinear elastic foundations: a GDQ solution for static analysis with a posteriori stress and strain recovery. J Indian Inst Sci 2013; 93: 635–688.

14.

Tornabene

Viola

Fantuzzi

. General higher-order equivalent single layer theory for free vibrations of doubly-curved laminated composite shells and panels. Compos Struct 2013; 104: 94–117.

15.

Tornabene

Fantuzzi

Bacciocchi

et al.

A new approach for treating concentrated loads in doubly-curved composite deep shells with variable radii of curvature. Compos Struct 2015; 131: 433–452.

16.

Tornabene

Fantuzzi

Bacciocchi

et al.

Accurate inter-laminar recovery for plates and doubly-curved shells with variable radii of curvature using layer-wise theories. Compos Struct 2015; 124: 368–393.

17.

Tornabene

Fantuzzi

Bacciocchi

et al.

Dynamic analysis of thick and thin elliptic shell structures made of laminated composite materials. Compos Struct 2015; 133: 278–299.

18.

Tornabene

Fantuzzi

Bacciocchi

et al.

Free vibrations of composite oval and elliptic cylinders by the generalized differential quadrature method. Thin-Wall Struct 2015; 97: 114–129.

19.

Tornabene

Brischetto

Fantuzzi

et al.

Numerical and exact models for free vibration analysis of cylindrical and spherical shell panels. Compos Part B-Eng 2015; 81: 231–250.

20.

Tornabene

Fantuzzi

Viola

et al.

Stress and strain recovery for functionally graded free-form and doubly-curved sandwich shells using higher-order equivalent single layer theory. Compos Struct 2015; 119: 67–89.

21.

Fantuzzi

Bacciocchi

Tornabene

et al.

Radial basis functions based on differential quadrature method for the free vibration of laminated composite arbitrary shaped plates. Compos Part B-Eng 2015; 78: 65–78.

22.

Tornabene

Fantuzzi

Viola

. Inter-laminar stress recovery procedure for doubly-curved, singly-curved, revolution shells with variable radii of curvature and plates using generalized higher-order theories and the local GDQ method. Mech Adv Mater Struct 2016; 23: 1019–1045.

23.

Tornabene

. General higher order layer-wise theory for free vibrations of doubly-curved laminated composite shells and panels. Mech Adv Mater Struct 2016; 23: 1046–1067.

24.

Tornabene

Fantuzzi

Bacciocchi

et al.

Effect of agglomeration on the natural frequencies of functionally graded carbon nanotube-reinforced laminated composite doubly-curved shells. Compos Part B-Eng 2016; 89: 187–218.

25.

Tornabene

Fantuzzi

Bacciocchi

. The local GDQ method for the natural frequencies of doubly-curved shells with variable thickness: a general formulation. Compos Part B Eng 2016; 92: 265–289.

26.

Bacciocchi

Eisenberger

Fantuzzi

et al.

Vibration analysis of variable thickness plates and shells by the generalized differential quadrature method. Compos Struct 2016; 156: 218–237.

27.

Carrera

. A refined multi-layered finite-element model applied to linear and non-linear analysis of sandwich plates. Compos Sci Technol 1998; 58: 1553–1569.

28.

Saha

Kalamkarov

. Micromechanical thermoelastic model for sandwich composite shells made of generally orthotropic materials. J Sandw Struct Mater 2009; 11: 27–56.

29.

Brischetto

Carrera

Demasi

. Improved response of unsymmetrically laminated sandwich plates by using zig-zag functions. J Sandw Struct Mater 2009; 11: 257–267.

30.

Brischetto

. Hygrothermoelastic analysis of multilayered composite and sandwich shells. J Sandw Struct Mater 2013; 15: 168–202.

31.

Maturi

Ferreira

AJM

Zenkour

et al.

Analysis of laminated shells by Murakami’s zig-zag theory and radial basis functions collocation. J Appl Math 2013; 2013: 1–14.

32.

Carrera

. Historical review of zig-zag theories for multilayered plates and shells. Appl Mech Rev 2003; 56: 287–308.

33.

Carrera

. On the use of the Murakami’s zig-zag function in the modeling of layered plates and shells. Comput Struct 2004; 82: 541–554.

34.

Ribeiro

Akhavan

Teter

et al.

A review on the mechanical behaviour of curvilinear fibre composite laminated panels. J Compos Mater 2014; 48: 2761–2777.

35.

Waldhart C. Design of tow-placed, variable-stiffness laminates. Master of Science in Engineering Mechanics, Virginia Polytechnic Institute and State University, USA, 1996.

36.

Blom AW. Structural performance of fiber-placed, variable-stiffness composite conical and cylindrical shells. PhD Thesis, Delft University of Technology, 2010.

37.

Sendeckyj

. Longitudinal shear modulus of filamentary composite containing curvilinear fibers. Fibre Sci Technol 1970; 2: 211–222.

38.

Hyer

Charette

. Use of curvilinear fiber format in composite structure design. AIAA J 1991; 29: 1011–1015.

39.

Parnas

Oral

Ceyhan

. Optimum design of composite structures with curved fiber courses. Compos Sci Tech 2003; 63: 1071–1082.

40.

Blom

Setoodeh

Hol

JMAM

et al.

Design of variable-stiffness conical shells for maximum fundamental eigenfrequency. Comput Struct 2008; 86: 870–870.

41.

Blom

Tatting

Hol

JMAM

et al.

Fiber path definitions for elastically tailored conical shells. Compos Part B Eng 2009; 40: 77–84.

42.

Blom

Abdalla

Gürdal

. Optimization of course locations in fiber-placed panels for general fiber angle distributions. Compos Sci Technol 2010; 70: 564–570.

43.

Rouhi

Ghayoor

Hoa

et al.

Stiffness tailoring of elliptical composite cylinders for axial buckling performance. Compos Struct 2016; 150: 115–123.

44.

Coburn

Weaver

. Buckling analysis, design and optimisation of variable-stiffness sandwich panels. Int J Solids Struct 2016; 96: 217–228.

45.

Nik

Fayazbakhsh

Pasini

et al.

Optimization of variable stiffness composites with embedded defects induced by Automated Fiber Placement. Compos Struct 2014; 107: 160–165.

46.

Akbarzadeh

Nik

Pasini

. The role of shear deformation in laminated plates with curvilinear fiber paths and embedded defects. Compos Struct 2015; 118: 217–227.

47.

Kim

Weaver

Potter

. Computer aided modelling of variable angle tow composites manufactured by continuous tow shearing. Compos Struct 2015; 129: 256–267.

48.

Groh

RMJ

Weaver

. Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects. Compos Struct 2014; 107: 0 482–493.

49.

Coburn

Weaver

. Buckling analysis of stiffened variable angle tow panels. Compos Struct 2014; 111: 259–270.

50.

Zucco

Groh

RMJ

Madeo

et al.

Mixed shell element for static and buckling analysis of variable angle tow composite plates. Compos Struct 2016; 152: 324–338.

51.

Raju

Weaver

. Postbuckling analysis of variable angle tow plates using differential quadrature method. Compos Struct 2013; 106: 74–84.

52.

White

Raju

Weaver

. Initial post-buckling of variable-stiffness curved panels. J Mech Phys Solids 2014; 71: 132–155.

53.

Raju

Weaver

. Buckling and postbuckling of variable angle tow composite plates under in-plane shear loading. Int J Solids Struct 2015; 58: 270–287.

54.

White

Weaver

. Post-buckling analyses of variable-stiffness composite cylinders in axial compression. Compos Struct 2015; 123: 190–203.

55.

Akhavan

Ribeiro

. Natural modes of vibration of variable stiffness composite laminates with curvilinear fibers. Compos Struct 2011; 93: 3040–3047.

56.

Yazdani

Ribeiro

. A layerwise p-version finite element formulation for free vibration analysis of thick composite laminates with curvilinear fibres. Compos Struct 2015; 120: 531–542.

57.

Tornabene

Fantuzzi

Bacciocchi

et al.

Higher-order theories for the free vibration of doubly-curved laminated panels with curvilinear reinforcing fibers by means of a local version of the GDQ method. Compos Part B Eng 2015; 81: 196–230.

58.

Honda

Narita

. Natural frequencies and natural modes of laminated composite plates reinforced with arbitrary curvilinear fiber shape paths. J Sound Vib 2012; 331: 180–191.

59.

Akhavan

Ribeiro

. Free geometrically nonlinear oscillations of perfect and imperfect laminates with curved fibres by the shooting method. Nonlinear Dynam 2015; 81: 949–965.

60.

Akhavan

Ribeiro

. Non-linear forced periodic oscillations of laminates with curved fibres by the shooting method. Int J Nonlinear Mech 2015; 76: 176–189.

61.

Ribeiro

. Non-linear modes of vibration of thin cylindrical shells in composite laminates with curvilinear fibres. Compos Struct 2015; 122: 184–197.

62.

Diaz

Fagiano

Abdalla

et al.

A study of interlaminar stresses in variable stiffness plates. Compos Struct 2012; 94: 1192–1199.

63.

Akhavan

Ribeiro

De Moura

MFSF

. Large deflection and stresses in variable stiffness composite laminates with curvilinear fibres. Int J Mech Sci 2013; 73: 14–26.

64.

Yazdani

Ribeiro

. Geometrically non-linear static analysis of unsymmetric composite plates with curvilinear fibres: p-version layerwise approach. Compos Struct 2014; 118: 74–85.

65.

Yazdani

Ribeiro

Rodrigues

. A p-version layerwise model for large deflection of composite plates with curvilinear fibres. Compos Struct 2014; 108: 181–190.

66.

Tornabene

Fantuzzi

Bacciocchi

. Higher-order structural theories for the static analysis of doubly-curved laminated composite panels reinforced by curvilinear fibers. Thin-Wall Struct 2016; 102: 222–245.

67.

Vescovini

Dozio

. A variable-kinematic model for variable stiffness plates: vibration and buckling analysis. Compos Struct 2016; 142: 15–26.

68.

Shu

. Differential quadrature and its application in engineering, Berlin, Germany: Springer, 2000.

69.

Tornabene

Fantuzzi

Ubertini

et al.

Strong formulation finite element method based on differential quadrature: a survey. Appl Mech Rev 2015; 67: 020801–020801.

70.

Viola E, Tornabene F, Fantuzzi N, et al. DiQuMASPAB Software, DICAM Department, Alma Mater Studiorum – University of Bologna, http://software.dicam.unibo.it/diqumaspab-project (accessed 4 June 2017).