Composite Structures: Francesco Tornabene, Matteo Viscoti, Rossana Dimitri, Luciano Rosati
Composite Structures: Francesco Tornabene, Matteo Viscoti, Rossana Dimitri, Luciano Rosati
                                                                           Composite Structures
                                                            journal homepage: www.elsevier.com/locate/compstruct
A R T I C L E I N F O A B S T R A C T
Keywords:                                                       In the present contribution a general formulation is proposed to account for general boundary conditions within
Anisotropic laminated structures                                the dynamic analysis of anisotropic laminated doubly-curved shell having arbitrary shape and variable thickness.
Doubly-curved shells                                            Different analytical expressions are considered for the shell thickness variation along the geometrical principal
Equivalent Single Layer
                                                                directions, and the distortion of the physical domain is described by a mapping procedure based on Non-Uniform
GDQ method
General boundary conditions
                                                                Rational Basis Spline (NURBS) curves. Mode frequencies and shapes are determined employing higher-order
                                                                theories within an Equivalent Single Layer (ESL) framework. The related fundamental relations are tackled
                                                                numerically by means of the Generalized Differential Quadrature (GDQ) method. The dynamic problem is
                                                                derived from the Hamiltonian Principle, leading to a strong formulation of the governing equations. General
                                                                external constraints are enforced along the edges of the shell employing a distribution of linear springs
                                                                distributed on the faces of the three-dimensional solid and accounting for a spatial coordinate-dependent stiffness
                                                                along both in-plane and out-of-plane directions. Moreover, a Winkler-type foundation with general distribution
                                                                of linear springs is modelled on the top and bottom surfaces of the shell. A systematic set of numerical examples
                                                                is carried out for the validation of the proposed theory by comparing mode frequencies with predictions from
                                                                refined three-dimensional finite element analyses. Finally, we perform a sensitivity analysis of the dynamic
                                                                response of mapped curved structures for different spring stiffnesses and general external constraints, according
                                                                to various kinematic assumptions.
1. Introduction                                                                                        specifically, moderately thick shells can be analysed according to the so-
                                                                                                       called Equivalent Single Layer (ESL) approach, whereas a Layer-Wise
    Recent advances in many engineering fields encourage the employ                                   (LW) formulation is required by thick shells [5–6]. Unlike LW, in the
ment of curved lightweight structures in the design process. Most ap                                  case of ESL theories, the number of field variables does not depend on
plications require an optimized static and dynamic behaviour of such                                   the number of layers within the laminate. For very thick structures, a
structural members [1], whose design and material properties can in                                   three-dimensional elasticity theory is required to provide accurate re
fluence their mechanical response. Among innovative materials, lami                                   sults [7], whereas moderately thick shells commonly rely on an ESL
nated composites are largely used in several engineering fields due to                                 model [8–11]. Based on this approach, the three-dimensional kinematic
their lightweight nature and high stiffness and strength [2]. An accurate                              and mechanical quantities defined within each layer are referred to the
superposition of differently-oriented layers provides the best solution for                            geometric reference surface, according to some a-priori assumptions,
the specific engineering problem object of analysis. However, for com                                 thus reducing the problem to two-dimensions.
plex structural shapes, several issues must be considered for a proper                                     At the same time, classical methods set the original hypotheses
mechanical modelling, primarily an efficient geometrical description of                                focusing on the in-plane behaviour of plates and shells [12–15]. In this
arbitrary shapes [3–4].                                                                                way, warping and in-plane coupling effects along each direction of the
    With particular reference to displacement-based elastic theories, a                                structure are well described, especially if refined hypotheses are
two-dimensional problem can be set for a three-dimensional structure                                   assumed. In addition, in such theories the out-of-plane field variable
with fixed geometrical proportions in the structural thickness. More                                   component is assumed to be constant along the entire thickness of the
    * Corresponding author.
      E-mail address: francesco.tornabene@unisalento.it (F. Tornabene).
https://doi.org/10.1016/j.compstruct.2022.116542
Received 27 September 2022; Accepted 26 November 2022
Available online 9 December 2022
0263-8223/© 2022 Elsevier Ltd. All rights reserved.
F. Tornabene et al.                                                                                                         Composite Structures 309 (2023) 116542
structure. As a consequence, there are effects that are not predictable if        analytical closed-form solution can be found only for some specific
the lamination scheme assumes a softcore configuration [16–18].                   boundary conditions [73–74]. The Generalized Differential Quadrature
Moreover, shear stresses through-the-thickness profiles are discontin            (GDQ) method [75–78] has been proven to be one of the most powerful
uous, so that the shear correction factor is required in the formulation          algorithms belonging to this class [79]. Based on the Weierstrass inter
[19]. In this way, the macromechanical contribution of the actual stress          polation theory, it provides a quadrature weighted expression of the
distribution is considered.                                                       derivative of the unknown function with respect to the values assumed
    A milestone in the analysis of shells can be traced to the introduction       by the function in a discrete set of points selected from the definition
of the generalized approach for the description of the displacement field         domain. The success of the method lies on the definition of the weighting
[20–23]. According to this methodology, each component of the actual              coefficients and the assessment of the computational domain [80–81].
field variable is described in a unified manner by means of the so-called         Very accurate results can be obtained if non-uniform grid distributions
thickness functions, embedding the shape variation for each field                 are employed [82–84]. Moreover, there are several studies in which
component. In this way, several theories can be included in a unified             GDQ has been applied for the static [85–86] and dynamic [87–89]
model and a systematic set of analyses can be easily performed by simply          analysis of shell structures made of advanced materials. Generally
setting a different analytical expression for the thickness functions.            speaking, GDQ-based problems show a great accordance with respect to
Accordingly, the accuracy of an ESL theory can be easily assessed                 very refined solutions, even though a significantly reduced number of
starting from a proper definition of the through-the-thickness displace          Degrees of Freedom (DOFs) is employed. The GDQ method has been,
ment field, taking into account both polynomial [24–25] and trigono              also, successfully applied to problems involving singularities [90–91] as
metric [26–27] functions. Thus, a smooth continuous variation can be              well as time integration simulations [92]. Moreover, moving from such a
obtained along the thickness of the shell. In laminated structures,               procedure, it is also possible to perform integral operations [93] by
however, very complicated coupling effects can occur between two                  means of the Generalized Integral Quadrature (GIQ).
adjacent layers, and the displacement field can show a singularity. Thus,             When the unknown field variable is expressed from the interpolation
some thickness function must be set in order to account for such phe             of the values assumed in some specific point of a pre-determined
nomena, such as the so-called zigzag functions [28–30]. At the same               computational grid, the so-called weak formulation of the differential
time, the refined theories [31–32] account for the actual shear stiffness         problem comes out. In this way, the differentiation order of the solution
distribution in each layer of the laminate. In addition, the thickness            decreases. The classical Finite Element Method (FEM) follows this
functions can be introduced to account for some delamination issues at            approach, accounting for a decomposition of the physical domain [94].
the interlaminar stage [33–36].                                                   Unlike FEM, when the GDQ method is applied, a higher order smooth
    In the general case of doubly-curved structures, the fundamental              variation of the solution is guaranteed along the entire physical domain
governing equations are written in a curvilinear principal reference              [95]. On the other hand, when a weak formulation of the structural
system, following the well-known isogeometric approach [37–38]. For               problem is adopted, boundary conditions are defined from the simple
shells with arbitrary geometries, a distortion algorithm is introduced to         elimination of the DOFs related to the applied kinematic constraint. In
analyse them in their unmapped original shape. In the case of singu              contrast, strong form theories require some analytical expressions of the
larities like holes, a domain decomposition is essential [39–42] since a          external constraints which should be numerically implemented by
parametrization in principal coordinates is difficult to be achieved. All         means of the GDQ approach [96–98]. As a consequence, non-
the above-mentioned procedures are based on an isogeometric mapping               conventional external constraints can be expressed by means of a se
of the physical domain based on Non-Uniform Rational Basis Spline                 ries of pre-determined distributions of linear springs [99–102], whereby
(NURBS) curves [43–46] employed to define the structural edges. It                several configurations can be defined once the governing parameters are
should be remarked that classical distortion algorithms are based on a            properly calibrated [103–104]. In the present work, a general formu
polynomial description of the boundary edges in each computational                lation for the enforcement of external boundary conditions is presented.
domain [47].                                                                      A distribution of linear elastic springs is applied along each edge of
    Another interesting aspect regarding shell structural theories is the         arbitrarily shaped anisotropic shell structures, whose vibration is here
proper description of the mechanical behaviour for each layer. Due to             analysed according to a unified formulation based on the ESL strategy.
relatively recent manufacturing issues, many theories have been pro              Unlike other ESL formulations for the dynamic analysis of shell struc
posed in the literature, accounting for isotropic, non-homogeneous and            tures, in the present work a surface distribution of linear elastic springs
orthotropic media [47–51]. Most innovative materials, however, exhibit            is associated with each boundary of the three-dimensional solid and
a completely anisotropic behaviour, and require the introduction of               accounts not merely for an in-plane general dispersion of the springs, as
some coupling terms within the formulation, together with a series of             other formulations, but also for a generalized out-of-plane distribution
homogenization methods for a proper definition of the equivalent me              of linear elastic springs. In this way, a three-dimensional configuration
chanical properties. The interested reader can refer to some works                of the boundary conditions has been modelled. Furthermore, the cali
regarding to the application of Functionally Graded Materials [52–58],            bration of the position and shape parameters of the adopted distribu
Carbon Nanotubes [59–62], honeycomb and grid lattices [63–65] to                  tions allows one to constrain the structure in a specific point of its edges,
spatial structures.                                                               even though a spectral collocation method is adopted for the numerical
    Referring to the dynamic analysis of shells, the structural response          implementation of the differential problem.
can be successfully oriented thanks to a redistribution of the inertial               Following the ESL approach, the reference middle-surface is
masses by properly setting a thickness variation along the parametric             described by means of a parametrization in curvilinear principal co
lines [66–68]. To this purpose, very interesting results can be read for          ordinates, and isogeometric NURBS-based mapping is used to account
both linear static problems and free vibration analysis [69–70].                  for any possible domain distortion. A general relation is provided for the
    The governing equations can be easily derived from the application            thickness variation of each layer along the entire physical domain. A
of energy principles, depending on the nature of the problem. For a free          two-dimensional model is assessed starting from the computation of the
vibration problem, the Hamiltonian Principle can be set in a variational          elastic strain energy. The equilibrium equations are obtained from the
form, computing the elastic strain energy and inertia energy contribu            Hamiltonian Principle. A completely anisotropic constitutive relation is
tions. Following the integration by parts rule, a strong formulation of the       assigned to each layer of the stacking sequence, and the generalized
structural problem is obtained, without any reduction of the differen            elastic coefficients are defined for the entire lamination scheme. The
tiability order. Then, a boundary value problem is obtained from the              fundamental equations are derived with respect to the unknown vector
energy stationary component [71], that can be easily tackled numeri              of the generalized displacement field. The final relations are solved
cally with a spectral collocation method as reliable tool [72], because an        numerically by means of the GDQ method, providing a solution directly
                                                                              2
F. Tornabene et al.                                                                                                                            Composite Structures 309 (2023) 116542
in a strong form. The proposed theory is validated by a free vibration                 dinate system O α 1 , α 2 , α 3 is introduced, where axes α 1 and α 2 are
                                                                                                              ′
analysis of a mapped laminated plate with different external boundary                  defined from the principal directions of the reference surface, and α 3 =
constraints. Mode frequencies and shapes are successfully compared to                  ζ is taken parallel to the outward normal direction of the shell surface
those provided by a three-dimensional finite element approach. A                       (Fig. 1). In other words, any arbitrary point P of the shell is univocally
parametric analysis is also performed to check for the influence of some               identified from its position vector R(α 1 , α 2 , ζ), and it can be obtained as
boundary conditions key parameters on the dynamic response of those                    a shift of its projection P’ located on the reference surface, whose po
structures. Finally, a comprehensive set of modal analyses is performed                sition vector is denoted as r(α 1 , α 2 ), leading to the following relation
on singly-curved and doubly-curved panels with different lamination                    [22]:
schemes and external constraints, employing various Higher Order
                                                                                                                            h(α 1 , α 2 )
Shear Deformation Theories (HSDTs). The importance of the field var                   R(α 1 , α 2 , ζ) = r(α 1 , α 2 ) +                 z n(α 1 , α 2 )                        (1)
iable expansion order is pointed out, as well as the introduction of a                                                          2
zigzag assumption, for accurate predictions of complicated warping and                 where h(α1 , α2 ) is the shell thickness expression and z = 2ζ/h is a
stretching effects. All the analyses presented in the article, and the                 dimensionless through-the-thickness coordinate, such that z ∈ [ − 1, 1].
related formulation, are implemented in the DiQuMASPAB software                        Moreover, n(α 1 , α 2 ) denotes the midsurface outward normal unit vec
[105].                                                                                 tor. All these quantities are computed referring to each point (α 1 , α 2 ) of
                                                                                       the reference surface. The latter can be expressed following the pro
2. Equivalent single layer theory for a doubly-curved shell                            cedure reported in Ref. [22], taking into account the first order de
structure                                                                              rivatives of r(α 1 , α 2 ) with respect to in-plane coordinates α i = α 1 , α 2 ,
                                                                                       marked with r,i = ∂r/∂αi :
    In the present section an ESL fomulation is provided for the dynamic
analysis of laminated doubly-curved shell structures made of generally                          r,1 × r,2
                                                                                       n = ⃒⃒             ⃒                                                                      (2)
anisotropic materials. The geometry of the structure is described by                            r,1 × r,2 ⃒
using the main outcomes of differential geometry, whereas the equi                        For the sake of clarity, the vector product in Eq. (2) has been indi
librium equations are derived from the Hamiltonian principle. The                      cated with the symbol × . In order to provide a physical meaning to the
constitutive relationship of each layer of the structure is intended to be             geometric representation reported in Eq. (1), the in-plane coordinates
linear elastic, accounting for generally anisotropic materials.                                                                               [       ] [
                                                                                       variation must be restricted to a closed interval α01 , α11 × α02 , α12 ,
                                                                                                                                                                 ]
Fig. 1. Representation of a doubly-curved shell with respect to the reference middle surface according to the ESL description of the displacement field based on the
introduction of the thickness functions set along the principal directions of the structure. The physical domain mapping for the geometric description of an arbitrary-
shaped structure stems from a distortion algorithm based on NURBS curves.
                                                                                   3
F. Tornabene et al.                                                                                                                                                                  Composite Structures 309 (2023) 116542
                  ∑
                  l                               ∑
                                                  l                                                                            √̅̅̅̅̅̅̅̅̅̅̅̅             √̅̅̅̅̅̅̅̅̅̅̅̅
h(α 1 , α 2 ) =            h k (α 1 , α 2 ) =       (ζ k+1 (α 1 , α 2 ) − ζ k (α 1 , α 2 ) )                 (3)       A1 =     r,1 ⋅r,1 ,      A2 =      r,2 ⋅r,2                                                     (9)
                    k=1                           k=1
2, defined as:                                                                                                          ⎢ ⎥ ∑
                                                                                                                        ⎢ ⎥
                                                                                                                                  N +1 ⎢
                                                                                                                                       ⎢
                                                                                                                                                           ⎥⎢
                                                                                                                                                           ⎥ (τ) ⎥
                                                                                                                        ⎢ U2 ⎥ =       ⎢ 0      Fατ 2   0 ⎥⎢   u2 ⎥ ⎥ ⇔
                                                                                                                        ⎣ ⎦ τ=0 ⎣                          ⎦⎢
                                                                                                                                                            ⎣       ⎦
        α i − α0i                                                                                                                                       α3                                       (10)
ᾱi =                                                                                                                     U3              0       0    Fτ       (τ)
                                                                                                                                                               u3
        α1i − α0i              for i = 1, 2                                                                  (5)
                                                                                                                                                ∑
                                                                                                                                                N +1
αi = 1 − ᾱi
̃                                                                                                                       U(α 1 , α 2 , ζ, t) =          F τ (ζ)u(τ) (α 1 , α 2 , t)
     In the following, we describe the thickness variation expressions
                                                                                                                                                τ=0
ϕ f (α1 , α 2 ) for f = 1, ..., 4, as adopted in the present work [22]: where the unknown three-dimensional field has been expressed, for each
                     {                                                               {
                                       ¯αp1 1                                                       ¯αp2 2
ϕ 1 (α1 , α2 ) =                                           p1   , ϕ 2 (α1 , α2 ) =
                          (sin(π(n1 ᾱ1 + α1m ) ) )                                      (sin(π(n2 ᾱ2 + α2m ) ) )p 2
                     ⎧                                                               {                                                                                                                                 (6)
                     ⎨                    αp1 3
                                          ̃                                                          αp2 4
                                                                                                     ̃
ϕ 3 (α1 , α2 ) =                                                , ϕ 4 (α1 , α2 ) =
                                  α1 + α3m ) ) )p 3
                     ⎩ (sin(π (n3 ̃                                                                α2 + α4m ) ) )p 4
                                                                                         (sin(π(n4 ̃
                                                                                                                   4
F. Tornabene et al.                                                                                                                                       Composite Structures 309 (2023) 116542
    For a steady identification of the displacement field assumption, a                        each τ = 0, ..., N + 1, see Ref. [22]:
structural theory nomenclature is adopted. More specifically, the symbol
                                                                                                    N +1 ∑
                                                                                                    ∑    3                              N +1 ∑
                                                                                                                                        ∑    3                            N +1 ∑
                                                                                                                                                                          ∑    3
E refers to the ESL approach, D refers to the displacement-oriented                            ε=                 D ζ DαΩi F τ u(τ) =              Z (τ)α i DαΩi u(τ) =              Z (τ)α i ε(τ)α i
formulation, whereas N is the maximum order of the kinematic expan                                 τ= 0 i=1                            τ= 0 i=1                          τ= 0 i=1
sion adopted in each principal direction. When Z is introduced in the                                                                                                                                   (16)
nomenclature, the τ = (N + 1)-th order is added to the model, and the
                                                                                                    Thus, the three-dimensional strains ε(α1 , α2 , ζ, t) are related to their
thickness function defined in Eq. (11) is taken into account.
                                                                                               corresponding generalized expressions, based on the introduction of a
    As mentioned before, in an ESL framework, the position vector for a
                                                                                               kinematic matrix collecting the through-the-thickness derivatives of the
doubly-curved shell can be expressed from the midsurface r(α 1 , α 2 ), as
                                                                                               thickness functions set Fατ i (ζ) for each τ-th expansion order and for each
already defined in Eq. (1). Starting from the three-dimensional elasticity
                                                                                               principal direction α i = α 1 , α 2 , α 3 :
problem,         the       kinematic          strain    vector     ε(α 1 , α 2 , ζ, t) =                 ⎡ αi                                                        ⎤
[ ε 1 ε 2 γ 12 γ13         γ 23    ε 3 ]T can be expressed in principal curvilinear                        Fτ
                                                                                                                 0     0       0       0    0     0      0      0 ⎥
                                                                                                         ⎢H
coordinates [22]:                                                                                        ⎢ 1                                                         ⎥
                                                                                                         ⎢        αi                                                 ⎥
                                                 (             )                                         ⎢      F                                                    ⎥
                                                      ∑
                                                      3                                                  ⎢ 0      τ
                                                                                                                       0       0       0    0     0      0      0    ⎥
                                                                                                         ⎢                                                           ⎥
ε = DU = Dζ (DαΩ1 + DαΩ2 + DαΩ3 )U = Dζ                     DαΩi U                 (13)                  ⎢      H2                                                   ⎥
                                                                                                         ⎢             α        α                                    ⎥
                                                      i=1                                                ⎢           F   i
                                                                                                                              F   i                                  ⎥
                                                                                                         ⎢ 0     0     τ        τ
                                                                                                                                       0    0     0      0      0    ⎥
                                                                                                         ⎢            H1      H2                                     ⎥
where D is a three-dimensional kinematic operator. This operator can be                                  ⎢                                                           ⎥
                                                                                               Z (τ)α i
                                                                                                        =⎢                             αi           αi               ⎥    (17)
                                                                                                         ⎢                           F          ∂ F                  ⎥
split in a through-the-thickness part denoted with Dζ , which includes the                               ⎢ 0     0     0       0       τ
                                                                                                                                            0       τ
                                                                                                                                                         0      0    ⎥
                                                                                                         ⎢                            H1          ∂ζ                 ⎥
outward normal direction derivatives, and three operators DαΩi , i = 1, 2, 3                             ⎢
                                                                                                         ⎢
                                                                                                                                                                     ⎥
                                                                                                                                                                     ⎥
                                                                                                                                            αi            αi
accounting for the in-plane derivatives. Recalling the curvature thickness                               ⎢
                                                                                                         ⎢ 0     0     0       0       0
                                                                                                                                           Fτ
                                                                                                                                                  0
                                                                                                                                                       ∂Fτ
                                                                                                                                                                0
                                                                                                                                                                     ⎥
                                                                                                                                                                     ⎥
                                                                                                         ⎢                                                           ⎥
H i = 1 +ζ/R i for i = 1, 2 defined in Eq. (8), the differential operators                               ⎢                                 H  2         ∂ ζ          ⎥
                                                                                                         ⎢                                                           ⎥
introduced in Eq. (13) can be written in an extended form as [22]:                                       ⎣                                                    ∂Fτα i ⎦
                                                                                                            0    0     0       0       0    0     0      0
       ⎡                 ⎤       ⎡                     ⎤                                                                                                       ∂ζ
            1 ∂                          1 ∂A 1               ⎡             ⎤
       ⎢             0 0⎥        ⎢0                  0⎥                1                          The generalized strain vector ε(τ) (α 1 , α 2 , t) introduced in Eq. (16) is
       ⎢ A 1 ∂α1         ⎥       ⎢     A  A   ∂α       ⎥      ⎢ 0 0         ⎥
       ⎢                 ⎥
                                                                                               expressed, for each τ = 0, ..., N + 1, as follows:
                                         1 2     2
       ⎢                 ⎥       ⎢                     ⎥      ⎢       R1 ⎥
       ⎢ 1 ∂A 2          ⎥       ⎢        1   ∂        ⎥      ⎢             ⎥
       ⎢             0 0⎥        ⎢                     ⎥      ⎢        1 ⎥
                                 ⎢0                  0⎥       ⎢0 0          ⎥
       ⎢ A 1 A 2 ∂α1
       ⎢
                         ⎥
                         ⎥       ⎢
                                 ⎢
                                         A 2 ∂α2       ⎥
                                                       ⎥
                                                              ⎢
                                                              ⎢       R2 ⎥
                                                                            ⎥                  ε(τ) αi = D αΩi u(τ) for τ = 0, 1, 2, …, N, N + 1 i = 1, 2, 3                                            (18)
       ⎢                 ⎥       ⎢                     ⎥      ⎢             ⎥
       ⎢     1 ∂A1       ⎥       ⎢        1   ∂        ⎥      ⎢0 0
       ⎢−            0 0⎥        ⎢0                  0⎥                0 ⎥
       ⎢ A 1 A2 ∂α2
       ⎢
                         ⎥
                         ⎥       ⎢       A 1 ∂α1       ⎥
                                                              ⎢
                                                              ⎢
                                                                            ⎥
                                                                            ⎥                  2.3. Anisotropic elastic constitutive relation
       ⎢                 ⎥       ⎢                     ⎥      ⎢0 0     0    ⎥
       ⎢ 1 ∂             ⎥       ⎢         1   ∂A      ⎥      ⎢             ⎥
                                 ⎢                 2   ⎥      ⎢             ⎥
  α1   ⎢
DΩ = ⎢ A 2 ∂α2
                     0 0⎥    α2
                         ⎥, DΩ = ⎢
                                 ⎢
                                   0 −
                                        A 1 A2 ∂α1 ⎥
                                                     0    α
                                                       ⎥ DΩ = ⎢
                                                           3
                                                              ⎢0 0
                                                                    1 ∂ ⎥
                                                                            ⎥                     Based on the kinematic relations (16) within an ESL framework, it is
       ⎢                 ⎥
       ⎢        1        ⎥       ⎢
                                 ⎢0
                                                       ⎥
                                                       ⎥
                                                              ⎢
                                                              ⎢     A 1 ∂α1 ⎥
                                                                            ⎥                  possible to derive the constitutive elastic equations starting from the
       ⎢                 ⎥                  0        0
       ⎢ −           0 0⎥        ⎢
                                 ⎢
                                                       ⎥
                                                       ⎥
                                                              ⎢
                                                              ⎢     1 ∂ ⎥
                                                                            ⎥                  generalized Hooke’s law for each layer of the stacking sequence.
       ⎢       R1        ⎥       ⎢           1         ⎥      ⎢0 0          ⎥
       ⎢
       ⎢
                         ⎥       ⎢0       −          0⎥       ⎢     A 2 ∂α2 ⎥                  Referring to the k-th layer, the so-called material reference system
       ⎢      0      0 0⎥⎥       ⎢           R2        ⎥      ⎢             ⎥
                                 ⎢                     ⎥      ⎢        0 ⎥                     O α 1 α 2 ζ(k) is rotated with respect to the geometric principal reference
                                                                                                ′   (k) (k)
       ⎢                 ⎥       ⎢0                    ⎥      ⎢0 0          ⎥
       ⎢      1      0 0⎥⎥       ⎢          0        0 ⎥      ⎢             ⎥
       ⎢                                                                                       system O α 1 α 2 ζ. Thus, each layer of the stacking sequence features the
                                                                                                              ′
       ⎢                 ⎥       ⎢                     ⎥      ⎢             ⎥
       ⎢                 ⎥       ⎢                     ⎥      ⎢0 0     0 ⎥
       ⎢      0      0 0 ⎥       ⎢
                                 ⎣
                                   0        1        0 ⎥
                                                       ⎦
                                                              ⎣             ⎦                  following anisotropic behaviour [22]:
       ⎣                 ⎦                                      0 0    1
              0      0 0           0        0        0
(14)
                                                                                           5
F. Tornabene et al.                                                                                                                                                               Composite Structures 309 (2023) 116542
⎡             ⎤         ⎡                                                    ⎤⎡          ⎤
                                                                                                                               The stiffness matrix C̄          is now referred to the three-dimensional
                                                                                                                                                                 (k)
        (k)         (k)               (k)  (k)  (k)  (k)  (k)                      (k)
⎢̂σ           ⎥ ⎢ C11                C12  C16  C14  C15  C13  ⎥⎢ ̂ε                      ⎥                                                                          [                         ]T
⎢ 1           ⎥ ⎢                                             ⎥⎢ 1                       ⎥                                  stress and strain vectors σ = σ 1 σ 2 τ 12 τ 13 τ 23 σ 3             and ε(k) =
                                                                                                                                                                      (k) (k) (k) (k) (k)
                                                                                                                                                                       (k)                (k)
⎢̂   (k)      ⎥ ⎢ (k)                 (k)  (k)  (k)  (k)  (k) ⎥⎢ (k)                     ⎥
⎢ σ2          ⎥ ⎢ C12                C22 C26 C24 C25 C23 ⎥⎢ ̂     ε                      ⎥                                  [                            ]T
⎢ (k)         ⎥ ⎢ (k)                                         ⎥⎢ 2                       ⎥
⎢ ̂τ
⎢ 12
              ⎥ ⎢C
              ⎥ ⎢ 16
                                      (k)
                                     C26   (k)
                                          C66   (k)
                                               C46   (k)
                                                    C56   (k) ⎥⎢ (k)
                                                         C36  ⎥⎢ ̂γ 12
                                                                                         ⎥
                                                                                         ⎥                  (k)               ε(k) (k) (k) (k) (k) (k)
                                                                                                                                1 ε 2 γ 12 γ 13 γ 23 ε 3    , respectively. Each component of Eq. (20) is
⎢ (k)         ⎥ = ⎢ (k)                                       ⎥⎢                             ̂ (k) = C(k) ̂
                                                                                         ⎥ ⇔ σ            ε
⎢ ̂τ          ⎥ ⎢C                                                                                                          denoted with¯E ij in order to easily refer to the rotated three-dimensional
                                      (k)  (k)  (k)  (k)  (k) ⎥⎢ (k)                     ⎥                                                         (k)
⎢ 13          ⎥ ⎢ 14                 C24  C46  C44  C45  C34  ⎥⎢ ̂γ 13                   ⎥
⎢ (k)         ⎥ ⎢ (k)                                     (k) ⎥⎢ (k)                     ⎥
                                                                                                                            stiffness coefficients ¯C(k) ij , and to reduced terms ¯Q ij valid in the case of
                                      (k)  (k)  (k)  (k)
⎢ ̂τ 23       ⎥ ⎢ C15                C25  C56  C45  C55  C35  ⎥⎢ ̂γ 23                   ⎥                                                                                                 (k)
⎢             ⎥ ⎢                                             ⎥⎢                         ⎥
⎣ (k)         ⎦ ⎣ (k)                                     (k) ⎦⎣ (k)                     ⎦
                                                                                                                            plane-stress assumption. Starting from the computation of the elastic
                                      (k)  (k)  (k)  (k)
  σ3
  ̂                C13               C23 C36 C34 C35 C33         ̂ε3
                                                                                                                            strain energy, taking into account all the constitutive relations reported
                                                                                                              (19)          in Eqs. (19)–(20), the generalized ESL constitutive relation is derived in
                                         [                                        ]T                                        terms of the generalized stress resultant vector S(τ)α i (α1 , α2 , t) =
being             ̂ (k) = ̂
                  σ          σ (k)  σ (k) τ (k) τ (k) τ (k) σ (k) and       ̂ε (k) =                                         [                                                                 ]T
                                1 ̂    2 ̂   12 ̂  13 ̂  23 ̂  3                                                              N(1) i N(2) i N(12) i N(21) i T(1) i T (2) i P(1) i P(2) i S(3) i , as detailed in
                                                                                                                                 τα    τα     τα      τα       τα     τα     τα     τα     τα
[                                     ]T
  ε (k) ε (k)                              the three-dimensional stress and strain                                          Ref. [22]:
               (k) (k) (k) (k)
  ̂  1 ̂   2 ̂γ 12 ̂γ 13 ̂γ 23 ̂ε 3
vectors, respectively, and                             C ij   for i, j = 1, ..., 6 the elastic stiffness con
                                                        (k)                                                                            N +1 ∑
                                                                                                                                       ∑    3
                                                                                                                            S(τ)αi =               A(τη)αi αj ε(η)αj for τ = 0, 1, 2, ..., N, N + 1, αi = α1 , α2 , α3
stants of the k-th layer referred to O                   Once the elastic      1 α2 ζ .
                                                                             α(k)
                                                                         ′        (k) (k)
                                                                                                                                       η=0   j=1
constitutive relationship is expressed for the entire laminate, a homog                                                                                                                                             (21)
enization technique should be developed in order to obtain a single
relation valid for the entire structure. In the present work, setting ζ(k) =                                                      In Eq. (21), a generalized 9 × 9 stiffness matrix A      comes out for
                                                                                                                                                                                                 (τη)α i α j
                                                                                                                            each τ, η-th order of the kinematic assumption (10) and for each α i α j =
ζ a rotation transformation T , characterized by an angle θ k between
                                                          (k)
                                                                                                                            α 1 , α 2 , α 3 principal direction. In an extended matrix form, one gets:
  1 and α 1 , can be defined. In this way, the stiffness matrix C
α(k)                                                                 in Eq.
                                                                 (k)
(22)
O α 1 α 2 ζ defined for the shell under consideration [22]:                                                                     If we introduce the notations Fταi = ∂0 Fταi /∂ζ0 and Fη j = ∂0 Fη j /∂ζ0 that
    ′                                                                                                                                                                                               α          α
                        ⎡                                 ⎤                                                                 allow us to indicate and identify that the 0-th order derivatives of the
                                                                                                                            thickness functions Fταi , Fτ j coincide with the thickness functions them
                                                                                                                                                        α
                                           ⎢¯E(k)       ¯E(k)   (k)   (k)   (k)   (k)
                                                          12 ¯E 16 ¯E 14 ¯E 15 ¯E 13 ⎥
                                           ⎢ 11                                       ⎥                                     selves, the generic component of the generalized stiffness matrix
                                           ⎢ (k)                                      ⎥
                                           ⎢¯E12        ¯E22 ¯E26 ¯E24 ¯E25 ¯E(k)
                                                          (k)   (k)   (k)   (k)
                                                                                  23 ⎥
                                           ⎢ (k)
                                           ⎢¯E
                                                                                      ⎥
                                                                                  (k) ⎥
                                                                                                                            A(τη)α i α j can be expressed by the following condensed relation, accord
                                           ⎢            ¯E(k)   (k)   (k)   (k)
                                                          25 ¯E 66 ¯E 46 ¯E 56 ¯E 36 ⎥                                      ing to the nomenclature introduced in Eq. (22):
                                                                                                              (20)
     (k)
C̄         =T C T (k)       (k)   (k)T
                                         = ⎢ 16                                       ⎥
                                           ⎢¯E(k)       ¯E(k)   (k)   (k)   (k)   (k) ⎥
                                           ⎢ 14           24 ¯E 46 ¯E 44 ¯E 45 ¯E 34 ⎥
                                           ⎢ (k)                                  (k) ⎥
                                           ⎢¯E15
                                           ⎢            ¯E(k)   (k)   (k)   (k)
                                                          25 ¯E 56 ¯E 45 ¯E 55 ¯E 35 ⎥
                                                                                      ⎥
                                           ⎣ (k)          (k)   (k)   (k)   (k)   (k) ⎦
                                            ¯E13        ¯E23 ¯E36 ¯E34 ¯E35 ¯E33
                                                                                               for τ, η = 0, 1, 2, ..., N, N + 1
                                                                                               for n, m = 1, 2, 3, 4, 5, 6
                                         ∫
                             ∑l              ζk+1         ∂f Fηα j ∂g Fτα i H1 H2
                                                                                                                                                                                                                     (23)
 (τη)[fg]α α
A nm (pq)i j            =                           ¯B(k)                    p    q dζ         for p, q = 0, 1, 2
                                   k=1
                                          ζk
                                                       nm
                                                           ∂ζf ∂ζg (H1 ) (H2 )
                                                                                               for αi , αj = α1 , α2 , α3
                                                                                               for f , g = 0, 1
                                                                                                                      6
F. Tornabene et al.                                                                                                                                                                            Composite Structures 309 (2023) 116542
where H k for k = 1, 2 is defined in Eq. (8). In Eq. (23) one can find the                                  f (τη)α i        (− )                                  (− )      (− )          (+)                        (+)   (+)
generic τ, η-th order of the field variable expansion, as well as the indexes                           L ii             = kif f (− ) Fηαi (− ) Fταi (− ) H1 H2 + kif f (+) Fηαi (+) Fταi (+) H1 H2                                  for
n, m of the stiffness coefficients referred to the k-th layer. The latter has                           i = 1, 2, 3
been denoted with ¯B(k)
                    nm so that Eij                   =¯Cij , or Eij         =¯Qij . According to                                                                                                                                       (28)
                                               (k)        (k)        (k)         (k)
 Eq. (23), also the shear correction factor κ(ζ) can be embedded in the                                    According to Eq. (28), a variation of linear springs can be assessed on
 model when lower order theories cannot predict the actual through-the-                                 the top and bottom surfaces of the shell, from the enforcement of two
 thickness shear stresses distribution [22], namely:                                                                             (         )
                                                                                                        bivariate functions f (±) ξ 1 , ξ 2 . If a constant Winkler distribution of
        {
           ¯E(k)    for n, m = 1, 2, 3, 6                                                               springs with stiffnesses k(±)
                                                                                                                                  if is defined on the top and bottom surface, the
¯B(k) =       nm
                                                                   (24)
  nm
          κ(ζ)¯E(k)
                 nm
                    for n, m = 4, 5                                                                     function at issue should be defined at each point of the physical domain,
                                                                                                        such that:
    The generalized ESL constitutive relation reported in Eq. (21) can be                                    (         )
rearranged in such a way that each component of the stress resultant                                    f (±) ξ 1 , ξ 2 = 1                                                  (29)
vector S(τ)α i (α1 , α2 , t), for each τ = 0, ...N + 1, can be directly expressed                          Moreover, a Gaussian two-dimensional spring distribution can be
in terms of the generalized kinematic vector u(η) (α1 , α2 , t), with η = 0,...,                                                                                          (          )
                                                                                                        applied to the external surfaces of the shell. In this case, f (±) ξ 1 , ξ 2
N + 1:                                                                                                  reads as follows:
⎡       ⎤           ⎡                                ⎤
                                                                                                                                   ((      )2 (       )2 )
                                                                                                                                                                                  (±)                    (±)
⎢ (τ)αi     ⎥     ⎢ (τη)αi α1                                     ⎥                                               (         )                              −   1
                                                                                                                                                                          ξ1 − ξ
                                                                                                                                                                                   1m      +
                                                                                                                                                                                                 ξ2 − ξ
                                                                                                                                                                                                          2m
⎢N
⎢ 1
            ⎥
            ⎥
                  ⎢O
                  ⎢ 11              O(21τη)αi α2     O(31τη)αi α3 ⎥
                                                                  ⎥                                     f   (±)
                                                                                                                      ξ 1, ξ 2 =
                                                                                                                                            1
                                                                                                                                                       e
                                                                                                                                                               2
                                                                                                                                                                            Δ
                                                                                                                                                                                (±)
                                                                                                                                                                                 1
                                                                                                                                                                                                   Δ
                                                                                                                                                                                                       (±)
                                                                                                                                                                                                        2                              (30)
⎢ (τ)αi     ⎥     ⎢ (τη)αi α1                           τη α α ⎥                                                                        2πΔ(±) (±)
                                                                                                                                            1 Δ2
                                       τη αi α2
⎢ N2
⎢
            ⎥
            ⎥
                  ⎢ O12
                  ⎢                 O(22 )           O(32 ) i 3 ⎥ ⎥
⎢ (τ)αi     ⎥     ⎢ (τη)αi α1          τη αi α2         τη α α ⎥
⎢ N12       ⎥     ⎢ O13             O(23 )           O(33 ) i 3 ⎥⎡          ⎤                           where Δ 1 , Δ 2 > 0 and ξ 1m , ξ 2m > 0 are the position and scaling pa
                                                                                                                          (±)       (±)                        (±)         (±)
⎢           ⎥     ⎢                                               ⎥
⎢ N (τ)αi   ⎥     ⎢ O(τη)αi α1         τη α α
                                    O(24 ) i 2       O(34 ) i 3 ⎥
                                                        τη α α        (η)
⎢ 21
⎢ (τ)α
            ⎥ ∑   ⎢
            ⎥ N+1 ⎢ 14
                                                                  ⎥⎢ u1
                                                                  ⎥         ⎥                           rameters of the distribution, referred to the top ( + ) and bottom (− )
                                                       (τη)αi α3 ⎥⎢ (η)     ⎥
⎢T i
⎢ 1
            ⎥=
            ⎥
                  ⎢ O(τη)αi α1
                  ⎢ 15
                                      (τη)αi α2
                                    O25              O35          ⎥⎢ u      ⎥                (25)       surface of the shell, for each parametric line of the computational
⎢ (τ)αi     ⎥ η=0 ⎢ (τη)αi α1                                      ⎣ 2      ⎦
                                       τη α α           τη α α ⎥      (η)                               domain.
⎢ T2
⎢
            ⎥
            ⎥
                  ⎢ O16
                  ⎢                 O(26 ) i 2       O(36 ) i 3 ⎥ ⎥ u3
⎢ (τ)αi
⎢ P1
            ⎥
            ⎥
                  ⎢ (τη)αi α1
                  ⎢ O17               (τη)α α          (τη)α α    ⎥                                        Linear elastic boundary springs can be applied on the edge faces of
                                    O27 i 2          O37 i 3 ⎥
⎢
⎢ (τ)αi
            ⎥
            ⎥
                  ⎢
                  ⎢ (τη)αi α1         (τη)α α          (τη)α α ⎥
                                                                  ⎥                                     the structure by means of predefined distributions of normal and shear
⎢ P2        ⎥     ⎢ O18             O28 i 2          O38 i 3 ⎥                                                                                                                                                                    (k)αmj
⎢ (τ)α
⎢S i
            ⎥
            ⎥
                  ⎢
                  ⎢ O(τη)αi α1
                                                                  ⎥                                     stresses along each boundary edge (Fig. 2). Let us denote with k if                                                                for
⎣ 3         ⎦     ⎣ 19              O(29τη)αi α2     O(39τη)αi α3 ⎥
                                                                  ⎦
                                                                                                        i = 1, 2, 3, j = 1, 2 and m = 0, 1 the stiffness of the boundary linear
                                                                                                        springs at issue in each principal direction. The following definitions can
                                                                                                                                                           (      )
                                                                                                        be introduced, accounting for their variation f αm1 , α 2 along the shell
2.4. Winkler-type foundation and general boundary conditions                                                                                                     [       ]
                                                                                                        boundaries, setting α 1 = αm1 , with m = 0, 1 and α 2 ∈ α02 , α12 :
    In the present work, a general distribution of surface linear elastic                               ¯σ(k)
                                                                                                              ( m          )     (k)αm1 ( m         ) (
                                                                                                                                         f α 1 , α 2 U 1 αm1 , α 2 , ζ
                                                                                                                                                                       )
                                                                                                           1 α 1 , α 2 , ζ = − k 1f
springs is applied on the top and the bottom surface of the shell.                                            (            )              (         )   (              )
                                                                                                                                                                                                                                       (31)
                                                                                                                                      m
                                                                                                                                 (k)α 1
Furthermore, a general distribution of springs is considered for each                                   ¯τ(k)   m
                                                                                                           12 α 1 , α 2 , ζ = − k 2f     f αm1 , α 2 U 2 αm1 , α 2 , ζ                                         for m = 0, 1
                                                                                                              (            )          m (           ) (                )
lateral surface of the three-dimensional doubly-curved solid. In this way,                              ¯τ(k)   m                (k) α
                                                                                                                                         f αm1 , α 2 U 3 αm1 , α 2 , ζ
                                                                                                           13 α 1 , α 2 , ζ = − k 3f
                                                                                                                                       1
(27)
            f(τη)α i
where L ii             for i = 1, 2, 3 can be expressed as:
                                                                                                    7
F. Tornabene et al.                                                                                                                                                                 Composite Structures 309 (2023) 116542
  (k) (            )     (k)αm (           )∑N+1 α      (η) (        )                                          dimensional ESL model. Taking into account Eqs. (33)–(34), the
¯τ 12 α 1 , αm2 , ζ = − k 1f 2 f α 1 , αm2      F 1 (ζ)u 1 α 1 , αm2
                                             η=0 η
                                                                                                                following compact relation is assessed, following the procedure from
       (           )     (k)αm2 (          )∑N+1 α       η (           )                                        Ref. [22]. Thus, the generalized boundary conditions for α 1 = αm1 for
¯σ (k)       m
    2 α 1 , α 2 , ζ = − k 2f   f α 1 , αm2      F 2 (ζ)u(2 ) α 1 , αm2
                                             η=0 η
                                                                                             for m = 0, 1                           [            ]
                                                                                                                m = 0, 1 and α 2 ∈ α02 , α12 , are defined as:
   (k) (            )            (k)αm2    (             )∑N+1               (η) (       )
¯τ 23 α 1 , αm2 , ζ = − k 3f              f α 1, α 2   m
                                                                   Fαη 3 (ζ)u 3 α 1 , α 2
                                                                                       m
                                                                                                                                    ⎡                                        ⎤
                                                             η=0                                                ⎡          ⎤           (τη)α                                   ⎡      ⎤
                                                                                                                                      L f 1(2)α1 m      0            0
                                                                                                    (34)         ¯ N
                                                                                                                ⎢ 1 ⎥
                                                                                                                     (τ)α1          ⎢
                                                                                                                                    ⎢            1                           ⎥ u(1η)
                                                                                                                                                                             ⎥ ⎢      ⎥
                                                                                                                ⎢ (τ)α2 ⎥      ∑⎢
                                                                                                                               N +1
                                                                                                                                                                             ⎥⎢ (η) ⎥
                                                                                                                                                                                        (40)
                                                                                                                                                    (τη)α
    In the present work, the general boundary conditions along the edges                                        ⎢¯N 12 ⎥ = −
                                                                                                                ⎣          ⎦
                                                                                                                                    ⎢ 0
                                                                                                                                    ⎢
                                                                                                                                                   L f 2(2)α2 m
                                                                                                                                                              1
                                                                                                                                                                     0       ⎥⎢ u 2 ⎥
                                                                                                                                                                             ⎥⎣       ⎦
of the structure have been assessed in terms of two different distribu                                              (τ)α3      η=0 ⎣                                        ⎦    (η)
                                                                                                                  ¯T 1                                           (τη )α          u3
                                                                                                                                           0            0       L f 3(2)α2 m
tions of linear springs, namely a Double-Weibull and a super elliptic one.                                                                                                                 1
kinematic expansion order at the edges characterized by α 1 = αm1 , for                                            More specifically, the kinetic term δT is assessed taking into account
                     [      ]
m = 0, 1 and α 2 ∈ α02 , α12 . If l is the number of laminae occurring in the                                   the inertial contribution of each k-th layer of a l laminae stacking
stacking sequence, one has [22]:                                                                                sequence to the overall structure [22]:
                                 ∫                                                                                       l ∫
                                                                                                                         ∑        ζk+1   ∫ ∫
         (         ) ∑l               ζk+1
¯N (1τ)α1 αm1 , α 2 =                        ¯σ (k)   α1
                                                 1 λ̄F τ H2 dζ
                                                                                                                δT =                                ρ(k) (δU)T Ü H 1 H 2 A 1 A 2 dα 1 dα 2 dζ                       (44)
                           k=1       ζk                                                                                  k=1     ζk        α1 α2
                           l ∫
                           ∑          ζk+1
    τ α2 (          )
¯N (12)      αm1 , α 2 =                       ¯τ(k)    α2
                                                  12 λ̄F τ H2 dζ                                    (38)        where ρ(k) is the density associated with each k-th layer of the stacking
                           k=1       ζk                                                                                            [            ]T
                                 ∫                                                                              sequence and Ü = Ü 1 Ü 2 Ü 3 the second order time derivative of the
         (         ) ∑l               ζk+1
¯T (1τ)α3 αm1 , α 2 =                        ¯τ(k)    α3
                                                13 λ̄F τ H2 dζ
                           k=1     ζk                                                                           three-dimensional displacement field vector. The virtual variation δΦ of
                                                                                                                the elastic strain energy can be stored in Eq. (43) accounting for the
   In the same way, the boundary stresses can be computed on the edges
                                     [         ]                                                                generalized stress resultant S(τ)α i (α1 , α2 , t) and the higher order ESL
with α 2 = αm2 for m = 0, 1 and α 1 ∈ α01 , α11 as follows [22]:
                                                                                                                strains ε(τ)α i (α1 , α2 , t) for each τ = 0, ..., N +1 and αi = α 1 , α 2 , α 3 :
                                 ∫
    τ)α1 (          ) ∑l              ζk+1
                                                                                                                         ∑    3 ∫ ∫
                                                                                                                         N +1 ∑
¯N (21     α 1 , αm2 =                         ¯τ(k)    α1                                                                          (                          )T
                                                  12 λ̄Fτ H1 dζ                                                 δΦ =                                δε(τ)αi         S(τ)αi A1 A2 dα1 dα2                             (45)
                           k=1       ζk
                                                                                                                         τ=0     i=1     α1 α2
                           l ∫ ζk+1
                           ∑
   (τ)α2 (          )
¯N 2         α 1 , αm2 =                         (k)
                                             ¯σ 2 λ̄Fτα2 H1 dζ                                      (39)        The complete demonstration of Eq. (45) can be found in Ref. [22].The
                                     ζk
                                                                                                                energetic contribution δL f of the external loads consists of a part related
                           k=1
                                 ∫
         (         ) ∑l               ζk+1                                                                                                                     [                 ]T
¯T (τ2)α3 α 1 , αm2 =                        ¯τ(k)    α3
                                                23 λ̄Fτ H1 dζ                                                   to the top and bottom surface actions q fk = q 1 fk q 2 fk q 3 fk reported
                                                                                                                                                         (τ )    (τ ) (τ )  (τ )
                           k=1     ζk
                                                                                                                in Eq. (27) related to the Winkler foundation, and another one ac
       The generalized stress resultants of Eqs. (38)–(39) can be expressed                                     counting for boundary stresses (38), (39), according to Eq. (40) [22]:
in terms of the generalized displacement components u(iη) for i = 1, 2, 3
and η = 0,...,N + 1. In this way, they can be easily embedded in the two-
                                                                                                            8
F. Tornabene et al.                                                                                                                                                                            Composite Structures 309 (2023) 116542
          N +1 ∫ ∫ (
          ∑                                                       )                                                      ⎡                                                                    ⎤
                      τ        τ       τ       τ       τ       τ
δL f =              q(1 )fk δu(1 ) + q(2fk) δu(2 ) + q(3fk) δu(3 ) A 1 A 2 dα 1 dα 2 +                                  I 0(τη)α 1 α 1                     0                       0
          τ=0    α1 α2                                                                                                ⎢                                                                       ⎥
                                                                                                           M(τη)     =⎢       0                      I 0(τη)α 2 α 2                0          ⎥ for τ, η = 0, ..., N + 1           (49)
        N +1 ∮ (
        ∑                                               )                                                             ⎣                                                                       ⎦
                  (τ)α    (τ)   (τ)α   (τ)   (τ)α   (τ)                                                                                                                        0(τη)α 3 α 3
      +         ¯N 21 1 δu 1 +¯N 2 2 δu 2 +¯T 2 3 δu 3 A 1 dα 1 +                                                                   0                      0               I
                                                                                                (46)
          τ=0
                α1
          N +1 ∮ (
          ∑                                               )                                                where each component I0(τη)α i α j for i, j = 1, 2, 3 can be computed as:
                    (τ)α   (τ)   (τ)α    (τ)   (τ)α   (τ)
      +           ¯N 1 1 δu 1 +¯N 12 2 δu 2 +¯T 1 3 δu 3 A 2 dα 2
                                                                                                                               l ∫
                                                                                                                               ∑            ζk+1
          τ=0
                α2
                                                                                                           I 0(τη)α i α j =                           ρ(k) Fτα i Fηα j H1 H2 dζ                for α i , α j = α 1 , α 2 , α 3 ,
    In the previous relation, the two last integrals account for an inte                                                       k=1        ζk
gration along the curved edges of the physical domain characterized by                                     τ, η = 0, ..., N + 1
α 2 and α 1 constant, respectively, whereas the first one is a surface in                                                                                                                                                         (50)
tegral referred to the entire physical domain. Employing the virtual                                          Introducing in Eq. (47) the generalized constitutive relationship
variation of the kinetic energy derived in Eq. (44), the elastic contri                                   (21), as well as the ESL kinematic relation reported in Eq. (16), the
bution of Eq. (45), and boundary springs reported in Eq. (46) in the                                       application of the variational theorems leads to the final form of the
Hamiltonian Principle of Eq. (43), the equilibrium relations for an                                        fundamental equations:
anisotropic doubly-curved shell are derived in the following compact
form [22]:                                                                                                 ∑
                                                                                                           N +1
                                                                                                                (                  )       ∑
                                                                                                                                           N +1
                                                                                                                    L(τη) − Lf (τη) u(η) =      M(τη) ü(η)                                   for τ = 0, ..., N + 1                (51)
∑
3                              ∑
                               N +1                                                                        η=0                                              η=0
      D*Ωαi S(τ)αi + q fk =                                                                     (47)
                        (τ)
                                      M(τη) ü(η) for τ = 0, ..., N + 1
i=1                            η=0
                                                                                                           where Lf(τη) is the fundamental matrix related to the surface Winkler
             [                 ]T                                                                          elastic support, accounted in Eq. (27) for any τ, η-th kinematic expansion
where ü(η) = ü 1 ü 2 ü 3
                (η ) (η ) (η )
                                  is the second order time derivative of the                               order. The elastic fundamental differential matrix L(τη) , for each
                                                                                                           τ, η = 0, ..., N +1 order, assumes the following expanded form [22]:
generalized displacement field u(η) , for each η = 0, ..., N + 1, whereas
D*Ωαi for i = 1, 2, 3 stand for the following equilibrium operators of di
mensions 3 × 9, according to Ref. [22]:
           ⎡                                     ⎤T              ⎡                                ⎤T
                                                                                                                     ⎡                                                ⎤T
           ⎢                                    ⎥                ⎢                               ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢ 1         ∂       1 ∂A2            ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢              +                0   0⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢ A1      ∂α1 A1 A2 ∂α1              ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢0
           ⎢
           ⎢               1 ∂A 2
                                                ⎥
                                                ⎥                ⎢        −     1 ∂A 1
                                                                              A1 A2 ∂α2
                                                                                                0⎥
                                                                                                 ⎥
                                                                                                                 ⎢
                                                                                                                 ⎢0            0                 −    1
                                                                                                                                                      R1
                                                                                                                                                                     ⎥
                                                                                                                                                                     ⎥
           ⎢         −                     0   0⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢             A1 A2 ∂α1              ⎥                ⎢     1 ∂                ∂A1    ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢0
                                                                 ⎢     A2 ∂α2
                                                                                + A11A2   ∂α2   0⎥
                                                                                                 ⎥
                                                                                                                 ⎢0
                                                                                                                 ⎢             0                 −    1
                                                                                                                                                      R2
                                                                                                                                                                     ⎥
                                                                                                                                                                     ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢              1 ∂A1                 ⎥                ⎢0
           ⎢                               0   0⎥                ⎢
                                                                       1 ∂
                                                                                + A11A2   ∂A2
                                                                                                0⎥
                                                                                                 ⎥
                                                                                                                 ⎢0
                                                                                                                 ⎢             0                      0              ⎥
                                                                                                                                                                     ⎥
           ⎢           A1 A2 ∂α2                ⎥                ⎢     A1 ∂α1             ∂α1
                                                                                                 ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢0            0                      0              ⎥
           ⎢                                    ⎥                ⎢0             1 ∂A2
                                                                                                0⎥               ⎢                                                   ⎥
           ⎢ 1         ∂       1 ∂A1            ⎥                ⎢            A1 A2 ∂α1          ⎥               ⎢                                                   ⎥
           ⎢              +                0   0⎥                                                                                                                                                                                  (48)
D*Ωα 1   = ⎢ A2      ∂α2 A1 A2 ∂α2              ⎥      D*Ωα 2   =⎢
                                                                 ⎢0              0
                                                                                                 ⎥
                                                                                                0⎥     D*Ωα 3   =⎢
                                                                                                                 ⎢0            0      1 ∂
                                                                                                                                      A1 ∂α1
                                                                                                                                                 + A11A2         ∂A2 ⎥
                                                                                                                                                                 ∂α1 ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢             1                      ⎥                ⎢                               ⎥               ⎢0                                              ∂A1 ⎥
           ⎢                                    ⎥                ⎢0              1
                                                                                                0⎥               ⎢             0      1 ∂
                                                                                                                                                 + A11A2
           ⎢                               0   0⎥                                                                                     A2 ∂α2                     ∂α2 ⎥
           ⎢             R1                     ⎥                ⎢               R2              ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                                     ⎢0              0              0⎥               ⎢0            0                      0              ⎥
           ⎢              0                0   0⎥
                                                ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                                     ⎢0             − 1             0⎥               ⎢0            0                      0              ⎥
           ⎢             − 1               0   0⎥
                                                ⎥                ⎢
                                                                 ⎢
                                                                                                 ⎥
                                                                                                 ⎥
                                                                                                                 ⎢
                                                                                                                 ⎢
                                                                                                                                                                     ⎥
                                                                                                                                                                     ⎥
           ⎢                                    ⎥                ⎢0              0              0⎥               ⎢0            0                 − 1                 ⎥
           ⎢              0                0   0⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢              0                0   0⎥                ⎢                               ⎥               ⎢                                                   ⎥
           ⎢                                    ⎥                ⎢                               ⎥               ⎣                                                   ⎦
           ⎢                                    ⎥                ⎣                               ⎦
           ⎣                                    ⎦
                                                                                                       9
F. Tornabene et al.                                                                                                                                                          Composite Structures 309 (2023) 116542
                                                                               ⎡                                                                                    ⎤⎞
⎛⎡                                                                        ⎤  ⎢                                                                                      ⎥ ⎟⎡ (0)       ⎤
       L(00)          L(01)       ⋯       L(0(N) )       L(0(N+1) )          ⎢ Lf (00)                Lf (01)       ⋯         Lf (0(N) )           Lf (0(N+1) )     ⎥⎟ u
⎜⎢                                                                        ⎥ ⎢                                                                                       ⎥ ⎟⎢ (1)       ⎥
⎜⎢ L(10)              L(11)       ⋯       L(1(N) )       L(1(N+1) )       ⎥ ⎢   Lf (10)               Lf (11)       ⋯         Lf (1(N) )           Lf (1(N+1) )     ⎥ ⎟⎢ u         ⎥
⎜⎢                                                                        ⎥ ⎢                                                                                       ⎥ ⎟⎢           ⎥
⎜⎢     ⋮               ⋮          ⋱          ⋮                ⋮           ⎥− ⎢       ⋮                   ⋮          ⋱             ⋮                     ⋮
                                                                                                                                                                    ⎥ ⎟⎢
                                                                                                                                                                    ⎥ ⎟⎢ ⋮         ⎥=
⎜⎢                                                                        ⎥ ⎢                                                                                                      ⎥
⎜⎢ ((N)0 )                                                                ⎥ ⎢                                                                                       ⎥ ⎟⎢           ⎥
⎝⎣ L                L((N)1 )      ⋯      L((N)(N) )     L((N)(N+1) )      ⎦ ⎢⎢ L
                                                                                  f ((N)0 )
                                                                                                     Lf ((N)1 )     ⋯        Lf ((N)(N) )         Lf ((N)(N+1) )    ⎥ ⎟⎣ u(N)
                                                                                                                                                                    ⎥⎟             ⎦
                                                                             ⎢ f ((N+1)0 )                                                                          ⎥ ⎟ (N+1)
   L((N+1)0 )      L((N+1)1 )     ⋯     L((N+1)(N) )   L((N+1)(N+1) )        ⎣L                  Lf ((N+1)1 )       ⋯       Lf ((N+1)(N) )     Lf ((N+1)(N+1) )     ⎦⎠ u
                                                                                                                                                                                                             (53)
                                                                      ⎡                                                                                  ⎤⎡              ⎤
                                                                     M(00)               M(01)            ⋯         M(0(N) )           M(0(N+1) )           ü(0)
                                                                   ⎢                                                                                     ⎥⎢ (1)          ⎥
                                                                   ⎢ M(10)               M    (11)
                                                                                                          ⋯         M   (1(N) )
                                                                                                                                       M    (1(N+1) )
                                                                                                                                                         ⎥⎢ ü           ⎥
                                                                   ⎢                                                                                     ⎥⎢              ⎥
                                                                  =⎢
                                                                   ⎢   ⋮                   ⋮              ⋱             ⋮                     ⋮          ⎥⎢ ⋮
                                                                                                                                                         ⎥⎢
                                                                                                                                                                         ⎥
                                                                                                                                                                         ⎥
                                                                   ⎢                                                                                     ⎥⎢ (N)          ⎥
                                                                   ⎣ M((N)0 )           M((N)1 )          ⋯        M((N)(N) )          M((N)(N+1) )      ⎦⎣ ü           ⎦
                                                                          M((N+1)0 )   M((N+1)1 )         ⋯       M((N+1)(N) )        M((N+1)(N+1) )          ü(N+1)
3. Governing equations for distorted domains                                                                    being w i (for each i = 0,..., n) a coefficient defining the influence of the
                                                                                                                corresponding control point P i on the i-th basis function Ni,p (u), con
    In the previous sections, the fundamental governing equations have
                                                                                                                sisting of a piecewise polynomial. The latter is a univariate B-Spline
been derived for a laminated doubly curved shell in a principal coor
                                                                                                                curve which has been defined setting a = 0 and b = 1. In particular, the
dinate parametrization setting. However, in cases of shells with arbi
                                                                                                                extended clamped knot vector Ω has been set so that Ω =
trary geometries, a key issue consists in finding a curvilinear principal                                       [                                               ]
coordinate system, according to differential geometry principles. For                                             a, …, a , up+1 , …, um− p− 1 , b, …, b , being p +1 the first and the last
                                                                                                                  ⏟̅̅̅̅̅⏞⏞̅̅̅̅̅⏟                 ⏟̅̅̅̅̅⏞⏞̅̅̅̅̅⏟
this reason, a mapping procedure of the physical domain is adopted                                                  p+1                                            p+1
(Fig. 1). Based on the domain distortion, it is possible to apply Eq. (53)                                      knot multiplicity and m the number of the selected breakpoints. To sum
also for any arbitrary shaped structure. To this end, the blending func                                        up, the well-known recursive procedure [22] has been adopted for the
tions α1 = α1 (ξ1 , ξ2 ) and α2 = α2 (ξ1 , ξ2 ) should be assessed first [22]:                                  derivation of Ni,p (u) due to its computationally-friendly nature. More
                                                                                                                specifically, it provides a piecewise polynomial, taking into account a
                           1(                                                                                              )
        α 1 (ξ1 , ξ2 ) =      (1 − ξ2 )ᾱ1(1) (ξ1 ) + (1 + ξ1 )ᾱ1(2) (ξ2 ) + (1 + ξ2 )ᾱ1(3) (ξ1 ) + (1 − ξ1 )ᾱ1(4) (ξ2 ) +
                           2
                                                                                                                                                                                                             (54)
               1(                                                                                                       )
           −      (1 − ξ1 )(1 − ξ2 )α1(1) + (1 + ξ1 )(1 − ξ2 )α1(2) + (1 + ξ1 )(1 + ξ2 )α1(3) + (1 − ξ1 )(1 + ξ2 )α1(4)
               4
                           1(                                                                                              )
        α 2 (ξ1 , ξ2 ) =      (1 − ξ2 )ᾱ2(1) (ξ1 ) + (1 + ξ1 )ᾱ2(2) (ξ2 ) + (1 + ξ2 )ᾱ2(3) (ξ1 ) + (1 − ξ1 )ᾱ2(4) (ξ2 ) +
                           2
                                                                                                                                                                                                             (55)
               1(                                                                                                       )
           −      (1 − ξ1 )(1 − ξ2 )α2(1) + (1 + ξ1 )(1 − ξ2 )α2(2) + (1 + ξ1 )(1 + ξ2 )α2(3) + (1 − ξ1 )(1 + ξ2 )α2(4)
               4
    As it can be seen, Eqs. (54)–(55) are based on a NURBS description of                                       step function for p = 0, and a linear combination of Ni,p− 1 (u) and
                                  (            )
the p = 1, ..., 4 external edge ᾱ1(p) , ᾱ2(p) of the structure within the                                     Ni+1,p− 1 (u) for an arbitrary p. This means that
                                    (            )
computational domain, whereas α1(q) , α2(q) for q = 1, ..., 4 represents                                                    {
                                                                                                                               1 if ui ⩽u < ui+1
the domain corners. In this way, the mapped shell geometry is                                                   Ni,0 (u) =
completely defined in a unified manner in a rectangular parent element.                                                        0 otherwise                                       (57)
When blending functions are assessed, a new set of natural coordinates                                                       u − ui                  ui+p+1 − u
                                                                                                                Ni,p (u) =            Ni,p− 1 (u) +               Ni+1,p− 1 (u)
ξ 1 ∈ [ − 1, +1] and ξ 2 ∈ [ − 1, +1] is introduced for the computational                                                   ui+p − ui               ui+p+1 − ui+1
domain mapping. The parameters ξ 1 and ξ 2 are defined from the
                                                                                                                    The mapping algorithm of the physical domain reported in Eqs. (54)–
normalization of the NURBS parametrization C(u) of the structural
                                                                                                                (55) is employed for the coordinate transformation α i (ξ1 , ξ2 ) for i = 1,2.
edges, being u ∈ [a, b] the general variable for the description of the
                                                                                                                In this way, the fundamental relations of Eq. (53) are suitable for arbi
curve, and a, b ∈ R. If we denote by p the degree of the curve, the
                                                                                                                trarily shaped domains. As a consequence, the first order derivatives of
function C(u) reads as follows [22]:
                                                                                                                the above-described blending functions are required. In a compact form,
         ∑n
               Ni,p (u)wi Pi                                                                                    they are accounted as follows [22]:
C(u) = ∑i=0  n                                                         (56)
             i=0 Ni,p (u)wi
                                                                                                        10
F. Tornabene et al.                                                                                                                                        Composite Structures 309 (2023) 116542
Fig. 2. Schematic representation of a general boundary distribution of linear elastic springs along the edges of a doubly-curved shell, for α i = α 1 , α 2 or ξ i = ξ 1 , ξ 2 .
A general surface variation is set along the three-dimensional thickness of the shell by properly setting a boundary distribution of elastic springs. A constant, linear
and parabolic distribution can be adopted along the thickness of the laminated structure.
⎡         ⎤   ⎡               ⎤⎡         ⎤
   ∂      ∂ξ1             ∂ξ2     ∂
⎢ ∂α ⎥ ⎢ ∂α                                                                                                ∂ξ1   1 ∂α2                                ∂ξ1      1 ∂α1
⎢ 1⎥ ⎢ 1                  ∂α1 ⎥⎢     ⎥
                              ⎥⎢ ∂ξ1 ⎥                                                           ξ1,α1 =       =          ,                ξ1,α2 =        =−
⎢     ⎥=⎢                     ⎥⎢     ⎥                                            (58)                     ∂α1 det(J) ∂ξ2                             ∂α2    det(J) ∂ξ2
⎣ ∂ ⎦ ⎣ ∂ξ1               ∂ξ2 ⎦⎣  ∂ ⎦                                                                                                                                                      (63)
                                                                                                         ∂ξ2      1 ∂α2                                   ∂ξ2   1 ∂α1
  ∂α2     ∂α2             ∂α2    ∂ξ2                                                           ξ2,α1 =       =−            ,                    ξ2,α2 =       =
                                                                                                         ∂α1    det(J) ∂ξ1                                ∂α2 det(J) ∂ξ1
   By employing the inverse transformation, a unified expression of the
                                                                                                   In this way, the multiplication by J− 1 of the first order derivatives
derivatives with respect to the natural coordinates ξ1 , ξ2 can be easily
                                                                                               with respect to the natural coordinates ξ1 , ξ2 leads to a direct compu
provided:
                                                                                               tation, within Eq. (53), of the derivatives with respect to α 1 , α 2 for an
⎡    ⎤ ⎡             ⎤⎡     ⎤    ⎡     ⎤
     ∂            ∂α1     ∂α2     ∂                ∂                                           arbitrarily shaped geometry.
⎢ ∂ξ ⎥ ⎢ ∂ξ
⎢ 1⎥ ⎢ 1                  ∂ξ1 ⎥⎢     ⎥
                              ⎥⎢ ∂α1 ⎥
                                          ⎢ ∂α ⎥
                                          ⎢ 1⎥                                                     Starting from Eq. (61), the second-order derivatives with respect to
     ⎥=⎢                             ⎥ = J⎢                                       (59)
⎢
⎣ ∂ ⎦ ⎣ ∂α1
                              ⎥⎢
                          ∂α2 ⎦⎣ ∂ ⎦      ⎣ ∂ ⎦
                                                ⎥                                              the principal coordinates α 1 , α 2 can also be performed in terms of ξ1 ,ξ2 ,
    ∂ξ2           ∂ξ2     ∂ξ2    ∂α2        ∂α2                                                as shown in Ref. [22]:
                                                                                          11
F. Tornabene et al.                                                                                                                                                                              Composite Structures 309 (2023) 116542
                          τ
u(nτ) = u(sτ) = u(ζ ) = 0              for τ = 0, 1, 2, ..., N, N + 1,              at ξ 1 = − 1 or ξ 1 = 1,     − 1⩽ξ 2 ⩽1
                          τ
                                                                                                                                                                                                                                                (69)
u(nτ)   =   u(sτ)   =   u(ζ )   =0     for τ = 0, 1, 2, ..., N, N + 1,              at ξ 2 = − 1 or ξ 2 = 1, − 1⩽ξ 1 ⩽1
   In the case of free edge (F), the following relations must be assessed                                        relationship presented in Ref. [22], in turn derived from the properties
on the interested edges:                                                                                         of the Lagrange polynomials L , valid for q⩾1:
                                          τ
Nn(τ) = 0,           (τ)
                    Nns  = 0, Tζ( ) = 0                 for τ = 0, 1, 2, ..., N, N + 1,        at ξ 1 = − 1 or ξ 1 = 1, − 1⩽ξ 2 ⩽1
                                                                                                                                                                                                                                                (70)
 (τ)
Nns     = 0,        Nn(τ)     = 0,    Tζ(τ)   =0        for τ = 0, 1, 2, ..., N, N + 1,        at ξ 2 = − 1 or ξ 2 = 1, − 1⩽ξ 1 ⩽1
                                                                                                            12
F. Tornabene et al.                                                                                                                                                Composite Structures 309 (2023) 116542
with k = 1, ..., I N I M according to the following relation:                                             points if compared to classical approaches. Based on the numerical
        (     )                                                                                           GDQ method of Eq. (73) for derivatives of a generic order, the pro
                             i = 1, ..., I N j = 1, ..., I M
A k = A ij            for                                                                   (76)          cedure allows a quadrature-based computation of integrals of a generic
                  k              k = i + (j − 1)I N
                                                                                                          thickness-dependent function f(ζ) according to the following expres
   Nevertheless, for the generic two-dimensional function f(ξ 1 , ξ 2 ), its                              sion [22]:
arbitrary (n + m)-th derivation order can be numerically assessed ac                                     ∫    ζk+1                 ∑
                                                                                                                                    IT
cording to the following matrix expression:                                                                             f (ζ)dζ =
                                                                                                                                       (                ) ( )
                                                                                                                                           w(k+1)g − wkg f ζ g                                     (81)
                                                                                                              ζk                    g=1
→(n+m)                  →
f      = ϛξ 1 (n)ξ 2 (m) f                                                                  (77)
                                                                                                          being I T the number of discrete points selected in the closed interval
being f and f (n+m) the matrices of the values assumed by the function f                                  [ζ k , ζ k+1 ], identified with ζ g for g = 1, ..., I T . In the present work, each
and its (n + m)-th order derivative, respectively, at each point of the                                   k-th layer of the stacking sequence of an l laminate has been discretized
computational domain, vectorized according to Eq. (75). Starting from                                     by means of the CGL distribution, as follows [22]:
Eq. (73), the GDQ weighting coefficients of Eq. (74) for the derivatives                                                    (        )
                                                                                                                              g− 1
along ξ 1 and ξ 2 are collected in a I N I M × I N I M matrix ϛξ 1 (n)ξ 2 (m) ,                           ξkg = − cos               π , k = 1, ..., l, g = 1, ..., IT for ξkg ∈ [ − 1, 1]
                                                                                                                             IT − 1
developed employing the well-known Kronecker product [22], as
                                                                                                                                                                                         (82)
follows:
                                                                                                              The discretization employed in Eq. (82), defined in the closed in
ϛξ 1 (n)ξ 2 (m) = ϛξ 2 (m) ⊗ ϛξ 1 (n)                                                       (78)
                                                                                                          terval [ − 1, 1], can be shifted to the physical interval [ζ k , ζ k+1 ] via the
                                                                                                          following coordinate transformation [22]:
being ϛξ r (q) for r = 1, 2 and q = n, m a squared matrix of order I Q = I N ,
I M . As also observed in the previous sections, when n = 0 or m = 0 the                                  ζkg =
                                                                                                                ζk+1 − ζk (        )
                                                                                                                            ξkg + 1 + ζk , k = 1, ..., l,              g = 1, ..., IT   for
identity matrix should be introduced, setting ϛξ r (0) = I.                                                         2                                                                              (83)
     When the discrete form of the fundamental governing equation for a                                   ξkg ∈ [ − 1, 1]
doubly-curved shell is considered according to Eq. (53), for each τ, η = 0,                                   Since the GIQ method reported in Eq. (81) comes from the general
..., N + 1, partial derivatives with respect to α 1 and α 2 principal di                                 properties of the fundamental theorem of integrals, the weighting co
rections can be performed starting from the derivation with respect to ξ 1                                efficients wkg and w(k+1)g are computed according to the recursive
and ξ 2 , according to Eq. (77). This means that matrices ϛα1 (1)α2 (0) and                               formulation set in Eq. (74). The interested reader can find more details
ϛα1 (0)α2 (1) must be introduced for n, m = 0, 1, as follows:                                             in Ref. [22]. In the present article, a free vibration analysis of doubly-
                → ◦                  → ◦                                                                  curved shells of arbitrary geometries is performed. Based on Eq. (51),
ϛα1 (1)α2 (0) = P 11 ϛξ1 (1)ξ2 (0) + P 21 ϛξ1 (0)ξ2 (1)
                →   ◦                → ◦                                                    (79)          for each η-th order of the kinematic expansion, we define the generalized
ϛα1 (0)α2 (1) = P 12 ϛξ1 (1)ξ2 (0) + P 22 ϛξ1 (0)ξ2 (1)
                                                                                                          displacement field vector u(η) (α 1 , α 2 , t) by applying a separation vari
                                                                                                          able method among the time-dependent and space-dependent quantities
where Pri for r = 1, 2 and i = 1, 2, collects the coefficients of the first
                                                                                                          [22]:
order derivatives of the blending transformations of Eq. (54)–(55), ac
cording to Eq. (63), and the symbol ◦ denotes the well-known Hadamard                                     u(η) (α 1 , α 2 , t) = U(η) (α 1 , α 2 )eiωt                                             (84)
product [22]. As far as the second order derivatives with respect to α 1
and α 2 are concerned, the weighting coefficient matrices are calculated                                  being ω the circular frequency and U(η) (α 1 , α 2 ) the corresponding modal
starting from Eqs. (64)–(66), providing the following equations:                                          shape vector, for each η = 0, ..., N + 1. Employing Eq. (84) in the
                →o2                   →o2                    → →                           →                      →
ϛα1 (2)α2 (0) = P 11 ◦ϛξ1 (2)ξ2 (0) + P 21 ◦ϛξ1 (0)ξ2 (2) + 2 P 11 ◦ P 21 ◦ϛξ1 (1)ξ2 (1) + P 111 ◦ϛξ1 (1)ξ2 (0) + P 211 ◦ϛξ1 (0)ξ2 (1)
                →o2                   →o2                    → →                           →                      →
ϛα1 (0)α2 (2) = P 12 ◦ϛξ1 (2)ξ2 (0) + P 22 ◦ϛξ1 (0)ξ2 (2) + 2 P 12 ◦ P 22 ◦ϛξ1 (1)ξ2 (1) + P 122 ◦ϛξ1 (1)ξ2 (0) + P 222 ◦ϛξ1 (0)ξ2 (1)                                                             (80)
                → →                          → →                          (→ →              → →          )                →                 →
ϛα1 (1)α2 (1) = P 11 ◦ P 12 ◦ϛξ1 (2)ξ2 (0) + P 21 ◦ P 22 ◦ϛξ1 (0)ξ2 (2) + P 11 ◦ P 22 + P 12 ◦ P 21 ◦ϛξ1 (1)ξ2 (1) + P 112 ◦ϛξ1 (1)ξ2 (0) + P 212 ◦ϛξ1 (0)ξ2 (1)
   As it can be seen, the coordinate transformation has been performed                                    fundamental relations (51), the free vibration governing equations for
thanks to the computation of the matrices Prij , with r = 1, 2 and i,j = 1,                               each τ = 0, ..., N +1 assume the following form:
2. They collect all the second order derivatives ξr,αi ,αj of the natural
                                                                                                           ∑
                                                                                                           N +1
                                                                                                                (                      )          ∑
                                                                                                                                                  N +1
coordinates ξr with respect to αi and αj , which have been performed                                                    L(τη) − Lf (τη) U(η) + ω2      M(τη) U(η) = 0        for
following the Eqs. (64)–(66). The numerical integration of the gener                                      η=0                                     η=0                                             (85)
alized stiffness coefficients A        , for τ, η = 0, ..., N +1 and α i , α j =
                                          (τη)αi αj                                                       τ = 0, 1, 2, ..., N, N + 1
α 1 , α 2 , α 3 , follows the Generalized Integral Quadrature (GIQ) algo                                     It is useful to rearrange the DOFs of the discretized system from a
rithm. The use of such a numerical technique for the computation of                                       separation within the computational grid defined in Eq. (72), of DOFs
integrals lies in the fact that Eq. (23) accounts for an integration on a                                 denoted with δb related to the boundaries (“b” nodes) from the other
variable interval. Furthermore, a general expression is considered here                                   unknown variables δd referred to the inner points (“d” points) [22]:
for the thickness functions. Moreover, the GIQ method allows one to                                       [                         ][ ]      [       ][ ]
compute integrals also for the case of a general through-the-thickness                                               f           f
                                                                                                            L̄bb − L̄bb L̄bd − L̄bd   δb     2 0   0     δb
variation of the material properties within each layer, as happens for                                                                   =  ω                                 (86)
                                                                                                            L̄ − L̄
                                                                                                                     f
                                                                                                                         L̄ − L̄
                                                                                                                                 f    δd        0 M̄dd δd
example when FGM laminae are considered in the lamination scheme.
                                                                                                                   db       db      dd      dd
Furthermore, the GIQ has been demonstrated to converge more steadily                                         The final form of the discrete stiffness matrix of Eq.(86), as well as
to stable and accurate solutions with a reduced number of integration                                     the mass one, follows the DOFs rearrangement at issue. In this way, it is
                                                                                                     13
F. Tornabene et al.                                                                                                               Composite Structures 309 (2023) 116542
Fig. 3. Geometric input data for the dynamic analysis of a rectangular plate and a catenoid mapped with a domain subjected to general external boundary conditions
and characterized by variable thickness along the principal lines. The equations of each reference surface have been written in curvilinear principal coordinates [22].
The computational domain mapping has been performed starting from the NURBS description of the structural edges.
                                                                                  14
F. Tornabene et al.                                                                                                          Composite Structures 309 (2023) 116542
Fig. 4. Geometric input data for the dynamic analysis of an elliptic cone and a revolution paraboloid mapped with a general domain subjected to general external
boundary conditions and characterized by variable thickness along the principal lines. The equations of each reference surface have been written in curvilinear
principal coordinates [22]. The computational domain mapping has been performed starting from the NURBS description of the structural edges.
                                                                               15
F. Tornabene et al.                                                                                                          Composite Structures 309 (2023) 116542
Fig. 5. Geometric input data for the dynamic analysis of a hyperbolic paraboloid and a helicoid mapped with a domain subjected to general external boundary
conditions and characterized by variable thickness along the principal lines. The equations of each reference surface have been written in curvilinear principal
coordinates [22]. The computational domain mapping has been performed starting from the NURBS description of the structural edges.
                                                                               16
F. Tornabene et al.                                                                                                          Composite Structures 309 (2023) 116542
Fig. 6. Geometric input data for the dynamic analysis of a pseudosphere and a hyperbolic hyperboloid mapped with a domain subjected to general external boundary
conditions and characterized by variable thickness along the principal lines. The equations of each reference surface have been written in curvilinear principal
coordinates [22]. The computational domain mapping has been performed starting from the NURBS description of the structural edges.
                                                                               17
F. Tornabene et al.                                                                                                                                                                         Composite Structures 309 (2023) 116542
Table 1
First ten mode frequencies of a rectangular plate [22] with variable thickness and mapped geometry. It has been externally constrained with a three-dimensional set of
linear boundary springs characterized by different distributions along the external edges. The analysis refers to different higher order theories. The results have been
compared to the outcomes of a refined three-dimensional finite element model.
                             (             )
  Mapped Rectangular Plate BKDDD FBKDDD F
  Double - Weibull distribution – sensitivity analysis
possible to reduce the dimension of the discrete system of Eq. (86) taking                                          characterized by different syngonies within each layer. More specif
into account only the DOFs related to the inner part δd of the compu                                               ically, anisotropic, orthotropic and isotropic continua have been
tational domain, since the remaining ones δb are already subjected to the                                           considered. To the first class belongs the following triclinic material
                                                                                                                    ( (k)               )
external constraints. From the application of such DOFs condensation,                                                ρ = 7750 kg/m3 , whose stiffness matrix has been taken from the
one gets:                                                                                                           DiQuMASPAB database [105]. It is expressed according to the conven
(      ((
                 f )  (        f )(        f )− 1 (          f )
                                                                 )       )                                          tion adopted in Eq. (19) with respect to the k-th layer material reference
 ¯M−dd1 L̄dd − L̄dd − L̄db − L̄db L̄bb − L̄bb       L̄bd − L̄bd    − ω2 I δd
                                                                                                                    system O α 1 α 2 ζ(k) :
                                                                                                                                 ′   (k) (k)
  =0                                                                                                                      ⎡                                                                   ⎤
                                                                                                      (87)                      (k)
                                                                                                                             ⎢ C11       C12
                                                                                                                                            (k)       (k)
                                                                                                                                                    C16
                                                                                                                                                                 (k)
                                                                                                                                                               C14
                                                                                                                                                                             (k)
                                                                                                                                                                            C15
                                                                                                                                                                                        (k)
                                                                                                                                                                                       C13 ⎥
                                                                                                                             ⎢                                                              ⎥
                                                                                                                             ⎢ (k)                                                      (k) ⎥
being I the identity matrix.
                                                                                                                                            (k)       (k)        (k)         (k)
                                                                                                                             ⎢ C12       C22        C26        C24          C25        C23 ⎥
                                                                                                                             ⎢ (k)                                                          ⎥
                                                                                                                             ⎢C          C26
                                                                                                                                            (k)       (k)
                                                                                                                                                    C66
                                                                                                                                                                 (k)
                                                                                                                                                               C46
                                                                                                                                                                             (k)
                                                                                                                                                                            C56        C36 ⎥
                                                                                                                                                                                        (k)
                                                                                                                             ⎢                                                              ⎥
5. Applications and results                                                                                         C(k)   = ⎢ 16                                                           ⎥
                                                                                                                             ⎢ C(k)       (k)
                                                                                                                                         C24
                                                                                                                                                     (k)
                                                                                                                                                    C46
                                                                                                                                                                (k)
                                                                                                                                                               C44
                                                                                                                                                                             (k)
                                                                                                                                                                            C45        C34 ⎥
                                                                                                                                                                                        (k)
                                                                                                                             ⎢ 14                                                           ⎥
                                                                                                                             ⎢ (k)        (k)        (k)        (k)          (k)        (k) ⎥
                                                                                                                             ⎢ C15       C25        C56        C45          C55        C35 ⎥
    In the present section we analyze the modal response of some                                                             ⎢                                                              ⎥
                                                                                                                             ⎣ (k)                                                      (k) ⎦
laminated panels, involving zero-curved, single-curved and doubly-
                                                                                                                                          (k)        (k)        (k)          (k)
                                                                                                                               C13       C23        C36        C34          C35        C33
curved structures, characterized by an arbitrary shape, completely                                                           ⎡                                       ⎤
anisotropic layers and a general thickness variation along the principal                                                      98.84 53.92  0.03   1.05  − 0.1 50.78
                                                                                                                            ⎢ 53.92 99.19         0.55 − 0.18 50.87 ⎥
parametric lines. The physical domain is described, for each case,                                                          ⎢              0.03                      ⎥
                                                                                                                            ⎢ 0.03                             0.02 ⎥
assuming a parameterization in principal coordinates. In addition, the                                                      ⎢
                                                                                                                           =⎢
                                                                                                                                     0.03  22.55 − 0.04 0.25         ⎥ GPa                                                          (88)
                                                                                                                            ⎢ 1.05   0.55 − 0.04 21.1    0.07  1.03 ⎥
use of NURBS-based blending functions, as discussed in Eqs. (54)–(55),                                                      ⎣ − 0.1 − 0.18 0.25
                                                                                                                                                                     ⎥
                                                                                                                                                  0.07  21.14 − 0.18 ⎦
allows one to solve the free vibration linear discrete system reported in
                                                                                                                              50.78 50.87  0.02   1.03 − 0.18 87.23
Eq. (87) referred to a squared computational domain (Fig. 1). For each
case study, we describe the mapped domain, together with the                                                              The three-dimensional stiffness matrix of the trigonal material [105]
                                                                                                                    (                    )
geometrical description of the structures and thickness variation, as                                                   ρ(k) = 2649 kg/m3 is accounted, instead, as follows:
shown in Figs. 3-6. The lamination schemes involve materials
                                                                                                               18
F. Tornabene et al.                                                                                                                                                                                 Composite Structures 309 (2023) 116542
Table 2
Mode frequencies of a rectangular plate [22] with variable thickness and mapped geometry composed by three layers of general-oriented anisotropic and orthotropic
materials. The external West (W) and East (E) edges have been constrained with a super elliptic distribution of linear springs with a pre-determined three-dimensional
translational stiffness in each principal direction. The influence of the extension of the spring distribution on the modal response has been pointed out.
                             (             )
  Mapped Rectangular Plate BKDDD FBKDDD F
  Double - Weibull distribution – sensitivity analysis
            ⎡                                                       ⎤
                                                                                                                                  (k)
               (k)
            ⎢ C11          (k)
                          C12
                                    (k)
                                   C16
                                            (k)
                                           C14
                                                   (k)
                                                  C15
                                                          (k)
                                                         C13        ⎥                                                           E 1 = 137.90 GPa                 G(k)
                                                                                                                                                                   12 = 7.10 GPa                    ν(k)
                                                                                                                                                                                                      12 = 0.30
            ⎢                                                       ⎥
            ⎢ (k)          (k)      (k)     (k)    (k)    (k)       ⎥                                                             E(k)   = 8.96 GPa
                                                                                                                                                                  (k)
                                                                                                                                                                 G 13      = 7.10 GPa               ν(k)
                                                                                                                                                                                                      13 = 0.30
                                                                                                                                                                                                                                            (91)
            ⎢ C12         C22      C26     C24    C25    C23        ⎥                                                               2
            ⎢ (k)                                                   ⎥
            ⎢C             (k)
                          C26       (k)
                                   C66      (k)
                                           C46     (k)
                                                  C56     (k)
                                                         C36        ⎥                                                             E(k)
                                                                                                                                    3    = 8.96 GPa               (k)
                                                                                                                                                                 G 23      = 6.21 GPa               ν(k)
                                                                                                                                                                                                      23   = 0.49
            ⎢ 16                                                    ⎥
C(k) =      ⎢ (k)                                                   ⎥
            ⎢C             (k)
                          C24       (k)
                                   C46      (k)
                                           C44     (k)
                                                  C45     (k)
                                                         C34        ⎥
            ⎢ 14
            ⎢ (k)
                                                                    ⎥
                                                                    ⎥                                                           whereas ρ(k) = 1450 kg/m3 is the associated material density.
                           (k)      (k)     (k)    (k)    (k)
            ⎢ C15
            ⎢             C25      C56     C45    C55    C35        ⎥
                                                                    ⎥                                                               In some cases, an isotropic central core has been included within the
            ⎣ (k)                                                   ⎦                                                                                                                       ( (k)
                                                                                                                                lamination scheme. In particular, zirconia                   E = 168.00
                           (k)      (k)     (k)    (k)    (k)
              C13         C23      C36     C34    C35    C33
                                                                                                                                                                       )             (
            ⎡                                                                               ⎤                                   GPa, ν(k) = 0.3, ρ(k) = 5700 kg/m3 , steel             E(k) = 210.00 GPa,
            86.74                6.99        0      0                   − 17.91      11.91                                                                       )                                  ( (k)
         ⎢ 6.99                  86.74       0      0                    17.91       11.91 ⎥                                     ν(k) = 0.3, ρ(k) = 7800 kg/m3           and       aluminum          E =
         ⎢                                                                                  ⎥                                                                                 )
         ⎢ 0                       0       39.88 − 17.91                   0           0 ⎥                                      70.00 GPa, ν(k) = 0.3, ρ(k) = 2707 kg/m3 materials [105] have been
        =⎢                                                                                  ⎥ GPa
         ⎢ 0
         ⎢                         0      − 17.91 57.94                    0           0 ⎥  ⎥                                   adopted in different combinations.
         ⎣ − 17.91               17.91       0      0                    57.94         0 ⎦                                          The external boundary conditions are assessed employing the
            11.91                11.91       0      0                      0         107.20                                     generalized formulation derived in Eqs. (40)–(41) for each side of the
                                                                                                                  (89)          mapped shell object of analysis (Fig. 2). Besides, they are obtained from
   The orthotropic materials considered in the present work refer to the                                                        a correct definition of the governing parameters of the linear springs
graphite-epoxy and glass–epoxy composites, whose linear elastic                                                                 distributed along the structural edges, taking into account both in-plane
behavior is described by means of the well-known engineering con                                                               and out-of-plane external constraints thickness variation, as reported in
                                                                                                                                Eqs. (40)–(41). Moreover, a second set of linear springs have been set
stants, referred to the material reference system O α 1 α 2 ζ(k) , as follows
                                                                                ′   (k) (k)
                                                                                                                                along the shell top and bottom surfaces, following the Winkler gener
[105]:
                                                                                                                                alized model introduced in Eq. (27).
E(k)                               G(k)                      ν(k)                                                                   A validation analysis starts considering the vibration response of a
  1 = 53.78 GPa                      12 = 8.96 GPa             12 = 0.25
                                                                                                                                mapped rectangular plate, as provided by the present structural model
  (k)
E 2 = 17.93 GPa                    G(k)
                                     13 = 8.96 GPa           ν(k)
                                                               13 = 0.25
                                                                                                                  (90)
                                                                                                                                and a three-dimensional FEM simulation. The natural frequencies are
E(k)
  3     = 17.93 GPa                G(k)
                                     23   = 3.45 GPa         ν(k)
                                                               23   = 0.34                                                      summarized comparatively in Tables 1-2, and the related mode shapes
    The density of the material of interest has been assumed equal to                                                           are represented in Figs. 7-9. Moreover, a parametric investigation is
ρ(k) = 1900 kg/m3 . The orthotropic behavior with respect to                                                                    performed to check for the influence of the lateral constraints shape
                                                                                                                                parameter, as well as the stiffness of the selected boundary springs, as
O α(k)  (k) (k)
                of the graphite-epoxy layer [105] considers the following
  ′
     1 α2 ζ
                                                                                                                                summarized in Table 3. Then, we investigate the influence on the dy
elastic constants:
                                                                                                                                namic response of curved structures of the higher order assumption of
                                                                                                                                the field variable according to Eq. (10). In particular, the assumed field
                                                                                                                           19
F. Tornabene et al.                                                                                                                               Composite Structures 309 (2023) 116542
Fig. 7. Mode shapes of a mapped laminated rectangular plate with variable thickness and general boundary conditions. A linear thickness variation has been
selected. The external constraint has been obtained from a Double - Weibull distribution of translational springs distributed along the three principal directions so
that the four corners of the structure are bounded. Corresponding mode frequencies, as well as other analysis information, have been collected in Table 1.
tion (Z) of Eq. (11) is adopted too. When the shear correction factor κ(ζ)              U 2 (α 1 , α 2 , ζ, t) = u(0)                    (1)
                                                                                                                   2 (α 1 , α 2 , t) + ζu 2 (α 1 , α 2 , t)
                                                                                                                                                                                  (92)
is included in the model according to Eq. (24), it is declared as a su                                                       (0)
                                                                                                   U 3 (α 1 , α 2 , ζ, t) = u 3 (α 1 , α 2 , t)
perscript. Moreover, in cases of reduced elastic coefficients Eij = Qij
                                                                     (k)     (k)
                                                                                            Besides, when the TSDT model is employed, a third order expansion
within each stiffness matrix for each k-th layer according to Eq. (20), the
                                                                                        of the in-plane components of the field variable is considered as follows
symbol RS is introduced as a subscript. Some simulations rely on clas
                                                                                        [22]:
sical ESL theories. For the sake of completeness the FSDT displacement
field is defined as [22]:
                                                                                   20
F. Tornabene et al.                                                                                                                Composite Structures 309 (2023) 116542
Fig. 8. Mode shapes of a mapped laminated rectangular plate with variable thickness and general boundary conditions. A linear thickness variation has been
selected. The external constraint has been obtained from a super elliptic distribution of translational springs distributed along the three principal directions so that
the midpoints belonging to the four edges of the structure are bounded. Corresponding mode frequencies, as well as other analysis information, have been collected
in Table 1.
                                                                                  21
F. Tornabene et al.                                                                                                                Composite Structures 309 (2023) 116542
Fig. 9. Mode shapes of a mapped laminated rectangular plate with variable thickness and general boundary conditions. A linear thickness variation has been
selected. The external constraint has been obtained from a super elliptic distribution of translational springs distributed along the three principal directions applied
on the West (W) and East (E) sides of the structure. Corresponding mode frequencies, as well as other analysis information, have been collected in Table 1.
                                                                                  22
F. Tornabene et al.                                                                                                                                                           Composite Structures 309 (2023) 116542
     For an easy identification of the general external constraints within                                                         properly setting the parameters of the general thickness expression re
each structural side, we adopt the notation BKTTT , which means that the                                                           ported in Eq. (4). As it can be seen in Table 1, a lamination scheme made
side at issue is bounded (B) with a set of linear springs (K) oriented along                                                       of three generally oriented layers has been employed, characterized by a
α 1 , α 2 and ζ principal directions (Figs. 1-2). The symbol T stands for the                                                      central triclinic core covered by two external orthotropic graphite-epoxy
linear spring distribution along the edge, for each direction. More spe                                                           skins. Four different boundary conditions have been modelled, and the
cifically, T = D denotes the Double – Weibull distribution introduced in                                                           first ten natural frequencies have been provided for each case,
Eq. (36), whereas the super elliptic distribution of Eq. (37) is identified                                                        employing different higher order theories. The fully-clamped (CCCC)
with T = S.                                                                                                                        boundary condition has been obtained according to Eq. (69). Among all
     Natural frequencies of doubly-curved shells subjected to general                                                              the ESL higher order theories considered in the numerical investigation,
boundary conditions have been reported in Tables 4-10, whereas the                                                                 the most accurate results with respect to the 3D FEM outcomes are
mode shapes are collected in Figs. 10-22. The formulation presented in                                                             provided by the EDZ4 theory. This is related to the complex nature of the
this work provides very accurate results with respect to the FEM solu                                                             lamination scheme, with coupling issues between adjacent laminae, and
tions, despite the lower DOFs required by the GDQ approach. In addi                                                               complex in-plane and out-of-plane deformation effects within each
tion, the differential geometry description of shells allows one to easily                                                         layer. A fully-constrained condition at the corners of the structure in the
address very complex structures, as well as the higher order assumption                                                            mapped geometry has been modelled taking into account, for each side
                                                                                                                                                                                        (   )
is capable of considering several warping and coupling effects that                                                                of the structure, a Double – Weibull distribution BKDDD of linear springs
cannot be depicted by classical ESL approaches. In this way, a complete                                                            with reference stiffness equal to 1⋅1021 N/m3 in each principal direction,
three-dimensional capability can be achieved even in the case of a two-                                                            setting also ξ̄ m = ̃
                                                                                                                                                       ξ m = 0.0025 and p = 20. As it can be seen in Table 1,
dimensional model.                                                                                                                 the choice of a higher order model with N = 4 is essential to capture
                                                                                                                                   warping effects along the thickness, as well as the zigzag kinematic
                                                                                                                                   assumption of Eq. (11) for the N + 1-th expansion order. Furthermore,
5.1. Modal analysis of a mapped plate
                                                                                                                                   Fig. 7 shows that the symmetry of the boundary conditions does not
                                                                                                                                   induce symmetric modes due to the anisotropy of the layers in the
    We now assess the accuracy of the higher order theory presented in
                                                                                                                                   stacking sequence. Moreover, such very complex external constraints
this work, based on the free vibration analysis performed on a mapped
                                                                                                                                   induce, in higher modes, some localized deformations which are
rectangular plate. The main geometric features of the structure have
                                                                                                                                   completely disregarded by lower order ESL theories. Finally, the GDQ
been reported in Fig. 3, referred to the Domain-01. In both α 1 and α 2
                                                                                                                                   method has proved to be an accurate numerical capable to interpret such
principal directions, a linear thickness variation has been selected by
Table 3
Parametric analysis on first ten mode frequencies of a laminated rectangular plate [22] with variable thickness externally constrained by a set of boundary linear
translational springs applied on the East (E) and West (W) edges of the structure, whereas North (N) and South (S) sides of the shell are fully clamped. Reference
solutions have been calculated with both the ESL GDQ approach and the finite element method. As can be seen, natural frequencies decrease, for each mode number, as
the springs’ thickness decreases.
                            (            )
  Mapped Rectangular Plate CBKCCC CBKCCC
  Linear Boundary Springs – sensitivity analysis
1 2 3 4 5 6 7 8 9 10 DOFs
  3D FEM           0⋅10 0             189.811            309.828                   418.143          472.516              561.695                672.275   714.288   732.047   885.974       893.054      2131302
    (CCCC)
  CCCC             0⋅100              191.448            312.641                   424.668          477.000              570.809                679.497   729.446   743.686   904.394       905.330      16182
  CBCB             1⋅1021             191.448            312.641                   424.668          477.000              570.809                679.497   729.446   743.686   904.394       905.330      16182
  CBCB             1⋅1020             191.448            312.641                   424.668          477.000              570.809                679.497   729.446   743.686   904.394       905.330      16182
  CBCB             1⋅1015             191.461            312.661                   424.697          477.025              570.850                679.521   729.491   743.738   904.416       905.392      16182
  CBCB             1⋅1013             194.645            317.039                   431.657          482.025              580.121                684.848   737.350   759.753   908.817       918.267      16182
  CBCB             1⋅1012             172.031            172.596                   262.376          358.305              447.225                522.679   571.898   618.313   663.252       784.893      16182
  CBCB             1⋅1011             125.175            239.683                   246.290          300.771              352.251                414.991   467.260   573.874   594.966       658.730      16182
  CBCB             1⋅1010             119.689            178.693                   246.778          282.975              325.546                326.378   387.837   389.226   389.282       443.478      16182
  CBCB             1⋅105               63.861             76.178                   166.838          195.393              219.081                310.877   344.791   385.306   402.599       455.248      16182
  CBCB             1⋅100               63.864             76.175                   166.840          195.392              219.080                310.878   344.792   385.306   402.598       455.248      16182
  CFCF             0⋅100               63.864             76.175                   166.840          195.392              219.080                310.878   344.792   385.306   402.598       455.248      16182
  3D FEM           0⋅100               63.640             75.555                   166.031          194.315              217.638                309.013   342.572   384.917   399.353       451.292      2131302
    (CFCF)
                                                       (k)ξ01       (k)ξ11               (k)ξ01         (k)ξ11            (k)ξ01       (k)ξ11
  Boundary Springs: Constant distribution, k 1f                 = k 1f        = k f , k 2f        = k 2f         = k f , k 3f      = k 3f       = kf
  Lamination Scheme: (0/30/45),            h01
                                            =    h03
                                                   = 0.02 m,                 h02
                                                                        = 0.03 m
    Laminae Material Sequence: 1 graphite-epoxy, 2nd triclinic, 3rd graphite-epoxy
                                          st
  ESL Analysis: CGL computational grid with IN = 31 and IM = 33 discrete points, EDZ4 displacement field
  Thickness Variation: ϕ 1 (α 1 , α 2 ) =¯αp11 , ϕ 2 (α 1 , α 2 ) =¯αp21 , p 1 = p 2 = 1, δ 1 = δ 2 = 1.0
    Isogeometric Mapping: Domain-01
                                                                                                                             23
F. Tornabene et al.                                                                                                                                                                            Composite Structures 309 (2023) 116542
Table 4
Mode frequencies of a mapped catenoid [22] characterized by variable thickness and subjected to general boundary linear elastic constraints. The results have been
provided by the unified ESL approach employing different higher order theories. Elastic translational springs have been distributed along the structure in all the
principal directions in order to obtain different kinds of elastic supports. A sinusoidal thickness variation has been set along the physical domain.
  Mapped Catenoid (a = 0.85 m)
  Mode          FSDTκ=1.2
                    RS
                                         FSDTZ RS              TSDT RS             TSDTZ RS     ED1 RS            EDZ1 RS        ED2κ=1.2           EDZ2           ED3       EDZ3        ED4              EDZ4
  f [Hz]
  DOFs    5046        7569           10092                                         12615        5046              7569           7569               10092          10092     12615       12615            15138
                        (              )
  Clamped Opposite Edges FBKSSS FBKSSS
  1       79.650      79.280         78.155                                        78.530       78.682            79.518         78.912             77.713         78.732    79.054      78.531           79.135
  2       127.752     127.008        125.098                                       124.904      127.014           128.800        126.395            125.033        126.059   125.862     125.475          125.355
  3       157.365     156.683        154.962                                       155.160      156.450           158.050        156.876            154.514        156.428   156.594     155.846          156.055
  4       263.850     262.853        260.755                                       260.386      259.129           265.925        263.988            261.032        263.115   262.501     262.314          261.642
  5       297.436     296.274        291.856                                       290.822      294.327           297.688        294.612            290.528        293.975   292.703     292.989          291.523
  6       339.774     338.280        332.788                                       331.563      337.601           340.799        335.981            331.474        335.152   333.833     333.733          332.130
  7       454.421     452.713        448.363                                       447.777      445.382           458.147        453.663            448.677        452.069   451.044     450.512          449.519
  8       532.191     529.955        522.704                                       522.600      526.499           530.117        525.672            517.861        525.230   524.931     523.928          523.584
  9       538.690     534.236        526.735                                       525.436      533.769           534.916        531.932            522.755        529.638   528.354     527.963          526.251
  10      638.399     635.502        629.633                                       628.565      628.869           639.028        636.517            627.367        633.980   632.482     631.816          630.244
  Boundary Springs:
    North (N) edge: super elliptic distribution (ξ̄m = 0, ̃
                                                          ξm = 0.23, p = 1000),
    South (S) edge: super elliptic distribution (ξ̄m = 1, ̃
                                                          ξm = 0.23, p = 1000),
     (k)ξ01        (k)ξ11                                    (k)ξ0        (k)ξ11                       (k)ξ01       (k)ξ11
    k 1f        = 1⋅1021 N/m3 , k 2f 1 = k 2f
              = k 1f                                                               = 1⋅1021 N/m3 , k 3f         = k 3f       = 1⋅1021 N/m3
                      (                        )
  Clamped Half Corners BKSSS BKSSS BKSSS BKSSS
  1       164.778      163.454          160.020                                    158.397      164.520           162.408        160.971            158.084        160.051   158.552     159.590          157.391
  2       202.815      201.954          199.338                                    198.970      200.193           202.595        201.580            198.448        200.741   200.227     200.151          199.564
  3       260.273      259.581          257.358                                    259.958      257.071           258.071        257.894            253.486        258.880   261.287     258.321          262.303
  4       330.007      328.372          325.096                                    325.875      326.995           330.995        327.801            323.192        327.712   328.291     326.479          327.929
  5       353.012      353.718          350.375                                    351.048      350.334           353.467        352.351            346.776        353.292   353.803     351.908          352.907
  6       405.357      403.741          399.330                                    397.378      404.219           410.293        403.374            398.637        403.101   401.303     400.845          398.325
  7       533.028      530.502          522.581                                    521.224      527.680           528.418        525.524            516.959        524.318   522.738     523.060          521.100
  8       570.632      568.689          565.682                                    564.924      553.821           575.986        576.232            569.951        573.826   571.881     572.961          571.824
  9       586.363      584.852          577.487                                    577.198      581.094           587.051        582.986            573.805        582.466   581.942     580.571          580.191
  10      633.816      630.781          622.343                                    621.141      622.975           631.284        628.401            618.961        625.670   623.977     624.335          622.535
                                                                                                                             0        1         0              1                     0         1          0           1
  Boundary Springs: super elliptic distribution (ξ̄m = 0.5, ̃ξm = 0.01, p = 40), k 1f 1 = k 1f 1 = k 1f 2 = k 1f 2 = 1⋅1021 N/m3 , k 2f 1 = k 2f 1 = k 2f 2 = k 2f 2 = 1⋅1021 N/m3 ,
                                                                                  (k)ξ     (k)ξ     (k)ξ     (k)ξ                   (k)ξ     (k)ξ     (k)ξ     (k)ξ
  1            203.957                   203.065               201.325             201.045      196.960           203.547        203.637            200.980        202.553   201.904     202.374          201.714
  2            221.594                   220.819               219.334             219.066      216.283           224.309        222.379            220.082        221.531   220.955     220.781          220.196
  3            258.021                   256.607               252.585             252.059      255.560           256.056        254.583            250.428        253.691   253.046     253.134          252.402
  4            287.970                   286.107               281.520             280.957      285.309           286.594        284.164            279.613        283.115   282.482     282.306          281.510
  5            376.633                   375.318               372.418             371.854      365.134           376.304        376.197            371.577        374.630   373.478     374.080          372.897
  6            463.161                   461.607               458.101             457.416      448.495           465.191        464.340            458.486        462.046   460.604     461.065          459.721
  7            496.165                   493.589               487.655             486.840      487.056           493.750        491.749            484.890        489.880   488.603     489.046          487.702
  8            565.496                   563.577               556.512             555.617      559.286           564.546        560.276            551.932        558.877   557.785     557.543          556.251
  9            568.825                   565.720               556.730             555.740      562.673           566.093        561.631            553.701        559.801   558.574     558.056          556.756
  10           572.058                   569.720               564.830             563.867      566.090           579.792        572.050            565.801        570.817   569.542     567.801          566.507
                                                                                                                                  0         1              0         1                     0          1           0       1
  Boundary Springs: Double - Weibull distribution (ξ̄m = 0.05, ̃ξm = 0.05, p = 100), k 1f 1 = k 1f 1 = k 1f 2 = k 1f 2 = 1⋅1021 N/m3 , k 2f 1 = k 2f 1 = k 2f 2 = k 2f 2 = 1⋅
                                                                                      (k)ξ     (k)ξ     (k)ξ     (k)ξ                   (k)ξ     (k)ξ     (k)ξ     (k)ξ
phenomena with a reduced number of discrete points.                                                                                   dimensional model provides bending modes coupled with out-of-plane
    It must be remarked that the choice of different boundary conditions                                                              stretching effects. Such a behavior can be depicted only by the EDZ4
does not affect the accuracy of the solution provided by the proposed                                                                 kinematic assumption. As also shown in Fig. 9, mode shapes exhibit a
                                                        (                      )
model, as it can be seen from the second investigation BKSSS BKSSS BKSSS BKSSS                                                        complex warping phenomenon along the thickness of the plate in higher
in which the structure has been fixed at the middle of each edge. For                                                                 modes. Lower modes show a non-standard bending direction, because of
each principal direction the external constraints consider a super elliptic                                                           the complete anisotropy of the lamination scheme.
distribution (S) of linear springs, as described in Eq. (37), with a refer                                                               In the last investigation carried out on the mapped rectangular plate,
ence stiffness equal to 1⋅1021 N/m3 . For this case, the governing pa                                                                the structure is clamped (C) in two opposite straight edges only at half of
                                                                                                                                      the side. All the remaining parts are left free (F). Such a general
rameters have been set as ξ̄ m = 0.5, ̃ξ m = 0.01 and p = 40. It has been
                                                                                                                                      boundary condition has been set employing a super elliptic distribution
shown that, even though classical approaches reported in Eqs. (92)–(93)
                                                                                                                                      of linear springs. Also in this case, 3D FEM-based outcomes are well
are modelled together with the zigzag thickness function of Eq. (11),
                                                                                                                                      predicted only for higher order assumptions of the displacement field,
lower modes are not suitably described, since the refined three-
                                                                                                                                      according to Eq. (10), due to the complexity of the deformation for each
                                                                                                                                 24
F. Tornabene et al.                                                                                                                                                                                                        Composite Structures 309 (2023) 116542
Table 5
Modal analysis of a mapped elliptic cone [22] with variable thickness laminated with anisotropic, orthotropic and isotropic layers enforced with boundary and surface
elastic supports, employing different higher order theories. Boundary linear springs have been set so that the structure is constrained in the four outer corners. Gaussian
distribution has been employed for the description of the concentrated springs in the three principal directions located on the bottom surface. The influence of the
Winkler elastic surface support on the modal response of the shell has been shown.
  Mapped Elliptic Cone (a = 0.1 m, b = 0.15 m, α = 30 deg)
  Mode          FSDTκ=1.2
                    RS
                                          FSDTZ RS           TSDT RS           TSDTZ RS               EDZ1 RS           ED2κ=1.2              EDZ2                  ED3                  EDZ3                    ED4                EDZ4
  f [Hz]
  DOFs     5394        8091           10788          13485         8091                                                 8091                  10788                 10788                13485                   13485              16182
                                              (             )
  Winkler Foundation and Boundary    Springs BKSSS FBKSSS F
  1        1550.048    1545.468       1541.775       1534.978      1549.612                                             1543.381              1534.738              1546.916             1537.616                1545.171           1533.852
  2        1627.817    1615.560       1610.704       1608.351      1632.221                                             1613.510              1603.508              1613.805             1610.424                1611.858           1607.209
  3        1715.475    1699.338       1693.167       1690.052      1732.036                                             1699.060              1686.131              1700.041             1695.187                1697.357           1691.272
  4        1846.731    1829.989       1822.411       1819.417      1847.999                                             1830.364              1818.318              1834.821             1830.254                1832.642           1827.439
  5        1865.565    1856.857       1852.206       1848.151      1872.092                                             1847.142              1833.096              1849.995             1844.443                1846.460           1840.175
  6        2080.372    2070.928       2066.398       2063.146      2099.703                                             2082.738              2069.970              2082.549             2077.185                2080.411           2073.921
  7        2530.336    2510.172       2503.223       2497.348      2560.202                                             2508.428              2489.979              2508.469             2501.640                2504.137           2495.010
  8        2690.442    2678.320       2673.642       2666.058      2693.575                                             2685.293              2671.239              2685.107             2674.717                2683.452           2670.426
  9        3212.145    3189.888       3181.010       3175.909      3221.689                                             3183.811              3162.042              3188.761             3180.752                3184.696           3175.511
  10       3536.852    3515.432       3506.816       3499.085      3544.882                                             3520.161              3500.669              3521.786             3512.885                3519.033           3506.613
  Boundary Springs: Double - Weibull distribution (ξ̄m = ̃ξm = 0.025, p = 20),
     (k)ξ01        (k)ξ11        (k)ξ02         (k)ξ12                          (k)ξ01       (k)ξ11        (k)ξ02       (k)ξ12                                  (k)ξ01      (k)ξ11         (k)ξ02       (k)ξ12
    k 1f      = k 1f        = k 1f        = k 1f         = 1⋅1021 N/m3 , k 2f            = k 2f        = k 2f       = k 2f       = 1⋅1021 N/m3 , k 3f                    = k 3f        = k 3f       = k 3f       = 1⋅1021 N/m3
                                                                                                                                                                                  21            3
    Winkler Foundation: Gaussian distribution (ξ1n                         = 0.5, ξ2m = 0.7, Δ 1 = Δ 2 = 0.01),                       k 1f         k 2f         k 3f     = 1⋅10         N/m
                                                                                                                                       (− )         (− )         (− )
                                                                                                                                               =            =
                   (             )
  Boundary Springs BKSSS FBKSSS F
  1        636.107      633.198       630.955                                  629.962                641.705           636.200               632.398               635.855              634.772                 635.178            633.553
  2        886.733      881.083       879.205                                  878.037                901.820           883.957               876.968               883.677              881.911                 882.051            879.824
  3        1553.910     1549.092      1545.770                                 1538.800               1552.113          1547.142              1538.698              1550.429             1540.490                1548.832           1537.011
  4        1607.850     1596.657      1591.971                                 1589.430               1611.463          1596.161              1586.441              1596.092             1592.252                1594.452           1589.054
  5        1759.077     1748.778      1744.921                                 1741.473               1787.387          1753.237              1739.955              1753.599             1748.714                1750.350           1744.474
  6        1981.477     1973.769      1969.986                                 1967.395               1992.751          1985.077              1974.443              1984.744             1980.897                1983.305           1977.997
  7        2524.178     2504.285      2497.607                                 2491.814               2553.681          2502.108              2483.982              2502.538             2495.800                2498.228           2489.272
  8        2653.713     2643.590      2640.444                                 2631.672               2661.960          2654.944              2642.003              2654.120             2641.280                2652.496           2638.032
  9        2686.892     2672.832      2664.732                                 2660.252               2709.553          2676.715              2659.103              2677.485             2671.440                2674.364           2665.847
  10       3071.885     3060.884      3057.445                                 3054.049               3097.019          3084.963              3069.800              3085.868             3079.313                3083.479           3076.355
                                                                                                                        (k)ξ01       (k)ξ11        (k)ξ02         (k)ξ12                                (k)ξ01         (k)ξ11       (k)ξ02       (k)ξ12
  Boundary Springs: Double - Weibull distribution (ξ̄m = ̃ξm = 0.025, p = 20), k 1f                                              = k 1f       = k 1f        = k 1f         = 1⋅1021 N/m3 , k 2f                   = k 2f        = k 2f       = k 2f       = 1⋅1021 N/m3 ,
       (k)ξ0        (k)ξ1         (k)ξ0          (k)ξ1          21         3
      k 3f 1   =   k 3f 1    =   k 3f 2     =   k 3f 2    = 1⋅10     N/m
  Lamination Scheme: (30/iso/70), h01 = h03 = 0.003 m, h02 = 0.005 m
    Laminae Material Sequence: 1st trigonal material, 2nd steel, 3rd glass–epoxy
  Numerical Issues: CGL computational grid with IN = 33 and IM = 31 discrete points
                                          αp24 , p 4 = 20, δ 4 = 0.7
  Thickness Variation: ϕ 4 (α 1 , α 2 ) = ̃
    Isogeometric Mapping: Domain-02
mode, as shown in Fig. 9. Since a general anisotropic behavior of the                                                                         distribution of springs oriented along each principal direction. The
homogenized model has been adopted, lower modes are not symmetric.                                                                            stiffness of the translational springs has been varied. For both CCCC and
In addition, in higher modes, very complex deformations can be seen in                                                                        CFCF configurations, a reference solution has been provided by means of
the region involved by a varying structural stiffness.                                                                                        a refined 3D FEM model: a perfect alignment with the results provided
    In Table 2, we summarize the results from a parametric analysis with                                                                      by the ESL formulation is seen. In the latter, external boundary condi
respect to a three-layer mapped rectangular plate. Both West (W) and                                                                          tions have been assigned by following the kinematic approach provided
East (E) edges of the structure have been clamped with a Double –                                                                             in Eq. (69). Free (F) edges have been modelled starting from Eq. (70). It
Weibull distribution of linear springs oriented along the α 1 , α 2 and ζ                                                                     has been shown that the decrease of the natural frequencies does not
principal directions, characterized by a reference stiffness value equal to                                                                   follow the decay of the springs stiffness, since the influence of such
1⋅1021 N/m3 , according to Eq. (36). More specifically, the influence of                                                                      design parameters is evident only between 1⋅1012 N/m3 and
the constraint edge parameter ̃ ξ m has been checked on the first twenty                                                                      1⋅105 N/m3 .
mode frequencies of the plate. A reference CFCF solution has been
provided by a refined 3D FEM simulation, and it has been compared to
                                                                                                                                              5.2. Free vibration analysis of doubly-curved shells
predictions from the ESL theory, setting ̃
                                         ξ m = 1.0, with an excellent
accuracy. After this preliminary validation, different values of ̃ξ m have                                                                        We now present a series of modal analyses carried out on various
been considered, defining a reduction in the clamped area of West (W)                                                                         doubly-curved panels of different curvatures and mapped geometries.
and East (E) shell edges. As expected, a reduction of all mode frequencies                                                                    Different thickness variations along the parametric lines have been
is pointed out. Each analysis of the considered parametric investigation                                                                      selected, together with general boundary conditions. Moreover, the
has required only 16182 DOFs, embedding the EDZ4 displacement field                                                                           adopted lamination schemes are always unsymmetric and they are ob
assumption according to Eq. (10) and a CGL computational grid char                                                                           tained by layers with different material symmetries.
acterized by IN × IM = 31 × 33 discrete points. On the other hand, the                                                                            The first set of simulations have been led on a catenoid, whose
reference 3D FEM model required a total number of DOFs equal to                                                                               reference surface r(α 1 , α 2 ), according to Eq. (1), has been described in
2131302.                                                                                                                                      Fig. 3. The Domain-01 has been selected for the isogeometric mapping,
    In Table 3, we report the first ten mode frequencies for the same                                                                         along with a two-dimensional sinusoidal thickness variation character
mapped plate, constrained in West (W) and East (E) edges by a constant                                                                        ized by a single wave. Two external layers of trigonal materials are
                                                                                                                                    25
F. Tornabene et al.                                                                                                                                                                                                   Composite Structures 309 (2023) 116542
Table 6
Mode frequencies of a mapped revolution parabolic shell [22] characterized by variable thickness and subjected to general boundary linear elastic constraints. The
results have been provided by the unified ESL approach employing different higher order theories. Elastic translational springs have been distributed along the
structure in all the principal directions in order to obtain different kinds of elastic supports. A sinusoidal and power thickness variation has been set along the physical
domain.
  Mapped Revolution Parabolic Shell (R b = 0.01 m)
  Mode          FSDTκ=1.2
                    RS
                                        FSDTZ RS        TSDT RS             TSDTZ RS            EDZ1 RS              ED2κ=1.2             EDZ2             ED3                    EDZ3                  ED4               EDZ4
  f [Hz]
  DOFs         5046        7569                         10092               12615               7569                 7569                 10092            10092                  12615                 12615             15138
                           (         )
  Clamped      Portion Edge FFBKDDD F
  1            254.228     253.442                      251.400             250.166             255.194              253.219              251.153          254.078                252.498               253.465           253.126
  2            580.178     578.198                      572.224             569.314             581.174              572.583              570.042          574.145                571.356               573.640           569.342
  3            847.640     848.718                      836.266             834.021             854.515              833.135              830.294          838.488                838.343               840.777           839.220
  4            1083.495    1079.466                     1075.448            1070.430            1082.068             1074.685             1067.297         1081.045               1074.336              1082.645          1076.299
  5            1452.672    1455.553                     1441.937            1438.497            1458.070             1436.217             1429.822         1446.389               1442.937              1446.453          1442.648
  6            1725.737    1721.237                     1708.694            1706.115            1727.932             1705.015             1697.224         1712.165               1706.621              1709.146          1702.255
  7            2067.198    2057.864                     2041.066            2033.351            2061.869             2041.630             2029.336         2053.000               2040.636              2045.520          2035.037
  8            2368.604    2368.044                     2351.474            2346.518            2370.475             2349.703             2335.060         2358.254               2348.471              2354.231          2343.512
  9            2540.819    2537.859                     2522.309            2517.584            2571.080             2529.919             2517.674         2535.891               2528.701              2531.729          2527.158
  10           2764.731    2757.516                     2738.926            2730.741            2768.193             2738.190             2722.105         2752.016               2739.674              2744.294          2733.419
  Boundary Springs: Double - Weibull distribution (ξ̄m = 0.05, ̃ξm = 0.05, p = 100),
     (k)ξ01       (k)ξ11       (k)ξ0         (k)ξ1                            (k)ξ01       (k)ξ11       (k)ξ02        (k)ξ12                             (k)ξ01        (k)ξ11        (k)ξ02        (k)ξ12
    k 1f        = k 1f 2 = k 1f 2 = 1⋅1021 N/m3 ,
              = k 1f                                                         k 2f      = k 2f       = k 2f        = k 2f       = 1⋅1021 N/m3 , k 3f               = k 3f        = k 3f        = k 3f        = 1⋅1021 N/m3
                  (                       )
  Clamped Corners BKDDD BKDDD BKDDD BKDDD
  1       608.367       609.302        602.091                   612.768    598.207
                                                                              601.065                                                     599.220          604.716                604.023               607.042           597.405
  2       929.892       927.987        923.403                   933.087    918.143
                                                                              925.507                                                     919.196          930.803                925.265               932.603           930.958
  3       1300.096      1301.242       1287.168                  1300.464   1283.497
                                                                              1285.496                                                    1279.117         1293.289               1290.338              1290.816          1275.145
  4       1470.400      1472.420       1458.586                  1476.755   1455.017
                                                                              1457.640                                                    1452.390         1467.566               1468.042              1470.046          1469.916
  5       1567.612      1563.456       1556.351                  1572.613   1551.423
                                                                              1568.118                                                    1555.110         1573.407               1562.748              1574.618          1564.626
  6       1875.472      1871.351       1852.889                  1888.042   1847.929
                                                                              1858.623                                                    1849.050         1864.120               1856.628              1857.742          1827.832
  7       2077.297      2071.921       2061.927                  2086.257   2049.940
                                                                              2063.500                                                    2049.097         2077.754               2066.284              2081.600          2059.658
  8       2268.594      2264.757       2244.110                  2276.776   2238.710
                                                                              2252.121                                                    2239.775         2258.516               2252.769              2252.071          2242.865
  9       2446.282      2451.378       2427.879                  2463.422   2423.470
                                                                              2424.020                                                    2415.430         2440.258               2437.783              2435.667          2427.462
  10      2753.342      2748.848       2722.923                  2768.261   2717.979
                                                                              2731.747                                                    2717.494         2744.200               2732.683              2743.236          2735.020
  Boundary Springs: Double - Weibull distribution (ξ̄m = 0.0025, ̃ξm = 0.0025, p = 20),
     (k)ξ01       (k)ξ11       (k)ξ02        (k)ξ12                           (k)ξ0        (k)ξ11       (k)ξ02        (k)ξ12                             (k)ξ01        (k)ξ11        (k)ξ02        (k)ξ12
    k 1f      = k 1f            = 1⋅1021 N/m3 , k 2f 1 = k 2f
                           = k 1f        = k 1f                                                     = k 2f        = k 2f       = 1⋅1021 N/m3 , k 3f               = k 3f        = k 3f        = k 3f        = 1⋅1021 N/m3
                                           (         )
  Winkler Foundation and Boundary Springs BKDDD FFF
  1        178.382     177.464     183.731     188.083                                          171.496              174.672              177.501          182.866                188.902               175.605           156.040
  2        473.301     472.646     470.966     472.156                                          475.700              467.528              466.652          474.333                476.959               464.768           465.420
  3        888.584     887.151     859.486     852.495                                          893.447              853.970              837.898          873.406                866.914               844.397           834.123
  4        944.212     941.109     930.064     927.308                                          939.047              923.002              918.436          931.090                929.096               919.872           918.572
  5        1163.179    1161.490    1141.928    1136.981                                         1176.680             1147.257             1130.792         1145.224               1134.002              1126.886          1125.027
  6        1287.009    1283.271    1259.685    1259.160                                         1287.335             1240.897             1201.309         1234.109               1208.801              1209.551          1214.152
  7        1569.070    1571.806    1516.877    1505.195                                         1570.060             1547.380             1525.760         1539.533               1536.325              1460.209          1415.666
  8        1772.136    1770.177    1669.708    1651.876                                         1777.856             1663.028             1603.487         1684.156               1665.109              1633.841          1610.008
  9        1981.432    1979.153    1789.322    1786.559                                         1997.844             1812.475             1778.347         1796.458               1798.085              1783.534          1749.421
  10       2274.132    2271.357    2044.332    2037.543                                         2279.667             2108.517             2020.118         2045.607               2032.883              1949.098          1918.887
  Winkler Foundation: Gaussian distribution (ξ1n = 0.5, ξ2m = 0.75, Δ 1 = Δ 2 = 0.01), k 1f = k 2f = k 3f = 1⋅1021 N/m3
                                                                                                                                   (− )      (− )        (− )
considered in the central core of isotropic zirconia. Three different                                                                     the analysis, since the higher order expansion of the out-of-plane field
boundary conditions have been considered, namely clamped opposite                                                                         variable influences the final outcomes of the simulation. In Figs. 10-11,
      (            )                          (                   )
edges FBKSSS FBKSSS , clamped half corners BKSSS BKSSS BKSSS BKSSS and clam                                                              we report the mode shapes for the catenoidal panel under consideration.
                ( K       K
                              )                                                                                                           More in detail, in Fig. 10 we represent the eigenvectors of the catenoidal
ped edge areas FBDDD FBDDD . The first ten mode frequencies have been
reported in Table 4. A convergence of the natural frequencies can be                                                                      panel externally bounded only at a portion of the mapped edges: in this
observed for an increasing order of the kinematic expansion. Moreover,                                                                    case, the asymmetric lamination generalized stiffness matrix for the
the adoption of the zigzag hypothesis according to Eq. (11) introduces                                                                    entire lamination scheme induces a series of coupled bending and
some coupling issues between laminae, that yield different values of                                                                      warping deformations for all mode shapes. Very complicated de
mode frequencies. The influence of anisotropic materials is evident in                                                                    formations can be seen in the case of localized pointed clamping, as
the case of clamped half corners. As it can be seen from the mode fre                                                                    occurs in Fig. 11, in line with predictions from refined 3D FEM models.
quencies, there are no coupled modes despite the symmetric geometry of                                                                        In Table 5, the proposed model has been applied to a laminated
the panel in terms of physical domain, thickness distribution, and                                                                        elliptic cone, whose geometry has been mapped by means of the so-
external constraints. In the last case, the stretching effect is involved in                                                              called Domain-02. All the useful information regarding the assumed
                                                                                                                                  26
F. Tornabene et al.                                                                                                                                                                                         Composite Structures 309 (2023) 116542
Table 7
Mode frequencies of a mapped hyperbolic paraboloid [22] with variable thickness subjected to general surface and boundary linear elastic constraints. The results have
been provided by the unified ESL approach employing different higher order theories. Elastic translational springs have been distributed along the structure in all the
principal directions in order to obtain different kinds of elastic supports. A sinusoidal thickness variation has been set along the physical domain.
  Mapped Hyperbolic Paraboloid (k 1 = 4, k 2 = 4.5)
  Mode          FSDTκ=1.2
                    RS
                                        FSDTZ RS        TSDT RS     TSDTZ RS             ED1 RS           EDZ1 RS          ED2κ=1.2           EDZ2               ED3                  EDZ3               ED4        EDZ4
  f [Hz]
  DOFs     5046        7569                10092     12615       5046                                     7569             7569               10092              10092                12615              12615      15138
                                                   (              )
  Winkler Foundation and Boundary Springs BKSSS BKSSS BKSSS BKSSS
  1        98.424      98.084              95.460    75.797      95.304                                   98.040           96.706             96.624             97.672               96.737             88.879     94.774
  2        99.190      98.810              99.328    87.872      95.944                                   98.769           97.519             97.377             123.343              99.145             103.919    96.254
  3        163.913     163.601             147.329   97.896      162.390                                  162.343          161.251            160.847            123.343              132.487            109.651    159.424
  4        169.209     169.532             166.041   150.970     167.312                                  168.197          167.079            166.684            155.630              149.178            109.651    169.169
  5        171.001     170.751             179.514   165.638     169.436                                  173.057          170.376            169.958            169.720              166.962            161.667    169.169
  6        191.128     191.268             185.523   180.985     182.402                                  184.954          183.318            182.640            180.203              174.576            170.077    183.868
  7        194.041     191.946             185.523   185.205     191.252                                  191.650          188.887            188.389            186.654              185.512            180.958    183.868
  8        205.932     206.703             209.112   186.772     204.266                                  204.689          202.916            202.499            199.172              190.758            187.373    212.019
  9        215.314     215.157             209.112   209.919     210.140                                  214.636          212.645            212.142            210.623              209.243            201.355    212.019
  10       228.335     226.797             259.658   209.919     225.353                                  225.610          222.570            222.171            254.525              211.336            215.253    230.849
                   ( K K                   )
  Boundary Springs B SSS B SSS BKSSS BKSSS
  1        76.576      76.442              74.662    44.699      74.505                                   77.776           76.359             76.181             87.275               80.469             87.028     71.960
  2        95.807      95.486              94.104    57.817      92.838                                   96.008           95.000             94.865             112.039              98.463             87.137     91.285
  3        128.009     127.616             126.553   89.712      123.600                                  128.418          127.053            126.978            112.039              121.625            106.558    126.054
  4        150.066     149.778             140.519   126.934     146.704                                  150.023          149.436            149.068            129.765              132.788            106.558    149.761
  5        165.417     165.672             156.843   143.635     163.470                                  165.899          164.337            163.985            155.831              144.151            121.517    165.858
  6        169.929     169.971             169.946   157.812     168.807                                  168.477          167.372            166.964            167.808              157.245            156.159    171.702
  7        175.655     175.740             180.235   169.448     174.360                                  176.331          174.848            174.497            170.861              170.657            174.317    176.504
  8        204.377     204.972             199.389   184.342     201.178                                  203.860          202.177            201.819            194.815              189.627            174.317    205.420
  9        208.168     207.637             208.398   201.188     203.513                                  207.790          206.203            205.790            207.867              197.333            194.845    205.420
  10       211.143     210.891             229.837   205.917     207.680                                  211.115          209.285            208.770            231.098              208.722            215.060    225.285
  Winkler Foundation: Gaussian distribution (ξ1n = ξ2m = 0.5, Δ 1 = Δ 2 = 0.01), k 1f = k 2f = k 3f = 1⋅1021 N/m3
                                                                                                                   (− )     (− )       (− )
                                                                                                                       0           1          0             1                                      0           1       0      1
    Boundary Springs: super elliptic distribution (ξ̄m = 0.5, ̃ξm = 0.01, p = 40), k 1f 1 = k 1f 1 = k 1f 2 = k 1f 2 = 1⋅1021 N/m3 , k 2f 1 = k 2f 1 = k 2f 2 = k 2f 2 = 1⋅1021 N/m3 ,
                                                                                    (k)ξ     (k)ξ     (k)ξ     (k)ξ                   (k)ξ     (k)ξ     (k)ξ     (k)ξ
geometry have been reported in Fig. 4. In this case, a power thickness                                                          not included in the model, the fourth order of the kinematic expansion
variation has been selected along the meridian direction of the shell, so                                                       allows for the convergence of results. Some mode shapes of the present
that a mass concentration is held near the West (W) side of the structure.                                                      investigations are collected in Fig. 12. Namely, the influence of the
The laminate consists of three layers of trigonal material (89), isotropic                                                      surface support is evident from the bending component of the structural
steel and orthotropic glass–epoxy (90). Two different cases for boundary                                                        deformation in all the reported modes.
conditions have been implemented. In the first one, identified by                                                                   Moreover, the isogeometric mapping by means of the Domain-02 has
( K            )
 B SSS FBKSSS F , the structure has been fixed in the four corners by setting                                                   been applied to a revolution parabolic shell. In this case, the shell
a super elliptic linear spring distribution along West (W) and East (E)                                                         thickness has been varied following a combination of power and one-
edges. Moreover, a Gaussian distribution of Winkler elastic support,                                                            wave sinusoidal analytical expressions. The lamination scheme em
defined according to Eq. (30), has been set, accounting for a point                                                             beds an orthotropic graphite-epoxy (91), a triclinic material (88) and a
support on the bottom surface. The second configuration of boundary                                                             trigonal continuum (89). Three different external constraints have been
conditions is exactly the same as the previous one, but does not involve                                                        considered, and the first ten mode frequencies have been computed by
any surface support. With particular reference to the latter constraint                                                         means of different displacement field assumptions in Table 6. In the first
                                                                                                                                          (         )
settings, lower mode frequencies are slightly influenced by the                                                                 analysis, FFBKDDD F boundary conditions are considered. The structure
assumption of the zigzag thickness function defined in Eq. (11), even                                                           has been fixed on two portions of the East (E) edge, thanks to the use of
though a higher order power expansion is held within the field variable                                                         the Double – Weibull distribution of linear springs, introduced in Eq.
features of Eq. (10). Moreover, when the Winkler elastic foundation is                                                          (36). The corresponding mode shapes are reported in Fig. 13. The out-of-
                                                                                                                           27
F. Tornabene et al.                                                                                                                                         Composite Structures 309 (2023) 116542
Table 8
Mode frequencies of a mapped helicoid [22] characterized by variable thickness and subjected to general boundary linear elastic constraints. The results have been
provided by the unified ESL approach employing different higher order theories. Elastic translational springs have been distributed along the structure in all the
principal directions in order to obtain different kinds of elastic supports. A sinusoidal thickness variation has been set along the physical domain.
  Mapped Helicoid (a = 1 m)
  Mode         FSDTκ=1.2
                   RS
                                FSDTZ RS         TSDT RS        TSDTZ RS         ED1 RS         EDZ1 RS         ED2κ=1.2       EDZ2             ED3          EDZ3         ED4          EDZ4
  f [Hz]
  DOFs       5046          7569           10092        12615          5046          7569                        7569           10092            10092        12615        12615        15138
                                                                        (            )
  In-Plane General Constraint – Spring Parabolic Thickness distribution BKDD FBKDD C
  1          122.010       122.281        122.064      122.056        113.414       125.041                     122.610        122.444          122.891      122.800      123.025      123.235
  2          144.034       144.421        144.191      144.061        133.676       150.790                     144.538        144.337          144.774      144.606      144.779      144.625
  3          155.531       155.870        155.539      155.359        145.639       163.489                     155.962        155.753          155.971      155.660      156.621      156.525
  4          182.544       183.399        182.722      181.901        175.474       195.157                     183.004        182.792          183.410      182.664      182.360      181.226
  5          223.767       224.527        223.942      223.649        209.370       232.842                     224.623        224.301          225.060      224.689      224.854      224.520
  6          246.950       247.376        246.869      246.533        225.446       256.115                     248.081        247.650          248.219      247.855      247.297      246.786
  7          258.643       259.449        258.932      258.834        244.601       273.562                     259.795        259.400          260.326      260.174      259.853      259.867
  8          309.549       310.914        309.924      309.426        298.628       322.230                     310.420        310.013          311.069      310.487      310.663      310.151
  9          344.798       346.161        345.140      344.588        327.748       365.952                     345.945        345.478          346.793      346.243      346.566      346.481
  10         360.214       361.418        360.662      360.653        343.471       383.167                     361.497        361.024          362.215      362.008      362.538      362.481
                                                                                                  0        1                           0         1
  Boundary Springs: Double - Weibull distribution (ξ̄m = 1, ̃ξm = 0.27, p = 1000), k 1f 1 = k 1f 1 = 1⋅1021 N/m3 , k 2f 1 = k 2f 1 = 1⋅1021                  N/m3
                                                                                    (k)ξ     (k)ξ                   (k)ξ     (k)ξ
                                                                  (           )
  General Constraint – Spring Parabolic Thickness distribution BKDDD FBKDDD C
  1         124.367         124.646        124.501        124.775       116.903       128.263       125.355     125.176        125.691                       125.841      125.867      126.450
  2         150.562         150.935        150.719        150.565       140.936       157.216       151.306     151.090        151.540                       151.358      151.507      151.262
  3         176.871         177.293        176.881        176.278       166.492       185.302       177.391     177.157        177.380                       176.740      177.675      176.750
  4         207.743         208.985        208.278        208.347       203.891       227.361       208.635     208.402        209.636                       209.687      208.593      208.682
  5         256.006         256.771        256.310        256.547       237.156       270.764       257.323     256.908        257.853                       257.924      257.196      257.212
  6         287.992         289.339        288.545        288.202       272.950       300.008       289.288     288.803        290.096                       289.691      290.024      289.494
  7         313.138         314.418        313.618        313.348       298.497       324.110       314.250     313.833        315.006                       314.707      314.673      314.346
  8         357.031         358.893        357.579        356.948       338.933       382.403       359.092     358.556        360.238                       359.471      360.031      359.168
  9         362.698         363.760        362.914        362.576       342.997       392.605       364.206     363.716        364.783                       364.316      364.836      364.302
  10        410.670         412.118        411.230        411.021       386.741       430.702       412.125     411.519        413.013                       412.636      412.756      412.164
                                                                                                  0        1                           0         1
  Boundary Springs: Double - Weibull distribution (ξ̄m = 1, ̃ξm = 0.27, p = 1000), k 1f 1 = k 1f 1 = 1⋅1021 N/m3 , k 2f 1 = k 2f 1 = 1⋅1021 N/m3
                                                                                    (k)ξ     (k)ξ                   (k)ξ     (k)ξ
plane power ESL formulation provides a rapid convergence of both low                                   found in Table 8. Referring to the former, the super elliptic distribution
and high natural frequencies. In the second case study, the four corners                               of linear springs have been applied to the structure, together with a
have been fixed employing the above-mentioned Double – Weibull                                         concentrated Winkler support in the middle of the structure. The
distribution. Also in this case, a rapid convergence of results is noticed as                          importance of such elastic foundation on lower modes can be seen from
the kinematic order increases. The influence of the zigzag function                                    the natural frequencies values reported in Table 7. The first nine ei
introduced in Eq. (11), together with a higher order displacement field                                genvectors of this last example at issue, calculated by means of the EDZ4
assumption, is more evident in the last configuration identified as                                    theory, are reported in Figs. 15-16. In the case of clamped half edges
( K        )                                                                                           ( K           )
 B DDD FFF , where the shell is fixed at two extreme points of the West                                 B SSS FFBKSSS , the zigzag thickness function according to Eq. (11) and
(W) edge, together with a point constraint acting at the bottom surface,                               the use of the higher order theory is crucial, since it allows one to
obtained from a Gaussian distribution of Winkler springs in each prin                                 consider an additional vibration mode of the structure that cannot be
cipal direction. Actually, lower frequencies decrease as the displacement                              otherwise detected by lower order theories, as well as the coupling of
field unified formulation becomes more complicated. For this last case,                                some mode shapes. In Table 8, a free vibration analysis has been per
the first nine mode shapes have been reported in Fig. 14.                                              formed for the previously described helicoidal panel. Apart from the
    A mapping by means of the Domain-03 has been introduced in Fig. 5,                                 cantilevered boundary conditions configuration (FFFC), two arrange
characterized by four curved edges. It has been applied to a hyperbolic                                ments of general external constraints have been considered, taking into
paraboloid and helicoidal panel. These two structures contain a one-                                   account a parabolic-through-the-thickness distribution of linear elastic
                                                                                                                                               (           )
wave sinusoidal thickness variation, according to Eqs. (4), (6). The                                   springs. In the first case, denoted by BKDD FBKDD C , two portions of the
first ten mode frequencies of the above panels have been computed for                                  West (W) and East (E) side of the shell have been clamped by a proper
various external restraints taking into account different higher order                                 setup of a Double – Weibull distribution of linear springs, according to
theories. The results of the analysis held on the hyperbolic paraboloid                                Eq. (36). As it can be seen, general boundary conditions have been set
are arranged in Table 7, whereas those referred to the helicoid can be
                                                                                                  28
F. Tornabene et al.                                                                                                                                          Composite Structures 309 (2023) 116542
Table 9
Mode frequencies of a mapped pseudosphere [22] of revolution with a tractrix meridian characterized by variable thickness and subjected to general surface and
boundary linear elastic constraints. The results have been provided by the unified ESL approach employing different higher order theories. Elastic translational springs
have been distributed along the structure in all the principal directions in order to obtain different kinds of elastic supports. A power thickness variation has been set
along the physical domain.
  Mapped Pseudosphere (a = 1.5 m)
  Mode            FSDTκ=1.2
                      RS
                                     FSDTZ RS            TSDT RS           TSDTZ RS             EDZ1 RS             ED2κ=1.2     EDZ2         ED3          EDZ3         ED4            EDZ4
  f [Hz]
  DOFs       5046          7569                          10092             12615                7569                7569         10092        10092        12615        12615          15138
                         (              )
  Three-Points Constraint BKDDD FFBKDDD
  1          132.566       132.770                       132.200           131.874              132.902             131.802      130.959      132.356      131.615      133.053        132.249
  2          238.474       236.532                       234.838           233.291              236.528             233.619      231.397      235.150      232.996      234.770        233.047
  3          332.070       332.512                       330.315           329.418              335.442             332.325      330.500      333.388      331.697      332.671        331.470
  4          399.071       402.328                       399.491           398.348              405.907             399.097      396.198      406.167      403.191      404.538        402.620
  5          442.837       442.728                       438.437           437.777              447.405             442.220      440.167      442.492      440.587      441.789        440.638
  6          570.894       578.494                       572.572           569.979              582.100             573.912      570.202      581.800      577.484      580.863        578.567
  7          601.917       601.763                       595.587           592.836              606.063             597.130      592.946      608.541      604.379      606.347        604.568
  8          712.135       712.116                       706.418           705.596              719.744             714.024      711.528      713.190      711.618      712.641        711.929
  9          739.361       745.584                       739.310           737.535              750.355             741.149      737.377      748.287      744.376      748.139        746.395
  10         780.042       789.332                       781.887           781.829              800.816             790.836      788.772      800.102      799.030      795.918        795.055
  Boundary Springs: North (N) and West (W) edges: Double - Weibull distribution (ξ̄m = ̃
                                                                                       ξm = 0.0025, p = 20),
     (k)ξ01       (k)ξ02                     (k)ξ0       (k)ξ0                        (k)ξ0        (k)ξ0
    k 1f        = 1⋅1021 N/m3 , k 2f 1 = k 2f 2 = 1⋅1021 N/m3 , k 3f 1 = k 3f 2 = 1⋅1021 N/m3
              = k 1f
                         (             )
  Clamped Quarter Edges FBKSSS BKSSS F
  1         272.570         278.327           272.816         273.021       259.609    256.682                                   255.536      264.307      263.683      254.173        252.797
  2         309.561         311.467           308.358         307.537       313.686    308.425                                   306.092      312.665      310.500      311.611        310.251
  3         384.621         388.170           385.454         386.775       398.975    399.614                                   399.870      398.289      399.124      395.606        396.102
  4         466.920         480.373           477.290         476.160       488.540    479.243                                   478.067      494.980      492.393      490.706        491.321
  5         492.680         493.478           493.344         491.729       490.419    485.016                                   480.661      495.209      492.393      494.237        491.759
  6         561.090         563.696           559.288         557.991       563.670    551.607                                   547.505      562.646      558.771      562.490        559.953
  7         633.202         636.755           627.435         624.265       643.197    630.992                                   626.341      640.611      635.137      639.886        636.667
  8         780.971         787.750           780.682         778.090       798.389    786.858                                   781.564      795.631      789.857      791.929        788.251
  9         865.565         879.567           870.640         868.140       887.157    870.595                                   866.051      892.395      887.344      890.310        888.054
  10        898.530         908.870           903.850         908.581       924.708    918.438                                   918.602      928.270      930.543      923.326        922.774
  Boundary Springs: Super elliptic distribution (ξ̄m = 1.0, ̃ξm = 0.25, p = 20),
     (k)ξ11       (k)ξ12                     (k)ξ11      (k)ξ12                       (k)ξ11       (k)ξ12
    k 1f      = k 1f       = 1⋅1021 N/m3 , k 2f       = k 2f      = 1⋅1021 N/m3 , k 3f         = k 3f       = 1⋅1021 N/m3
  Two Clamped Edges (FCFF)
  1          312.388            318.739               317.535        317.052        298.024                         297.305      295.887      304.864      303.565      302.887        301.493
  2          501.598            516.300               512.712        512.098        526.253                         518.415      516.766      533.852      532.146      531.888        531.287
  3          560.495            560.956               558.135        558.001        567.142                         560.430      557.773      562.957      561.101      561.930        561.106
  4          657.637            657.953               650.755        649.248        666.575                         654.185      649.868      661.784      657.624      659.196        657.187
  5          822.181            837.681               829.859        828.466        845.241                         831.458      828.114      845.600      841.889      843.346        841.677
  6          895.778            897.798               890.405        888.528        907.678                         894.720      888.980      903.705      898.435      899.953        897.038
  7          1019.937           1024.393              1012.275       1009.327       1034.845                        1014.433     1007.568     1030.200     1023.254     1026.050       1022.544
  8          1051.266           1053.173              1038.904       1035.615       1064.888                        1041.884     1034.070     1053.890     1045.730     1049.198       1044.863
  9          1196.403           1212.604              1204.554       1201.236       1223.266                        1205.197     1197.853     1228.676     1220.439     1223.483       1219.211
  10         1289.146           1285.160              1269.122       1265.109       1303.028                        1275.534     1265.023     1281.371     1270.670     1275.454       1269.686
  Lamination Scheme: (30/iso/iso/45), h01 = h02 = h03 = h04 = 0.007 m
    Laminae Material Sequence: 1st trigonal material, 2nd aluminium, 3rd steel, 4th trigonal                        material
  Numerical Issues: CGL computational grid with IN = IM = 31 discrete points
  Thickness Variation: ϕ i (α 1 , α 2 ) =¯αpi i , p i = p 1 = p 2 = 10, δ 1 = δ 2 = 1.2
    Isogeometric Mapping: Domain-04
only along the in-plane principal directions α 1 , α 2 . A similar configu                                          have been reported in Table 9. In particular, three-points constraints
        (             )                                                                                              ( K            )                          (             )
ration BKDDD FBKDDD C has also been considered in which linear springs                                                B DDD FFBKDDD , clamped quarter edges FBKSSS BKSSS F and (FCFF) con
are distributed along α 1 , α 2 , ζ. As expected, mode frequencies do not                                            figurations have been studied. As it can be seen, if in the first configu
differ too much between these two case studies. However, as it can be                                                ration there is a very rapid convergence of results as the higher order
                                                                                                                                                      (            )
seen from a comparison between Figs. 17-18, the corresponding mode                                                   model gets complicated, in the FBKSSS BKSSS F stable results are reached
shapes differ due to the presence of the out-of-plane deformation in the                                             for higher and lower modes when the ED4 theory is assumed, regardless
former configuration.                                                                                                of the use of the zigzag behaviour set in Eq. (11). When the structure is
    The last two examples have been obtained from the physical domain                                                clamped at two edges (FCFF), results provided by higher order theories
mapping declared in Fig. 6, characterized by two curved sides and two                                                differ from those provided by classical ESL approaches, underlying the
straight edges. In particular, a pseudosphere and a revolution hyperbolic                                            importance of the through-the-thickness stretching effect for all mode
hyperboloid have been blended with the considered mapped domain. A                                                   shapes, as it appears in Figs. 19-20. In Table 10 one can find the modal
power thickness variation has been selected along all parametric lines.                                              response of a revolution hyperbolic hyperboloid mapped with the
The mapped pseudosphere has been laminated with four generally ori                                                  domain reported in Fig. 6, tackled by means of the ESL formulation
ented layers, accounting for two outer skins of trigonal material (89) and                                           presented in this work, and various higher order theories. The lamina
a non-homogeneous central core made of aluminium and steel isotropic                                                 tion scheme consists of two external layers of glass–epoxy (90) and a
sheet, whose elastic properties have been defined in the previous sec                                               central trigonal layer (89). The first boundary conditions layup con
tions. Mode frequencies for three different sets of boundary conditions                                              siders a Double – Weibull distribution of linear elastic springs accounted
                                                                                                               29
F. Tornabene et al.                                                                                                                                                         Composite Structures 309 (2023) 116542
Table 10
Mode frequencies of a mapped revolution hyperbolic hyperboloid [22] characterized by variable thickness and subjected to general boundary linear elastic constraints.
The results have been provided by the unified ESL approach employing different higher order theories. Elastic translational springs have been distributed along the
structure in all the principal directions in order to obtain different kinds of elastic supports. A sinusoidal thickness variation has been set along the physical domain.
  Mapped Revolution Hyperbolic Hyperboloid (a = 1 m, c = 2 m)
  Mode       FSDTκ=1.2
                 RS
                                 FSDTZ RS            TSDT RS          TSDTZ RS   ED1 RS    EDZ1 RS    ED2κ=1.2               EDZ2             ED3            EDZ3        ED4             EDZ4
  f [Hz]
  DOFs      5046        7569             10092                        12615      5046      7569       7569                   10092            10092          12615       12615           15138
                    (                       )
  Clamped   Corners BKDDD BKDDD BKDDD BKDDD
  1         161.440     161.651          161.160                      160.968    160.750   163.155    161.606                160.641          161.820        161.642     161.215         161.010
  2         202.591     202.995          202.509                      202.038    201.629   203.625    202.652                201.409          203.243        202.665     202.443         201.749
  3         209.595     209.913          209.088                      208.584    208.767   211.645    209.533                208.230          209.915        209.300     209.070         208.453
  4         259.500     259.707          258.768                      257.690    258.544   262.512    259.295                257.705          259.893        258.666     258.641         257.461
  5         285.687     286.041          285.011                      284.187    284.654   287.798    285.135                283.474          285.814        285.026     284.572         283.624
  6         316.785     317.860          316.827                      315.973    315.358   319.213    316.580                314.941          318.010        317.076     316.718         315.607
  7         349.227     349.622          348.556                      347.374    347.749   350.813    348.536                346.455          349.609        348.232     348.271         346.752
  8         380.155     380.374          378.157                      376.829    379.509   384.885    378.423                376.312          379.413        378.137     377.175         375.910
  9         433.661     434.217          432.253                      430.056    432.630   439.350    432.307                429.723          433.987        431.874     431.552         429.230
  10        482.159     483.225          482.213                      480.466    479.777   487.075    482.667                479.415          484.223        482.250     482.087         480.112
                                                                                                           (k)ξ01       (k)ξ11       (k)ξ02      (k)ξ12                     (k)ξ01       (k)ξ11       (k)ξ02       (k)ξ12
  Boundary Springs: Double - Weibull distribution (ξ̄m = 0.0025, ̃ξm = 0.0025,                p = 20), k 1f         = k 1f       = k 1f       = k 1f      = 1⋅1021 N/m3 , k 2f       = k 2f       = k 2f       = k 2f       = 1⋅
                      (k)ξ01       (k)ξ11       (k)ξ02       (k)ξ12
    1021 N/m3 , k 3f           = k 3f        = 1⋅1021 N/m3
                                            = k 3f       = k 3f
                                                           (              )
  Clamped Half Edges – Spring Linear Thickness distribution BKSSS FFBKSSS
  1       50.687        51.003      50.700     50.563        50.568         51.411                    50.378                 50.138           50.762         50.576      50.408          50.451
  2       51.793        52.186      51.820     51.830        51.736         52.307                    51.528                 51.464           51.866         51.855      51.647          51.654
  3       112.887       113.254     112.749    112.383       112.571        113.970                   112.490                111.960          112.945        112.703     112.410         112.073
  4       176.204       176.822     175.631    174.943       175.858        178.193                   175.138                174.557          176.096        175.364     175.105         174.346
  5       261.445       261.885     260.111    258.888       261.053        264.247                   259.636                258.482          260.839        259.460     259.119         257.792
  6       276.455       277.448     276.563    275.915       275.668        280.739                   276.439                275.232          277.742        276.975     276.456         275.568
  7       362.348       362.668     361.881    360.852       360.115        364.492                   362.502                359.796          363.030        361.889     361.531         360.406
  8       379.344       379.263     377.793    376.865       377.848        382.324                   378.193                376.235          378.737        377.713     377.033         375.951
  9       403.351       403.671     403.450    402.517       401.821        404.905                   403.822                401.086          404.521        403.410     403.411         402.223
  10      416.406       416.147     414.658    413.510       415.780        423.136                   415.703                413.749          416.201        414.933     413.981         412.810
                                                                                                  0           0                                   0           0                           0            0
  Boundary Springs: Super elliptic distribution (ξ̄m = 0, ̃ξm = 0.47, p = 1000), k 1f 1 = k 1f 2 = 1⋅1021 N/m3 , k 2f 1 = k 2f 2 = 1⋅1021 N/m3 , k 3f 1 = k 3f 2 = 1⋅1021 N/m3
                                                                                  (k)ξ     (k)ξ                   (k)ξ     (k)ξ                   (k)ξ     (k)ξ
in Eq. (36), for all the shell edges, leading to four clamped corners. In                                    convergence, since ED1 and EDZ1 theories do not provide results com
this case, the warping phenomenon occurs in bending modes and                                                parable with other higher order theories.
therefore a linear through-the-thickness variation of the displacement
field components is inadequate for a proper derivation of modal eigen                                       6. Conclusions
values. Such warping effect is not related to interlaminar issues, as it
appears in Fig. 21, since not only the ED1, but also the EDZ1 theory                                             In the present work, an ESL formulation for the free vibration anal
provide results not comparable with other higher order theories. The                                         ysis of anisotropic doubly-curved shells with arbitrary geometry sub
second simulation on the hyperbolic hyperboloid has been performed by                                        jected to general boundary conditions has been proposed. The mapping
                                                          (            )
clamping the area near a corner of the structure, setting BKSSS FFBKSSS . In                                 of the physical domain has been obtained by means of a NURBS-based
addition, a linear through-the-thickness distribution of linear springs has                                  set of blending functions. A general thickness variation has been set
been considered. It is worth noting that the use of a completely aniso                                      along the principal parametric lines of the structure, whereas a general
tropic stacking sequence induces unsymmetric modes, as it can be seen                                        anisotropic material has been considered in each layer of the lamination
from mode frequencies values, and mode shapes in Fig. 22. The last free                                      scheme. Unlike previous studies regarding general boundary conditions
vibration analysis accounts for cantilever boundary conditions (FFFC).                                       within an ESL model, non-conventional constraints have been modelled
In this case, the structure exhibits some stretching effects in all modes,                                   along each lateral surface of the structural solid starting from an arbi
with clear discrepancies with respect to classical ESL theories. We recall                                   trary distribution of linear springs not only alongside the boundary of
that they account for constant through-the-thickness out-of-plane                                            the two-dimensional physical domain, but on the entire lateral surface of
displacement components. Moreover, a parabolic (N = 2) displacement                                          the three-dimensional shell along both in-plane and out-of-plane di
field assumption is at least required to ensure a satisfactory                                               rections. In this way, the calibration of the distribution governing
                                                                                                      30
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 10. Mode shapes of a mapped catenoid constrained in two edges with general boundary conditions. A sinusoidal thickness variation has been selected. The
external constraint has been obtained from a super elliptic distribution of translational springs along the edges of the mapped geometry by properly selecting the
governing parameters. Corresponding mode frequencies have been collected in Table 4.
                                                                                31
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 11. Mode shapes of a mapped catenoid constrained in the middle point of each edge of the structure. A sinusoidal thickness variation has been selected. The
external constraint has been obtained from a super elliptic distribution of translational springs along the edges of the mapped geometry by properly selecting the
governing parameters. Corresponding mode frequencies have been collected in Table 4.
                                                                                32
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 12. Mode shapes of a mapped elliptic cone constrained in its four outer corners. The influence of a concentrated spring applied on the bottom surface has been
pointed out. Corresponding mode frequencies have been collected in Table 5.
                                                                                33
F. Tornabene et al.                                                                                                             Composite Structures 309 (2023) 116542
Fig. 13. Mode shapes of a mapped revolution parabolic shell constrained in two portions of an edge. A sinusoidal thickness variation has been selected. The external
constraint has been obtained from a Double - Weibull distribution of translational springs along the edges of the mapped geometry by properly selecting the gov
erning parameters. Corresponding mode frequencies have been collected in Table 6.
                                                                                34
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 14. Mode shapes of a mapped revolution parabolic shell constrained in two portions of an edge and with a surface spring applied on the bottom surface. A
sinusoidal thickness variation has been selected. The external constraint has been obtained from a Double - Weibull distribution of translational springs along the
edges of the mapped geometry by properly selecting the governing parameters, as well as a Gaussian distribution for the definition of the restraint on the bottom
surface. Corresponding mode frequencies have been collected in Table 6.
                                                                                35
F. Tornabene et al.                                                                                                               Composite Structures 309 (2023) 116542
Fig. 15. Mode shapes of a mapped laminated hyperboloid paraboloid with variable thickness and general boundary conditions. A sinusoidal thickness variation has
been selected. The external constraint has been obtained from a super elliptic distribution of translational springs distributed along the three principal directions
applied on all the sides of the structure, and from a Gaussian distribution of linear springs applied on the bottom surface. Corresponding mode frequencies, as well as
other analysis information, have been collected in Table 7.
                                                                                  36
F. Tornabene et al.                                                                                                              Composite Structures 309 (2023) 116542
Fig. 16. Mode shapes of a mapped laminated hyperbolic paraboloid with variable thickness and general boundary conditions. A sinusoidal thickness variation has
been selected. The external constraint has been obtained from a super elliptic distribution of translational springs distributed along the three principal directions
applied on all the sides of the structure. Corresponding mode frequencies, as well as other analysis information, have been collected in Table 7.
parameters has led to the definition of point constraints without recur              Hamiltonian Principle directly in the strong form and they have been
ring to the classical domain decomposition methodology. Furthermore,                  numerically tackled by means of the GDQ method. In this way, the so
a Winkler elastic foundation with general distribution has been taken                 lution has been derived in an efficient way even though a reduced mass
into account on the top and the bottom surfaces of the shell. The                     storage has been required. The model has been validated with respect to
fundamental governing equations have been derived employing the                       refined 3D FEM simulations, with a perfect match between the relevant
                                                                                 37
F. Tornabene et al.                                                                                                             Composite Structures 309 (2023) 116542
Fig. 17. Mode shapes of a mapped helicoid constrained with in-plane elastic springs on two sides of the physical domain. A sinusoidal thickness variation has been
selected. The external constraint has been obtained from a Double - Weibull distribution of translational springs along the edges of the mapped geometry by properly
selecting the governing parameters. Corresponding mode frequencies have been collected in Table 8.
solutions. Different physical domain mappings have been investigated,                obtain very accurate results, despite the limited reduced number of
for structures with various curvatures, lamination schemes and                       discrete points implemented in each example.
displacement field assumptions. The influence of the selection of gov
erning parameters for the assessment of surface and lateral boundary                 Declaration of Competing Interest
conditions has been pointed out. The proposed formulation provides an
ESL formulation with 3D capability for a structural problem character                   The authors declare that they have no known competing financial
ized by very complex geometry, stacking sequence and external con                   interests or personal relationships that could have appeared to influence
straints. To sum up, the advanced numerical algorithm allows one to                  the work reported in this paper.
                                                                                38
F. Tornabene et al.                                                                                                             Composite Structures 309 (2023) 116542
Fig. 18. Mode shapes of a mapped helicoid constrained with elastic springs on two sides of the physical domain. A sinusoidal thickness variation has been selected.
The external constraint has been obtained from a Double - Weibull distribution of translational springs along the edges of the mapped geometry by properly selecting
the governing parameters. Corresponding mode frequencies have been collected in Table 8.
                                                                                39
F. Tornabene et al.                                                                                                           Composite Structures 309 (2023) 116542
Fig. 19. Mode shapes of a mapped revolution pseudosphere constrained in the neighbourhood of a shell corner. A power thickness variation has been selected. The
external constraint has been obtained from a super elliptic distribution of translational springs along South (S) and East (E) edges of the mapped geometry. Cor
responding mode frequencies have been collected in Table 9.
                                                                               40
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 20. Mode shapes of a mapped revolution pseudosphere constrained in three corners. A power thickness variation has been selected. The external constraint has
been obtained from a Double - Weibull distribution of translational springs along the edges of the mapped geometry by properly selecting the governing parameters.
Corresponding mode frequencies have been collected in Table 9.
                                                                               41
F. Tornabene et al.                                                                                                          Composite Structures 309 (2023) 116542
Fig. 21. Mode shapes of a mapped revolution hyperbolic hyperboloid constrained in its corners. A sinusoidal thickness variation has been selected. The external
constraint has been obtained from a Double - Weibull distribution of translational springs along the edges of the mapped geometry by properly selecting the gov
erning parameters. Corresponding mode frequencies have been collected in Table 10.
                                                                               42
F. Tornabene et al.                                                                                                            Composite Structures 309 (2023) 116542
Fig. 22. Mode shapes of a mapped revolution hyperbolic hyperboloid constrained in its corners. A sinusoidal thickness variation has been selected. The external
constraint has been obtained from a super elliptic distribution of translational springs along the edges of the mapped geometry by properly selecting the governing
parameters. Corresponding mode frequencies have been collected in Table 10.
                                                                                43
F. Tornabene et al.                                                                                                                                   Composite Structures 309 (2023) 116542
Data availability
Appendix I
     Second order derivatives with respect to principal coordinates α 1 , α 2 should be assessed by means of their corresponding ones with respect to
natural coordinates ξ 1 ,ξ 2 ∈ [ − 1, 1] × [ − 1, 1], employing the NURBS-based blending function presented in Eqs. (54)–(55). For the direct computation
of second-order derivatives set in Eqs. (64)–(66), the following definitions are required:
                  (               ( )2                                        )
              1    ∂α2 ∂2 α2       ∂α2 det(J)ξ1 ∂α2 ∂2 α2 ∂α2 ∂α2 det(J)ξ2
ξ1,α1 α1 =                     −                   −          +                                                                                     (A.1)
           det(J)2 ∂ξ2 ∂ξ1 ∂ξ2     ∂ξ2 det(J) ∂ξ1 ∂ξ22 ∂ξ1 ∂ξ2 det(J)
                       (                   (         )2                                                )
               1           ∂α1 ∂2 α1           ∂α1        det(J)ξ1 ∂α1 ∂2 α1 ∂α1 ∂α1 det(J)ξ2
ξ1,α2 α2 =                             −                          −         +                                                                                                        (A.2)
             det(J)2       ∂ξ2 ∂ξ1 ∂ξ2         ∂ξ2         det(J)   ∂ξ1 ∂ξ22 ∂ξ1 ∂ξ2 det(J)
                       (                                                                                     )
               1               ∂α2 ∂2 α1    ∂α2 ∂α1 det(J)ξ1 ∂α2 ∂2 α1 ∂α2 ∂α1 det(J)ξ2
ξ1,α1 α2 =         2
                           −              +                 +         −                                                                                                              (A.3)
             det(J)            ∂ξ2 ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ2 det(J) ∂ξ1 ∂ξ22 ∂ξ1 ∂ξ2 det(J)
                       (                                                          (          )2              )
               1             ∂α2 ∂2 α2 ∂α2 ∂α2 det(J)ξ1 ∂α2 ∂2 α2                      ∂α2        det(J)ξ2
ξ2,α1 α1 =                 −          +                +          −                                                                                                                  (A.4)
             det(J)2         ∂ξ2 ∂ξ21 ∂ξ2 ∂ξ1 det(J) ∂ξ1 ∂ξ1 ∂ξ2                       ∂ξ1        det(J)
                       (                                                          (          )2              )
               1               ∂α1 ∂2 α1 ∂α1 ∂α1 det(J)ξ1 ∂α1 ∂2 α1                    ∂α1        det(J)ξ2
ξ2,α2 α2 =                 −            +                +          −                                                                                                                (A.5)
             det(J)2           ∂ξ2 ∂ξ21 ∂ξ2 ∂ξ1 det(J) ∂ξ1 ∂ξ1 ∂ξ2                     ∂ξ1         det(J)
                       (                                                                                     )
               1               ∂α2 ∂2 α1    ∂α1 ∂α1 det(J)ξ1 ∂α2 ∂2 α1 ∂α2 ∂α1 det(J)ξ2
ξ2,α1 α2 =         2
                           −              +                 +         −                                                                                                              (A.6)
             det(J)            ∂ξ1 ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ1 det(J) ∂ξ2 ∂ξ21 ∂ξ1 ∂ξ1 det(J)
where det(J)ξ1 and det(J)ξ2 are the derivatives of the determinant of the Jacobian matrix J of the blending functions reported in Eqns. (54)–(55) with
respect to the natural coordinates ξ 1 and ξ 2 , respectively:
                ∂α1 ∂2 α2    ∂α2 ∂2 α1   ∂α2 ∂2 α1 ∂α1 ∂2 α2
 det(J)ξ1 =                −           +          −
                ∂ξ1 ∂ξ1 ∂ξ2 ∂ξ1 ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ21 ∂ξ2 ∂ξ21
                                                                                                                                                                                     (A.7)
             ∂α1 ∂2 α2    ∂α2 ∂2 α1   ∂α2 ∂2 α1 ∂α1 ∂2 α2
det(J)ξ2 = −            +           −          +
             ∂ξ2 ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ1 ∂ξ2 ∂ξ1 ∂ξ22 ∂ξ1 ∂ξ22
Appendix II
   An extended version of the generalized mass coefficients of the higher order mass matrix M(τη) is now provided, together with the fundamental
coefficients of L(τη) occurring in Eq. (85), expressed in their discrete form for each τ, η = 0, ..., N + 1.
                   (      )
   All the values M(τη) ij , for i, j = 1, 2, 3 of the generalized mass matrix occurring in Eq. (53) are calculated in each point of the CGL discrete grid
                                                                                                                    (τη)αi αj
already defined in Eq. (72), and they are collected in a IN IM × IN IM matrix denoted with M
                                                                                           ̃
                                                                                             ij                                 for i, j = 1, 2, 3:
             ⎧
             ⎨              0                for i ∕
                                                   =j
̃ (τη)αi αj = →0(τη)α i α j
M                                                                                                                                                                                    (A.8)
  ij
             ⎩ I            ◦ϛα 1 (0)α 2 (0)
                                             for i = j
                 ij
      →f(τη)αi α j
where L ij         is the vectorized form of a IN × IM matrix collecting the values assumed at each point of the computational grid by Eq. (28).
                                      (    )
  On the other hand, each term L(τη) ij , for i, j = 1, 2, 3, of the fundamental matrix L(τη) has been computed in each discrete point of the already
                                                                                                           ̃ (τη)αi αj :
mentioned CGL grid according to the following expression, leading to the discrete fundamental coefficients L ij
                                                                                                      44
F. Tornabene et al.                                                                                                                        Composite Structures 309 (2023) 116542
    As can be seen, the formulation adopted in the previous equation embeds the previously introduced by-column vectorization of each term Lij(k) i
                                                                                                                                                                          (τη)α αj
with k = 1, ..., 6 related to the derivatives along the principal coordinates α i , α j = α 1 , α 2 and the employment of the well-known Hadamard’s Product,
identified with the symbol ◦ . Nevertheless, ϛα1 (n)α2 (m) are the weighting coefficient matrices assessed in Eqs. (79)–(80) for the derivatives of n-th and
m-th order calculated along α 1 and α 2 directions. All the terms Lij(k) i j occurring in Eq. (A.9) have been hereafter collected by row.
                                                                                       (τη)α α
→(τη)α1 α1  →o− 3 →(τη)[00]α1 α1 α1 (1)α2 (0) → →o− 2 →o− 1 →(τη)[00]α1 α1 α1 (1)α2 (0) →
L 11(4) = − A 1 ◦ A 11(20)      ◦ϛ            A 1 + A 1 ◦ A 2 ◦ A 11(20)  ◦ϛ            A2 +
                                                                                 (                             )
  →o− 2             →(τη)[00]α1 α1 →o− 1 →o− 1 α1 (0)α2 (1) →(τη)[00]α1 α1 →o− 1 →(τη)[01]α1 α1 →(τη)[10]α1 α1
+ A 1 ◦ϛα1 (1)α2 (0) A 11(20)     + A 1 ◦ A 2 ◦ϛ            A 16(11)      + A 1 ◦ A 14(10)     − A 14(10)
→(τη)α1 α1  →o− 3 →(τη)[00]α1 α1 α1 (0)α2 (1) → →o− 1 →o− 2 →(τη)[00]α1 α1 α1 (0)α2 (1) →
L 11(5) = − A 2 ◦ A 66(02)      ◦ϛ            A 2 + A 1 ◦ A 2 ◦ A 66(02)  ◦ϛ            A1 +
                                                                                 (                             )
  →o− 2             →(τη)[00]α1 α1 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α1 α1 →o− 1 →(τη)[01]α1 α1 →(τη)[10]α1 α1
+ A 2 ◦ϛα1 (0)α2 (1) A 66(02)     + A 1 ◦ A 2 ◦ϛ            A 16(11)      + A 2 ◦ A 46(01)     − A 46(01)
→(τη)α1 α1 →o− 2 →o− 1 →(τη)[00]α1 α1 ( α1 (2)α2 (0) →     →o− 1             →                  → ) →o− 2 →o− 1 →(τη)[00]α1 α1 (→o− 1 α1 (1)α2 (0) → α1 (0)α2 (1) →
L 11(6) = A 1 ◦ A 2 ◦ A 12(11)       ◦ ϛ             A 2 − A 1 ◦ϛα1 (1)α2 (0) A 1 ◦ϛα1 (1)α2 (0) A 2 + A 1 ◦ A 2 ◦ A 16(20)   ◦ A 1 ◦ϛ             A 1 ◦ϛ         A1
                              )
               α1 (1)α2 (1) →
           − ϛ              A1 +
             (                                                                         (                               )
  →o− 1 →o− 1 →o− 1             →(τη)[00]α1 α1 →o− 1 α1 (0)α2 (1) →(τη)[00]α1 α1 →o− 1   →(τη)[00]α1 α1 →(τη)[00]α1 α1
+ A 1 ◦ A 2 ◦ A 1 ◦ϛα1 (1)α2 (0) A 12(11)     + A 2 ◦ϛ            A 26(02)      + R 1 ◦ 2 A 24(11)     − A 14(20)        +
                                                   )
  →(τη)[01]α1 α1 →(τη)[01]α1 α1 →(τη)[10]α1 α1                       →
+ A 14(10)      − A 24(01)     − A 24(01)              ◦ϛα1 (1)α2 (0) A 2 +
             (                                                                         (                               )
  →o− 1 →o− 1 →o− 1             →(τη)[00]α1 α1 →o− 1 α1 (0)α2 (1) →(τη)[00]α1 α1 →o− 1   →(τη)[00]α1 α1 →(τη)[00]α1 α1
− A 1 ◦ A 2 ◦ A 1 ◦ϛα1 (1)α2 (0) A 16(20)     + A 2 ◦ϛ            A 66(11)      + R 1 ◦ 2 A 46(20)     + A 46(11)        +
                                                   )
  →(τη)[01]α1 α1 →(τη)[01]α1 α1 →(τη)[10]α1 α1                       →       →o− 2 →o− 2 →(τη)[00]α1 α1 α1 (0)α2 (1) → α1 (1)α2 (0) →
− A 46(10)      − A 46(01)     − A 46(10)              ◦ϛα1 (0)α2 (1) A 1 + 2 A 1 ◦ A 2 ◦ A 26(11)     ◦ϛ            A 1 ◦ϛ         A2 +
             (                                                                        )                              (                             )
  →o− 2 →o− 2 →(τη)[00]α1 α1 ( α1 (1)α2 (0) → )o2 →(τη)[00]α1 α1 ( α1 (0)α2 (1) → )o2     →o− 2 →(τη)[00]α1 α1 →o− 1 →(τη)[01]α1 α1 →(τη)[10]α1 α1     →(τη)[11]α1 α1
− A 1 ◦ A 2 ◦ A 22(02)      ◦ ϛ             A2   + A 66(20)     ◦ ϛ             A1      − R 1 ◦ A 44(20)      + R 1 ◦ A 44(10)     + A 44(10)        − A 44(00)       +
  →o− 1 →o− 1             →(τη)[00]α1 α1 →o− 1 α1 (1)α2 (0) →(τη)[01]α1 α1 →o− 1 →o− 1 α1 (0)α2 (1) →(τη)[00]α1 α1
− A 1 ◦ R 1 ◦ϛα1 (1)α2 (0) A 14(20)     + A 1 ◦ϛ            A 14(10)      − A 2 ◦ R 1 ◦ϛ            A 46(11)       +
  →o− 1             →(τη)[01]α1 α1 →o− 1 →o− 2 →(τη)[00]α1 α1 α1 (1)α2 (0) → →o− 1 →o− 2 →(τη)[00]α1 α1 α1 (0)α2 (1) →
+ A 2 ◦ϛα1 (0)α2 (1) A 46(01)     + A 1 ◦ R 1 ◦ A 14(20)     ◦ϛ            R 1 + A 2 ◦ R 1 ◦ A 46(11)  ◦ϛ            R1                                                  (A.11)
                                                               (                          )
→(τη)α1 α2  →o− 3 →(τη)[00]α1 α2 α1 (1)α2 (0) → →o− 2 →o− 1 →(τη)[00]α1 α2 →(τη)[00]α1 α2                 →
L 12(4) = − A 1 ◦ A 16(20)      ◦ϛ            A 1 + A 1 ◦ A 2 ◦ A 11(20)  + A 66(20)        ◦ϛα1 (0)α2 (1) A 1 +
                                                                                                 45
F. Tornabene et al.                                                                                                                                    Composite Structures 309 (2023) 116542
             (                                             )
  →o− 2 →o− 1 →(τη)[00]α1 α2 →(τη)[00]α1 α2 →(τη)[00]α1 α2                 →
+ A 1 ◦ A 2 ◦ A 16(20)      − A 16(11)     − A 26(11)        ◦ϛα1 (1)α2 (0) A 2 +
             (                                             )
  →o− 1 →o− 2 →(τη)[00]α1 α2 →(τη)[00]α1 α2 →(τη)[00]α1 α2                 →
+ A 1 ◦ A 2 ◦ A 26(02)      + A 16(11)     + A 26(11)        ◦ϛα1 (0)α2 (1) A 1 +
  →o− 2             →(τη)[00]α1 α2 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α1 α2 →o− 1 →(τη)[10]α1 α2 →o− 1 →o− 1 →(τη)[00]α1 α2
+ A 2 ◦ϛα1 (0)α2 (1) A 26(02)     + A 1 ◦ A 2 ◦ϛ            A 12(11)      − A 2 ◦ A 24(01)     + A 2 ◦ R 1 ◦ A 24(11)      +
             (                                                                                              (                             )
  →o− 1 →o− 1 →o− 1             →(τη)[00]α1 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α1 α2 →o− 1 →(τη)[00]α1 α2 →o− 1 →(τη)[00]α1 α2 →(τη)[00]α1 α2
+ A 1 ◦ A 2 ◦ A 1 ◦ϛα1 (1)α2 (0) A 11(20)     + A 2 ◦ϛ            A 16(11)      + R 1 ◦ A 14(20)     − R 2 ◦ A 56(11)     + A 56(02)        +
                                                     )                                   (
  →(τη)[01]α1 α2 →(τη)[01]α1 α2 →(τη)[10]α1 α2                         → →o− 1 →o− 1 →o− 1                  →(τη)[00]α1 α2   →o− 1             →(τη)[00]α1 α2
+ A 56(10)      + A 56(01)     − A 14(10)                ◦ϛα1 (0)α2 (1) A 1 − A 1 ◦ A 2 ◦ A 1 ◦ϛα1 (1)α2 (0) A 16(11)      + A 2 ◦ϛα1 (0)α2 (1) A 66(02)      +
                             (                             )                                                )
  →o− 1 →(τη)[00]α1 α2 →o− 1 →(τη)[00]α1 α2 →(τη)[00]α1 α2     →(τη)[01]α1 α2 →(τη)[01]α1 α2 →(τη)[10]α1 α2                 →
+ R 1 ◦ A 46(11)      + R 2 ◦ A 15(11)     − A 25(02)        + A 25(01)      − A 15(10)     − A 46(01)        ◦ϛα1 (1)α2 (0) A 2 +
             (
  →o− 2 →o− 2 →(τη)[00]α1 α2 ( α1 (1)α2 (0) → )o2 →(τη)[00]α1 α2 ( α1 (0)α2 (1) → )o2
+ A 1 ◦ A 2 ◦ A 26(02)      ◦ ϛ             A2   + A 16(20)     ◦ ϛ             A1    +
    (                                   )                                     )
        →(τη)[00]α1 α2 →(τη)[00]α1 α2                  →                  →      →o− 1 →(τη)[01]α1 α2 →o− 1 →(τη)[10]α1 α2 →o− 1 →o− 1 →(τη)[00]α1 α2
−       A 12(11)      + A 66(11)         ◦ϛα1 (0)α2 (1) A 1 ◦ϛα1 (1)α2 (0) A 2 + R 1 ◦ A 45(10)      + R 2 ◦ A 45(01)     − R 1 ◦ R 2 ◦ A 45(11)      +
                       (                                                               )        (                                                               )
  →(τη)[11]α1 α2 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α1 α2                →(τη)[01]α1 α2     →o− 1 →o− 1             →(τη)[00]α1 α2                →(τη)[01]α1 α2
− A 45(00)      − A 1 ◦ R 2 ◦ϛ            A 15(11)       − ϛα1 (1)α2 (0) A 15(10)        − A 2 ◦ R 2 ◦ϛα1 (0)α2 (1) A 56(02)      − ϛα1 (0)α2 (1) A 56(01)        +
       (                                                                             )
  →o− 2 →o− 1 →(τη)[00]α1 α2 α1 (1)α2 (0) →     →o− 1 →(τη)[00]α1 α2 α1 (0)α2 (1) →
+ R 2 ◦ A 1 ◦ A 15(11)      ◦ϛ            R 2 + A 2 ◦ A 56(02)      ◦ϛ            R2                                                                                                 (A.12)
                 (         (                               )
→(τη)α1 α3 →o− 1 →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3            →o− 1 →(τη)[00]α1 α3 →(τη)[01]α1 α3 →(τη)[10]α1 α3 →o− 2 →(τη)[00]α1 α3 α1 (1)α2 (0) →
L 13(4) = A 1 ◦ R 1 ◦ A 11(20)                + A 44(20)      + R 2 ◦ A 12(11)      + A 13(10)     − A 44(10)      − A 1 ◦ A 14(20)       ◦ϛ             A1
                                                                            (                             )
                →o− 1 →o− 1 →(τη)[00]α1 α3 α1 (0)α2 (1) → →o− 1 →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3                       → →o− 1                  →(τη)[00]α1 α3
           + + A 1 ◦ A 2 ◦ A 46(20)            ◦ϛ        A 1 + A 1 ◦ A 2 ◦ A 14(20)        − A 24(11)       ◦ϛα1 (1)α2 (0) A 2 + A 1 ◦ϛα1 (1)α2 (0) A 14(20)      +
                                              )
             →o− 1             →(τη)[00]α1 α3
           + A 2 ◦ϛα1 (0)α2 (1) A 46(11)
                                                                                                 46
F. Tornabene et al.                                                                                                                               Composite Structures 309 (2023) 116542
                (     (                              )
→(τη)α1 α3 →o− 1 →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3
L 13(5) = A 2 ◦ R 1 ◦ A 16(11)       + A 45(11)
                                                                                                                                                                            (
              →o− 1 →(τη)[00]α1 α3 →(τη)[01]α1 α3 →(τη)[10]α1 α3 →o− 2 →(τη)[00]α1 α3 α1 (0)α2 (1) →               →o− 1 →o− 1 →(τη)[00]α1 α3 α1 (1)α2 (0) → →o− 1 →o− 1 →(τη)[00]α1 α3
            + R 2 ◦ A 26(02)         + A 36(01)       − A 45(01)      − A 2 ◦ A 56(02)      ◦ϛ          A 2 + − A 1 ◦ A 2 ◦ A 25(02)         ◦ϛ            A 2 + A 1 ◦ A 2 ◦ A 56(11)
                             )                                                                                           )
              →(τη)[00]α1 α3                 →        →o− 1             →(τη)[00]α1 α3 →o− 1 α1 (1)α2 (0) →(τη)[00]α1 α3
            + A 56(02)         ◦ϛα1 (0)α2 (1) A 1 + + A 2 ◦ϛα1 (0)α2 (1) A 56(02)     + A 1 ◦ϛ            A 15(11)
                      (     (                              )        (                              )                                 )
→(τη)α1 α3 →o− 1 →o− 1 →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3     →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3     →(τη)[01]α1 α3 →(τη)[01]α1 α3                 →
L 13(6) = A 1 ◦ A 2 ◦ R 2 ◦ A 12(11)       − A 22(02)        + R 1 ◦ A 11(20)      − A 12(11)        + A 13(10)      − A 23(01)        ◦ϛα1 (1)α2 (0) A 2 +
                  (          (                          )               (                             )
  →o− 1 →o− 1 →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3               →o− 1 →(τη)[00]α1 α3 →(τη)[00]α1 α3          →(τη)[01]α1 α3
+ A 1 ◦ A 2 ◦ R 2 ◦ A 26(02)                 + A 26(11)    + R 1 ◦ A 16(11)          + A 16(20)           + A 36(01)
                   )                           (                                                                   )        (
    →(τη)[01]α1 α3                 → →o− 1 →o− 1                  →(τη)[00]α1 α3 →o− 1 α1 (1)α2 (0) →(τη)[00]α1 α3     →o− 1 →o− 1             →(τη)[00]α1 α3
  + A 36(10)         ◦ϛα1 (0)α2 (1) A 1 + A 1 ◦ R 1 ◦ϛα1 (1)α2 (0) A 11(20)     + R 2 ◦ϛ            A 12(11)         + A 2 ◦ R 1 ◦ϛα1 (0)α2 (1) A 16(11)      +
                                     )                                                                          (
  →o− 1             →(τη)[00]α1 α3     →o− 1            →(τη)[01]α1 α3 →o− 1 α1 (0)α2 (1) →(τη)[01]α1 α3 →o− 2 →(τη)[00]α1 α3 →o− 1 →(τη)[00]α1 α3 α1 (1)α2 (0) →
+ R 2 ◦ϛα1 (0)α2 (1) A 26(02)       + A 1 ◦ϛα1 (1)α2 (0) A 13(10)     + A 2 ◦ϛ            A 36(01)      + R 1 ◦ A 14(20)           − A 1 ◦ A 11(20) ◦ϛ          R1
                                         )               (                                                                             )
    →o− 1 →(τη)[00]α1 α3 α1 (0)α2 (1) →       →o− 2 →o− 1 →(τη)[00]α1 α3 α1 (1)α2 (0) →            →o− 1 →(τη)[00]α1 α3 α1 (0)α2 (1) →
  − A 2 ◦ A 16(11)      ◦ϛ            R 1 + − R 2 ◦ A 1 ◦ A 12(11)           ◦ϛ           R 2 + A 2 ◦ A 26(02)         ◦ϛ            R2 +
                                   (                             )
  →o− 1 →o− 1 →(τη)[00]α1 α3 →o− 1 →(τη)[01]α1 α3 →(τη)[10]α1 α3     →o− 1 →(τη)[10]α1 α3 →(τη)[11]α1 α3
+ R 1 ◦ R 2 ◦ A 24(11)      + R 1 ◦ A 34(10)     − A 14(10)        − R 2 ◦ A 24(01)      − A 34(00)                                                                             (A.13)
                                                               (                          )
→(τη)α2 α1  →o− 3 →(τη)[00]α2 α1 α1 (1)α2 (0) → →o− 2 →o− 1 →(τη)[00]α2 α1 →(τη)[00]α2 α1                 →
L 21(4) = − A 1 ◦ A 16(20)      ◦ϛ            A 1 − A 1 ◦ A 2 ◦ A 11(20)  + A 66(20)        ◦ϛα1 (0)α2 (1) A 1 +
             (                                             )
  →o− 2 →o− 1 →(τη)[00]α2 α1 →(τη)[00]α2 α1 →(τη)[00]α2 α1                 →
+ A 1 ◦ A 2 ◦ A 16(20)      + A 16(11)     + A 26(11)        ◦ϛα1 (1)α2 (0) A 2 +
                                                                                 (                             )
  →o− 2             →(τη)[00]α2 α1 →o− 1 →o− 1 α1 (0)α2 (1) →(τη)[00]α2 α1 →o− 1 →(τη)[01]α2 α1 →(τη)[10]α2 α1
+ A 1 ◦ϛα1 (1)α2 (0) A 16(20)     + A 1 ◦ A 2 ◦ϛ            A 12(11)      + A 1 ◦ A 46(10)     − A 15(10)        +
       (                                          )
  →o− 1 →o− 1 →(τη)[00]α2 α1 →o− 1 →(τη)[00]α2 α1
+ A 1 ◦ R 2 ◦ A 15(11)      − R 1 ◦ A 46(20)
                                                               (                          )
→(τη)α2 α1  →o− 3 →(τη)[00]α2 α1 α1 (0)α2 (1) → →o− 1 →o− 2 →(τη)[00]α2 α1 →(τη)[00]α2 α1                 →
L 21(5) = − A 2 ◦ A 26(02)      ◦ϛ            A 2 + A 1 ◦ A 2 ◦ A 22(02)  + A 66(02)        ◦ϛα1 (1)α2 (0) A 2 +
             (                                             )
  →o− 1 →o− 2 →(τη)[00]α2 α1 →(τη)[00]α2 α1 →(τη)[00]α2 α1                 →
+ A 1 ◦ A 2 ◦ A 26(02)      − A 16(11)     − A 26(11)        ◦ϛα1 (0)α2 (1) A 1 +
                                                                                 (                             )
  →o− 2             →(τη)[00]α2 α1 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α1 →o− 1 →(τη)[01]α2 α1 →(τη)[10]α2 α1
+ A 2 ◦ϛα1 (0)α2 (1) A 26(02)     + A 1 ◦ A 2 ◦ϛ            A 66(11)      + A 2 ◦ A 24(01)     − A 56(01)        +
             (
  →o− 1 →o− 1 →(τη)[01]α2 α1 →(τη)[01]α2 α1 →(τη)[10]α2 α1 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α1
+ A 1 ◦ A 2 ◦ A 46(10)      + A 46(01)     − A 25(01)     + A 1 ◦ϛ            A 26(11)       +
                                                              (                             ))
  →o− 1             →(τη)[00]α2 α1 →o− 1 →(τη)[00]α2 α1 →o− 1 →(τη)[00]α2 α1 →(τη)[00]α2 α1                  →
+ A 2 ◦ϛα1 (0)α2 (1) A 22(02)     + R 2 ◦ A 25(02)     − R 1 ◦ A 46(20)     + A 46(11)         ◦ϛα1 (1)α2 (0) A 2 +
                                                                                            47
F. Tornabene et al.                                                                                                                                Composite Structures 309 (2023) 116542
              (
   →o− 1 →o− 1 →(τη)[01]α2 α1 →(τη)[01]α2 α1 →(τη)[10]α2 α1 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α1
 − A 1 ◦ A 2 ◦ A 14(10)      − A 24(01)     − A 56(10)     + A 1 ◦ϛ            A 66(20)       +
                                                               (                             ))
   →o− 1             →(τη)[00]α2 α1 →o− 1 →(τη)[00]α2 α1 →o− 1 →(τη)[00]α2 α1 →(τη)[00]α2 α1                  →
 + A 2 ◦ϛα1 (0)α2 (1) A 26(11)     + R 2 ◦ A 56(11)     + R 1 ◦ A 24(11)     − A 14(20)         ◦ϛα1 (0)α2 (1) A 1 +
              (                                                                          (                               )                                      )
   →o− 2 →o− 2 →(τη)[00]α2 α1 ( α1 (1)α2 (0) → )o2 →(τη)[00]α2 α1 ( α1 (0)α2 (1) → )o2     →(τη)[00]α2 α1 →(τη)[00]α2 α1                 →                  →
 + A 1 ◦ A 2 ◦ A 26(02)      ◦ ϛ             A2   + A 16(20)     ◦ ϛ             A1    −   A 12(11)      + A 66(11)        ◦ϛα1 (0)α2 (1) A 1 ◦ϛα1 (1)α2 (0) A 2 +
(A.14)
→(τη)α2 α2  →o− 3 →(τη)[00]α2 α2 α1 (1)α2 (0) → →o− 2 →o− 1 →(τη)[00]α2 α2 α1 (1)α2 (0) →
L 22(4) = − A 1 ◦ A 66(20)      ◦ϛ            A 1 + A 1 ◦ A 2 ◦ A 66(20)  ◦ϛ            A2 +
        (                                                                                                  )
   →o− 1 →(τη)[01]α2 α2 →(τη)[10]α2 α2 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α2 α2
 + A 1 ◦ A 56(10)      − A 56(10)     + A 1 ◦ϛ            A 66(20)      + A 2 ◦ϛ            A 26(11)
                     (
    τη α α2     →o− 1 →(τη)[01]α2 α2 →(τη)[10]α2 α2 →o− 2 →(τη)[00]α2 α2 α1 (0)α2 (1) →
¯L(22(5)
      ) 2
              = A 2 ◦ A 25(01)      − A 25(01)     − A 2 ◦ A 22(02)     ◦ϛ            A2 +
                                                                                                                      )
   →o− 1 →o− 1 →(τη)[00]α2 α2 α1 (0)α2 (1) →     →o− 1             →(τη)[00]α2 α2 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α2
 + A 1 ◦ A 2 ◦ A 22(02)      ◦ϛ            A 1 + A 2 ◦ϛα1 (0)α2 (1) A 22(02)     + A 1 ◦ϛ            A 26(11)
   (τη)α α2       →o− 2 →o− 1 →(τη)[00]α2 α2 ( α1 (2)α2 (0) →     →o− 1             →                  → )      →o− 2 →o− 1 →(τη)[00]α2 α2 (→o− 1 α1 (1)α2 (0) → α1 (0)α2 (1) →
¯L22(6)2      = − A 1 ◦ A 2 ◦ A 66(11)      ◦ ϛ             A 2 − A 1 ◦ϛα1 (1)α2 (0) A 1 ◦ϛα1 (1)α2 (0) A 2 + − A 1 ◦ A 2 ◦ A 16(20)      ◦ A 1 ◦ϛ             A 1 ◦ϛ         A1
                                   )
                    α1 (1)α2 (1) →
                − ϛ              A1 +
              (      (                               )
   →o− 1 →o− 1 →o− 1   →(τη)[00]α2 α2 →(τη)[00]α2 α2     →o− 1             →(τη)[00]α2 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α2 α2
 − A 1 ◦ A 2 ◦ R 2 ◦ 2 A 56(02)      + A 56(11)        + A 1 ◦ϛα1 (1)α2 (0) A 66(11)     + A 2 ◦ϛ            A 26(02)       +
                                                   )
   →(τη)[01]α2 α2 →(τη)[01]α2 α2 →(τη)[10]α2 α2                      →
 − A 56(10)      − A 56(01)     − A 56(01)             ◦ϛα1 (1)α2 (0) A 2 +
              (      (                               )
   →o− 1 →o− 1 →o− 1   →(τη)[00]α2 α2 →(τη)[00]α2 α2     →o− 1             →(τη)[00]α2 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α2 α2
 + A 1 ◦ A 2 ◦ R 2 ◦ 2 A 15(11)      − A 25(02)        + A 1 ◦ϛα1 (1)α2 (0) A 16(20)     + A 2 ◦ϛ            A 12(11)       +
                                                   )
   →(τη)[01]α2 α2 →(τη)[01]α2 α2 →(τη)[10]α2 α2                   →
 + A 25(01)      − A 15(10)     − A 15(10)          ◦ϛα1 (0)α2 (1) A 1 +
                   (
   →o− 2 →o− 2        →(τη)[00]α2 α2 α1 (0)α2 (1) → α1 (1)α2 (0) →     →(τη)[00]α2 α2 ( α1 (1)α2 (0) → )o2
 + A 1 ◦ A 2 ◦ 2 A 16(11)             ◦ϛ          A 1 ◦ϛ         A 2 − A 66(02)      ◦ ϛ             A2
                     (                  ) )
     →(τη)[00]α2 α2                → o2             →o− 1 →o− 1             →(τη)[00]α2 α2 →o− 1 α1 (0)α2 (1) →(τη)[01]α2 α2 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α2
   − A 11(20)       ◦ ϛα1 (0)α2 (1) A 1      + − A 2 ◦ R 2 ◦ϛα1 (0)α2 (1) A 25(02)        + A 2 ◦ϛ            A 25(01)      − A 1 ◦ R 2 ◦ϛ            A 56(11)       +
   →o− 1             →(τη)[01]α2 α2 →o− 1 →o− 2 →(τη)[00]α2 α2 α1 (0)α2 (1) → →o− 1 →o− 2 →(τη)[00]α2 α2 α1 (1)α2 (0) →
 + A 1 ◦ϛα1 (1)α2 (0) A 56(10)     + A 2 ◦ R 2 ◦ A 25(02)     ◦ϛ            R 2 + A 1 ◦ R 2 ◦ A 56(11)  ◦ϛ            R2 +
                              (                             )
   →o− 2 →(τη)[00]α2 α2 →o− 1 →(τη)[01]α2 α2 →(τη)[10]α2 α2     →(τη)[11]α2 α2
 − R 2 ◦ A 55(02)      + R 2 ◦ A 55(01)     + A 55(01)        − A 55(00)                                                                                                           (A.15)
                                                                                             48
F. Tornabene et al.                                                                                                                        Composite Structures 309 (2023) 116542
                                                                (                                    (                               )                   )       (
→(τη)α2 α3   →o− 3 →(τη)[00]α2 α3 α1 (1)α2 (0) → →o− 2 →o− 1        →(τη)[00]α2 α3 α1 (0)α2 (1) →      →(τη)[00]α2 α3 →(τη)[00]α2 α3                 →      →o− 1 →(τη)[01]α2 α3
L 23(4) = − A 1 ◦ A 46(20)       ◦ϛ            A 1 + A 1 ◦ A 2 ◦ − A 14(20)       ◦ϛ            A 1 + A 46(20)       + A 46(11)        ◦ϛα1 (1)α2 (0) A 2 + A 1 ◦ A 36(10)
                                                           (                            ))
             →(τη)[10]α2 α3 →o− 1 →(τη)[00]α2 α3 →o− 1 →(τη)[00]α2 α3 →(τη)[00]α2 α3
           − A 45(10)      + R 1 ◦ A 16(20)        + R 2 ◦ A 26(11)  + A 45(11)                 +
             (     (                              )        (                              )                                 )
  →o− 1 →o− 1 →o− 1 →(τη)[00]α2 α3 →(τη)[00]α2 α3     →o− 1 →(τη)[00]α2 α3 →(τη)[00]α2 α3     →(τη)[01]α2 α3 →(τη)[01]α2 α3                →
+ A 1 ◦ A 2 ◦ R 1 ◦ A 12(11)      − A 11(02)        + R 2 ◦ A 22(02)      − A 12(11)        + A 23(01)      − A 13(10)       ◦ϛα1 (0)α2 (1) A 1 +
       (                                                                                                 )
  →o− 1 →o− 1             →(τη)[00]α2 α3 →o− 1 α1 (1)α2 (0) →(τη)[00]α2 α3                →(τη)[01]α2 α3
+ A 1 ◦ R 1 ◦ϛα1 (1)α2 (0) A 16(20)     + R 2 ◦ϛ            A 26(11)       + ϛα1 (1)α2 (0) A 36(10)        +
       (                                                                                                 )
  →o− 1 →o− 1             →(τη)[00]α2 α3 →o− 1 α1 (0)α2 (1) →(τη)[00]α2 α3                →(τη)[01]α2 α3
+ A 2 ◦ R 1 ◦ϛα1 (0)α2 (1) A 12(11)     + R 2 ◦ϛ            A 22(02)       + ϛα1 (0)α2 (1) A 23(01)        +
        (                                                                        )         (
  →o− 1 →o− 2 →(τη)[00]α2 α3 α1 (1)α2 (0) → →o− 2 →(τη)[00]α2 α3 α1 (1)α2 (0) →      →o− 1 →o− 2 →(τη)[00]α2 α3 α1 (0)α2 (1) →
− A 1 ◦ R 1 ◦ A 16(20)         ◦ϛ         R 1 + R 2 ◦ A 26(11)    ◦ϛ          R 2 − A 2 ◦ R 1 ◦ A 12(11)        ◦ϛ           R1 +
                                        )
   →o− 2 →(τη)[00]α2 α3 α1 (0)α2 (1) →    →o− 1 →o− 1 →(τη)[00]α2 α3 →o− 2 →(τη)[00]α2 α3 →o− 1 →(τη)[01]α2 α3     →o− 1 →(τη)[10]α2 α3 →o− 1 →(τη)[10]α2 α3 →(τη)[11]α2 α3
 + R 2 ◦ A 22(02)      ◦ϛ            R 2 + R 1 ◦ R 2 ◦ A 15(11)     + R 2 ◦ A 25(02)     + R 2 ◦ A 35(01)      + − R 1 ◦ A 15(10)      − R 2 ◦ A 25(01)     − A 35(00)
                                                                                                                                                                          (A.16)
    Third Row of the fundamental operator
→(τη)α3 α1 →o− 2 →(τη)[00]α3 α1
L 31(1) = A 1 ◦ A 14(20)
                                                                ((                               )                                                      )       (
→(τη)α3 α1   →o− 3 →(τη)[00]α3 α1 α1 (1)α2 (0) → →o− 2 →o− 1       →(τη)[00]α3 α1 →(τη)[00]α3 α1                →        →(τη)[00]α3 α1 α1 (0)α2 (1) →     →o− 1 →(τη)[10]α3 α1
L 31(4) = − A 1 ◦ A 14(20)       ◦ϛ            A1 + A1 ◦A2 ◦       A 14(20)      + A 24(11)       ◦ϛα1 (1)α2 (0) A 2 + − A 46(20)      ◦ϛ            A 1 − A 1 ◦ A 13(10)
                                    (                           )                         )
             →(τη)[01]α3 α1 →o− 1 →(τη)[00]α3 α1 →(τη)[00]α3 α1     →o− 1 →(τη)[00]α3 α1
           − A 44(10)      + R 1 ◦ A 11(20)         + A 44(20)    + R 2 ◦ A 12(11)          +
                                                     (                               )               (
  →o− 1 →o− 2 →(τη)[00]α3 α1 α1 (1)α2 (0) → →o− 1 →(τη)[10]α3 α1 →(τη)[01]α3 α1          →o− 1 →o− 1 →(τη)[00]α3 α1
+ A 1 ◦ A 2 ◦ A 25(02)      ◦ϛ            A 2 − A 2 ◦ A 36(01)   − A 45(01)            − A 2 ◦ R 1 ◦ A 16(11)
                   )
    →(τη)[00]α3 α1       →o− 1 →o− 1 →(τη)[00]α3 α1 →o− 2 α1 (0)α2 (1) →(τη)[00]α3 α1 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α1
  + A 45(11)         + − A 2 ◦ R 2 ◦ A 26(02)       + A 2 ◦ϛ           A 56(02)       + A 1 ◦ A 2 ◦ϛ           A 46(11)
                                                                                           49
F. Tornabene et al.                                                                                                                                    Composite Structures 309 (2023) 116542
                                 (          (                                 )
→(τη)α3 α1       →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α1 →(τη)[00]α3 α1                    →o− 1 →(τη)[00]α3 α1 →(τη)[01]α3 α1 →(τη)[10]α3 α1 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α1
L 31(6)      = − A 1 ◦ A 2 ◦ R 1 ◦ A 44(20)                  + A 12(11)           + R 2 ◦ A 22(02)        + A 44(10)    − A 23(01)    − A 1 ◦ϛ               A 24(11)      +
                                                  )                                    (           (                           )
                 →o− 1             →(τη)[00]α3 α1                 → →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α1 →(τη)[00]α3 α1                          →(τη)[00]α3 α1       →(τη)[01]α3 α1 →(τη)[10]α3 α1
               − A 2 ◦ϛα1 (0)α2 (1) A 25(02)        ◦ϛα1 (1)α2 (0) A 2 + A 1 ◦ A 2 ◦ R 1 ◦ A 16(20)              − A 45(11)        + Ro−2 1 ◦ A 26(11)      + + A 45(01)        + A 36(10)
                                                                                     )                                   (
                 →o− 1             →(τη)[00]α3 α1 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α1                  → →o− 1 →o− 1 →o− 2 →(τη)[00]α3 α1 α1 (1)α2 (0) → α1 (0)α2 (1) →
               − A 1 ◦ϛα1 (1)α2 (0) A 46(20)      − A 2 ◦ϛ               A 56(11)      ◦ϛα1 (0)α2 (1) A 1 + A 1 ◦ A 2 ◦ A 1 ◦ A 46(20)         ◦ϛ           A 1 ◦ϛ          A1
                                                                                                                          )
  →o− 1 →(τη)[00]α3 α1 α1 (1)α2 (1) →     →o− 1 →(τη)[00]α3 α1 α1 (2)α2 (0) →     →o− 1 →(τη)[00]α3 α1 α1 (0)α2 (2) →
+ A 2 ◦ A 25(02)      ◦ϛ            A 2 + A 1 ◦ A 24(11)      ◦ϛ            A 2 − A 2 ◦ A 56(11)      ◦ϛ            A1        +
          (                                                                     )         (
  →o− 1           →o− 1            →(τη)[00]α3 α1                →(τη)[01]α3 α1     →o− 1     →o− 1              →(τη)[00]α3 α1
+ A 1 ◦ − R 1 ◦ϛα1 (1)α2 (0) A 44(20)             + ϛα1 (1)α2 (0) A 44(10)        + A 2 ◦ − R 1 ◦ϛα1 (0)α2 (1) A 45(11)
                                )          (                                                                                                   )
                 →(τη)[01]α3 α1     →o− 2 →o− 1 →(τη)[00]α3 α1 α1 (1)α2 (0) →           →o− 1 →(τη)[00]α3 α1 α1 (0)α2 (1) →     →(τη)[00]α3 α1
  + ϛα1 (0)α2 (1) A 45(01)        + R 1 ◦ A 1 ◦ A 44(20)              ◦ϛ          R 1 + A 2 ◦ A 45(11)      ◦ϛ            R 1 + A 14(20)         +
  →o− 1 →(τη)[01]α3 α1 →o− 1 →o− 1 →(τη)[00]α3 α1 →o− 1 →(τη)[10]α3 α1 →o− 1 →(τη)[01]α3 α1 →(τη)[11]α3 α1
− R 1 ◦ A 14(10)      + R 1 ◦ R 2 ◦ A 24(11)     + R 1 ◦ A 34(10)     − R 2 ◦ A 24(01)     − A 34(00)                                                                                  (A.17)
                 (
→(τη)α3 α2  →o− 1 →o− 2 →(τη)[00]α3 α2 α1 (1)α2 (0) →     →o− 1 →o− 1 →(τη)[00]α3 α2 α1 (0)α2 (1) →
L 32(4) = − A 1 ◦ A 1 ◦ A 46(20)      ◦ϛ            A 1 − A 1 ◦ A 2 ◦ A 14(20)      ◦ϛ            A1 +
             (                              )
  →o− 1 →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2                → →(τη)[10]α3 α2 →(τη)[01]α3 α2 →o− 1 →(τη)[00]α3 α2
+ A 1 ◦ A 2 ◦ A 46(11)      − A 46(20)       ◦ϛα1 (1)α2 (0) A 2 + A 36(10) − A 45(10)     + R 1 ◦ A 16(20)      +
       (                              )                                                                      )
  →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2     →o− 1             →(τη)[00]α3 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α2
+ R 2 ◦ A 26(11)      + A 45(11)        − A 1 ◦ϛα1 (1)α2 (0) A 46(20)     − A 2 ◦ϛ            A 56(11)
                 (
→(τη)α3 α2  →o− 1 →o− 2 →(τη)[00]α3 α2 α1 (0)α2 (1) →     →o− 1 →o− 1 →(τη)[00]α3 α2 α1 (1)α2 (0) →
L 32(5) = − A 2 ◦ A 2 ◦ A 25(02)      ◦ϛ            A 2 + A 1 ◦ A 2 ◦ A 56(02)      ◦ϛ            A2 +
               (                             )                                                                           (                             )
  →o− 1 →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2                   → →(τη)[10]α3 α2 →(τη)[01]α3 α2 →o− 1 →(τη)[00]α3 α2 →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2
− A 1 ◦ A 2 ◦ A 25(02)          + A 15(11)     ◦ϛα1 (0)α2 (1) A 1 + A 23(01)  − A 55(01)     + R 1 ◦ A 12(11)     + R 2 ◦ A 22(02)     + A 55(02)
                                                                       )
    →o− 1             →(τη)[00]α3 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α2
  − A 1 ◦ϛα1 (1)α2 (0) A 24(11)     − A 2 ◦ϛ            A 25(02)
                      (                           (                             )
→(τη)α3 α2 →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α2 →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2     →(τη)[01]α3 α2
L 32(6) = A 1 ◦ A 2 ◦ R 1 ◦ A 16(11)       + R 2 ◦ A 26(02)     − A 45(11)        + A 45(10)       +
                                                                                        )
  →(τη)[10]α3 α2 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α2 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α2                   →
+ A 36(01)      − A 1 ◦ϛ            A 46(11)      − A 2 ◦ϛ            A 56(02)           ◦ϛα1 (1)α2 (0) A 2 +
               (                                      (                )
  →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α2 →o− 1 →(τη)[00]α3 α2 →(τη)[00]α3 α2     →(τη)[10]α3 α2 →(τη)[01]α3 α2 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α2
− A 1 ◦ A 2 ◦ R 1 ◦ A 11(20)             + R 2 ◦ A 12(11)   + A 55(02)   − A 13(10)      − A 55(01)     − A 1 ◦ϛ            A 14(20)       +
                                     )
    →o− 1             →(τη)[00]α3 α2                 →
  − A 2 ◦ϛα1 (0)α2 (1) A 15(11)        ◦ϛα1 (0)α2 (1) A 1 +
              (
  →o− 1 →o− 1   →o− 2 →(τη)[00]α3 α2 α1 (1)α2 (0) → α1 (0)α2 (1) →     →o− 2 →(τη)[00]α3 α2 α1 (1)α2 (0) → α1 (1)α2 (0) →
+ A 1 ◦ A 2 ◦ − A 1 ◦ A 14(20)      ◦ϛ            A 1 ◦ϛ         A 1 + A 1 ◦ A 46(11)      ◦ϛ            A 1 ◦ϛ         A2
  →o− 2 →(τη)[00]α3 α2 α1 (0)α2 (1) → α1 (1)α2 (0) → →o− 1 →(τη)[00]α3 α2 α1 (1)α2 (1) → →o− 1 →(τη)[00]α3 α2 α1 (1)α2 (1) →       →o− 1 →(τη)[00]α3 α2 α1 (0)α2 (2) →
+ A 2 ◦ A 56(02)      ◦ϛ            A 2 ◦ϛ         A 2 + A 1 ◦ A 14(20)  ◦ϛ            A 1 − A 2 ◦ A 56(02)  ◦ϛ            A 2 + + A 2 ◦ A 15(11)      ◦ϛ            A1
                                           )
    →o− 1 →(τη)[00]α3 α2 α1 (2)α2 (0) →
  − A 1 ◦ A 46(11)       ◦ϛ            A2 +
                                                                                                50
F. Tornabene et al.                                                                                                                        Composite Structures 309 (2023) 116542
          (                                                                     )         (
  →o− 1           →o− 1            →(τη)[00]α3 α2                →(τη)[01]α3 α2     →o− 1     →o− 1              →(τη)[00]α3 α2
+ A 1 ◦ − R 2 ◦ϛα1 (1)α2 (0) A 45(11)             + ϛα1 (1)α2 (0) A 45(10)        + A 2 ◦ − R 2 ◦ϛα1 (0)α2 (1) A 55(02)         +
                                )          (                                                                                )
                 →(τη)[01]α3 α2     →o− 2 →o− 1 →(τη)[00]α3 α2 α1 (1)α2 (0) →           →o− 1 →(τη)[00]α3 α2 α1 (0)α2 (1) →
  + ϛα1 (0)α2 (1) A 55(01)        + R 2 ◦ A 1 ◦ A 45(11)              ◦ϛ          R 2 + A 2 ◦ A 55(02)      ◦ϛ            R2 +
  →o− 1 →o− 1 →(τη)[00]α3 α2 →o− 1 →(τη)[01]α3 α2 →o− 2 →(τη)[00]α3 α2 →o− 1 →(τη)[01]α3 α2 →o− 1 →(τη)[10]α3 α2 →(τη)[11]α3 α2
+ R 1 ◦ R 2 ◦ A 15(11)      − R 1 ◦ A 15(10)     + R 2 ◦ A 25(02)     − R 2 ◦ A 25(01)     + R 2 ◦ A 35(01)     − A 35(00)                                               (A.18)
                                                                           )
  →o− 1             →(τη)[00]α3 α3 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α3
+ A 1 ◦ϛα1 (1)α2 (0) A 44(20)     + A 2 ◦ϛ            A 45(11)
                 (
→(τη)α3 α3 →o− 1   →o− 2 →(τη)[00]α3 α3 α1 (0)α2 (1) →     →o− 1 →o− 1 →(τη)[00]α3 α3 α1 (0)α2 (1) →     →(τη)[01]α3 α3 →(τη)[10]α3 α3
L 33(5) = A 2 ◦ − A 2 ◦ A 55(02)       ◦ϛ            A 2 + A 1 ◦ A 2 ◦ A 55(02)      ◦ϛ            A 1 + A 35(01)      − A 35(01)      +
                                                                           )
  →o− 1             →(τη)[00]α3 α3 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α3
+ A 2 ◦ϛα1 (0)α2 (1) A 55(02)     + A 1 ◦ϛ            A 45(11)
                      (                                                         )
→(τη)α3 α3 →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α3 →o− 1 →(τη)[00]α3 α3 →(τη)[01]α3 α3                 →
L 33(6) = A 1 ◦ A 2 ◦ R 1 ◦ A 14(20)       + R 2 ◦ A 24(11)     + A 34(10)        ◦ϛα1 (1)α2 (0) A 2 +
             (                                                         )
  →o− 1 →o− 1 →o− 1 →(τη)[00]α3 α3 →o− 1 →(τη)[00]α3 α3 →(τη)[01]α3 α3                →
+ A 1 ◦ A 2 ◦ R 1 ◦ A 15(11)      + R 2 ◦ A 25(02)     + A 35(01)       ◦ϛα1 (0)α2 (1) A 1 +
                                                                           (
  →o− 1             →(τη)[01]α3 α3 →o− 1 α1 (0)α2 (1) →(τη)[01]α3 α3 →o− 1 →o− 1 α1 (1)α2 (0) →(τη)[00]α3 α3
+ A 1 ◦ϛα1 (1)α2 (0) A 34(10)     + A 2 ◦ϛ            A 35(01)      + A 1 ◦ R 1 ◦ϛ            A 14(20)       +
                                     )          (                                                                   )
  →o− 1             →(τη)[00]α3 α3         →o− 1 →o− 1             →(τη)[00]α3 α3 →o− 1 α1 (0)α2 (1) →(τη)[00]α3 α3
+ R 2 ◦ϛα1 (1)α2 (0) A 24(11)            + A 2 ◦ R 1 ◦ϛα1 (0)α2 (1) A 15(11)     + R 2 ◦ϛ            A 25(02)         +
        (                                                                             )         (
  →o− 2 →o− 1 →(τη)[00]α3 α3 α1 (1)α2 (0) →      →o− 1 →(τη)[00]α3 α3 α1 (0)α2 (1) →     →o− 2 →o− 1 →(τη)[00]α3 α3 α1 (1)α2 (0) →
− R 1 ◦ A 1 ◦ A 14(20)         ◦ϛ          R 1 + A 2 ◦ A 15(11)      ◦ϛ            R 1 − R 2 ◦ A 1 ◦ A 24(11)        ◦ϛ          R2 +
                                         )                                                                                      (
    →o− 1 →(τη)[00]α3 α3 α1 (0)α2 (1) →      →o− 2 →(τη)[00]α3 α3 →o− 2 →(τη)[00]α3 α3     →o− 1 →o− 1 →(τη)[00]α3 α3 →o− 1 →(τη)[01]α3 α3
  + A 2 ◦ A 25(02)      ◦ϛ            R 2 − R 1 ◦ A 11(20)       − R 2 ◦ A 22(02)       − 2 R 1 ◦ R 2 ◦ A 12(11)      − R 1 ◦ A 13(10)
                   )              (                             )
    →(τη)[10]α3 α3       →o− 1 →(τη)[01]α3 α3 →(τη)[10]α3 α3         →(τη)[11]α3 α3
  + A 13(10)         + − R 2 ◦ A 23(01)         + A 23(01)        − A 33(00)                                                                                             (A.19)
Appendix III
   We now present the discrete version of each component of the generalized operator accounted in Eq. (25) for each τ, η-th order of the field variable
kinematic expansion. They are expressed in terms of the generalized GDQ coefficient matrix ϛα1 (n)α2 (m) introduced in Eqs. (79) and (80). Once again, we
recall that ◦ stands for the Hadamard product [22]:
                                           →(τη)α i α j
    In the following, each vectorized term O jk(s)      for s = 1, 2, 3 has been reported, setting i, j = 1, 2, 3, k = 1, ..., 9 and τ, η = 0, ..., N + 1.
    First column of the generalized operator
→(τη)α i α 1 →o− 1 →(τη)[00]αi α1
O 11(1)     = A 1 ◦ A 11(20)
                                                                                               51
F. Tornabene et al.                                                                                                                Composite Structures 309 (2023) 116542
                                                                                         52
F. Tornabene et al.                                                                                                                   Composite Structures 309 (2023) 116542
                         (                                                                   )
→(τη)α i α 2 →o− 1 →o− 1    →(τη)[00]αi α2 α1 (1)α2 (0) →     →(τη)[00]αi α2 α1 (0)α2 (1) →     →o− 1 →(τη)[00]αi α2 →(τη)[01]αi α2
O 21(3)     = A 1 ◦ A 2 ◦ − A 16(11)      ◦ϛ            A 2 + A 11(20)      ◦ϛ            A 1 − R 2 ◦ A 15(11)      + A 15(10)                                      (A.30)
                                                                                        53
F. Tornabene et al.                                                          Composite Structures 309 (2023) 116542
                                                                        54
F. Tornabene et al.                                                                                                                                Composite Structures 309 (2023) 116542
                                                                                             55
F. Tornabene et al.                                                                                                                                  Composite Structures 309 (2023) 116542
        thickness by the local generalized differential quadrature method. Appl Sci 2017;           [87] Du H, Lim MK, Lin RM. Application of generalized differential quadrature to
        7:131.                                                                                           vibration analysis. J Sound Vib 1995;181:279–93.
 [70]   Alinaghizadeh F, Shariati M. Static analysis of variable thickness two-directional          [88] Tornabene F. On the critical speed evaluation of arbitrarily oriented rotating
        functionally graded annular sector plates fully or partially resting on elastic                  doubly-curved shells made of functionally graded materials. Thin-Walled Struct
        foundations by the GDQ method. J Braz Soc Mech Sci Eng 2015;37:1819–38.                          2019;140:85–98.
 [71]   Tornabene F. Free vibrations of anisotropic doubly-curved shells and panels of              [89] Tornabene F, Viscoti M, Dimitri R. Generalized higher order layerwise theory for
        revolution with a free-form meridian resting on Winkler-Pasternak elastic                        the dynamic study of anisotropic doubly-curved shells with a mapped geometry.
        foundations. Compos Struct 2011;94:186–206.                                                      Eng Anal Bound Elem 2022;134:147–83.
 [72]   Hussaini MY, Kopriva DA, Patera AT. Spectral collocation methods. Appl Numer                [90] Dimitri R, Tornabene F, Reddy JN. Numerical study of the mixed-mode behavior
        Math 1989;5:177–208.                                                                             of generally-shaped composite interfaces. Compos Struct 2020;237:111935.
 [73]   Sofiyev AH, Turan F, Kadıoglu F, Aksogan O, Hui D. Influences of two-parameter              [91] Dimitri R, Tornabene F. Numerical Study of the Mixed-Mode Delamination of
        elastic foundations on nonlinear free vibration of anisotropic shallow shell                     Composite Specimens. J Compos Sci 2018;2:30.
        structures with variable parameters. Meccanica 2022;57:401–14.                              [92] Tornabene F, Dimitri R, Viola E. Transient dynamic response of generally-shaped
 [74]   Sofiyev AH. A new approach to solution of stability problem of heterogeneous                     arches based on a GDQ-time-stepping method. Int J Mech Sci 2016;114:277–314.
        orthotropic truncated cones with clamped edges within shear deformation theory.             [93] Shu C, Chew YT, Khoo BC, Yeo KS. Solutions of three-dimensional boundary layer
        Compos Struct 2023;301:116209.                                                                   equations by global methods of generalized differential-integral quadrature. Int J
 [75]   Shu C. Differential quadrature and its application in engineering. Berlin                        Numer Meth Heat Fluid Flow 1996;227:108899.
        Heidelberg: Springer Science & Business Media; 2012.                                        [94] Adomian G. A review of the decomposition method in applied mathematics.
 [76]   Nejati M, Asanjarani A, Dimitri R, Tornabene F. Static and free vibration analysis               J Math Anal Appl 1988;135:501–44.
        of functionally graded conical shells reinforced by carbon nanotubes. Int J Mech            [95] Tornabene F, Fantuzzi N, Ubertini F, Viola E. Strong formulation finite element
        Sci 2017;130:383–98.                                                                             method based on differential quadrature: a survey. ASME Appl Mech Rev 2015;
 [77]   Du H, Lim MK, Lin R. Application of generalized differential quadrature method                   67:020801.
        to structural problems. Int J Numer Meth Eng 1994;37:1881–96.                               [96] Shu C, Du H. Implementation of clamped and simply supported boundary
 [78]   Bert CW, Malik M. The differential quadrature method for irregular domains and                   conditions in the GDQ free vibration analysis of beams and plates. Int J Solids
        application to plate vibration. Int J Mech Sci 1996;38:589–606.                                  Struct 1997;34:819–35.
 [79]   Fazzolari FA, Viscoti M, Dimitri R, Tornabene F. 1D-Hierarchical Ritz and 2D-               [97] Brischetto S, Tornabene F, Fantuzzi N, Bacciocchi M. Interpretation of boundary
        GDQ Formulations for the free vibration analysis of circular/elliptical cylindrical              conditions in the analytical and numerical shell solutions for mode analysis of
        shells and beam structures. Compos Struct 2021;258:113338.                                       multilayered structures. Int J Mech Sci 2017;122:18–28.
 [80]   Shu C, Chew YT, Richards BE. Generalized differential and integral quadrature               [98] Asadi E, Qatu MS. Static analysis of thick laminated shells with different
        and their application to solve boundary layer equations. Int J Numer Meth Fluids                 boundary conditions using GDQ. Thin-Walled Struct 2012;51:76–81.
        1995;21:723–33.                                                                             [99] Qin Z, Safaei B, Pang X, Chu F. Traveling wave analysis of rotating functionally
 [81]   Shu C, Richards BE. Application of generalized differential quadrature to solve                  graded graphene platelet reinforced nanocomposite cylindrical shells with
        two-dimensional incompressible Navier-Stokes equations. Int J Numer Meth                         general boundary conditions. Results Phys 2019;15:102752.
        Fluids 1992;15:791–8.                                                                      [100] Wang Q, Shao D, Qin B. A simple first-order shear deformation shell theory for
 [82]   Shu C, Du H. A generalized approach for implementing general boundary                            vibration analysis of composite laminated open cylindrical shells with general
        conditions in the GDQ free vibration analysis of plates. Int J Solids Struct 1997;               boundary conditions. Compos Struct 2018;184:211–32.
        34:837–46.                                                                                 [101] Singh AV, Shen L. Free vibration of open circular cylindrical composite shells with
 [83]   Shu C, Chen W, Xue H, Du H. Numerical study of grid distribution effect on                       point supports. J Aerosp Eng 2005;18:120–8.
        accuracy of DQ analysis of beams and plates by error estimation of derivative              [102] Li H, Pang F, Chen H. A semi-analytical approach to analyze vibration
        approximation. Int J Numer Meth Eng 2001;51:159–79.                                              characteristics of uniform and stepped annular-spherical shells with general
 [84]   Shu C, Chen W. On optimal selection of interior points for applying discretized                  boundary conditions. Eur J Mech-A/Solids 2019;74:48–65.
        boundary conditions in DQ vibration analysis of beams and plates. J Sound Vib              [103] Chen Y, Jin G, Liu Z. Free vibration analysis of circular cylindrical shell with non-
        1999;222:239–57.                                                                                 uniform elastic boundary constraints. Int J Mech Sci 2013;74:120–32.
 [85]   Tornabene F, Liverani A, Caligiana G. Laminated composite rectangular and                  [104] Qin Z, Zhao S, Pang X, Safaei B, Chu F. A unified solution for vibration analysis of
        annular plates: a GDQ solution for static analysis with a posteriori shear and                   laminated functionally graded shallow shells reinforced by graphene with general
        normal stress recovery. Compos B Eng 2012;43:1847–72.                                            boundary conditions. Int J Mech Sci 2020;170:105341.
 [86]   Tornabene F, Viscoti M, Dimitri R. Static analysis of anisotropic doubly-curved            [105] Tornabene F, Fantuzzi N, Bacciocchi M. DiQuMASPAB: Differential Quadrature
        shells with arbitrary geometry and variable thickness resting on a Winkler-                      for Mechanics of Anisotropic Shells. Arches and Beams, Esculapio, Bologna:
        Pasternak support and subjected to general loads. Eng Anal Bound Elem 2022;                      Plates; 2018.
        140:618–73.
56