Next Article in Journal
Magnetron-Sputtered Lead Titanate Thin Films for Pyroelectric Applications: Part 2—Electrical Characteristics and Characterization Methods
Next Article in Special Issue
The Beneficial Effect of a TPMS-Based Fillet Shape on the Mechanical Strength of Metal Cubic Lattice Structures
Previous Article in Journal
Low-Dimensional Vanadium-Based High-Voltage Cathode Materials for Promising Rechargeable Alkali-Ion Batteries
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Importance of the Recovery Procedure in the Semi-Analytical Solution for the Static Analysis of Curved Laminated Panels: Comparison with 3D Finite Elements

Department of Innovation Engineering, School of Engineering, University of Salento, 73100 Lecce, Italy
*
Author to whom correspondence should be addressed.
Materials 2024, 17(3), 588; https://doi.org/10.3390/ma17030588
Submission received: 19 December 2023 / Revised: 18 January 2024 / Accepted: 22 January 2024 / Published: 25 January 2024

Abstract

:
The manuscript presents an efficient semi-analytical solution with three-dimensional capabilities for the evaluation of the static response of laminated curved structures subjected to general external loads. A two-dimensional model is presented based on the Equivalent Single Layer (ESL) approach, where the displacement field components are described with a generalized formulation based on a higher-order expansion along the thickness direction. The fundamental equations are derived from the Hamiltonian principle, and the solution is found by means of Navier’s approach. Then, an efficient recovery procedure, derived from the three-dimensional elasticity equations and based on the Generalized Differential Quadrature (GDQ) method, is adopted for the derivation of the three-dimensional solution. Some examples of investigation are presented, where the numerical predictions of refined three-dimensional Finite-Element-based models are matched with a high level of accuracy. The model is validated for both straight and curved panels, taking into account different lamination schemes and load shapes. Furthermore, it is shown that the numerical solution to the elasticity problem in the recovery procedure is determining and accurately predicting the three-dimensional static response of the doubly-curved shell solid.

1. Introduction

Recent advancements in engineering applications require innovative strategies to accurately the static and dynamic responses of structural components of very complex shapes [1,2,3,4,5]. For this reason, advanced models are necessary to describe the geometry and the mechanical characteristics of these structures with reduced computational costs [6,7]. Three-dimensional approaches based on the elasticity equations provide highly accurate predictions of the structural response of a doubly-curved solid, but they can be computationally expensive [8,9]. On the other hand, two-dimensional formulations that consider a higher-order description of the field variable can be a valid alternative to three-dimensional models [10,11,12,13]. In a two-dimensional theory, a doubly-curved surface is considered to have equivalent mechanical properties instead of a doubly-curved three-dimensional solid [14,15,16]. The unknown field variables are described using either the Equivalent Single Layer (ESL) or the Layer-Wise (LW) approaches [17,18,19,20,21]. More specifically, in ESL theories, a higher-order expansion is established along the thickness direction by means of the so-called thickness functions. However, for moderately thick and thick panels, a LW description of the field variable may be more accurate. According to the LW methodology, the fundamental relations are written in each layer of shell, and the interaction between adjacent laminae is taken into account through compatibility conditions. Referring to two-dimensional ESL theories, when the shell laminated structure is made of advanced materials, it is very likely that a higher-order expansion of the unknown field variable is required [22,23,24], leading to the well-known Higher-Order Shear Deformation Theories (HSDTs). The higher-order expansion in hand can be performed with both polynomial and non-polynomial thickness functions [25,26,27,28]. Besides, the coupling between two-adjacent layers, which is known as zigzag effects, can be easily modeled in ESL theories with the so-called zigzag thickness functions, as shown in Refs. [29,30,31]. In this way, a piecewise inclination of the displacement field profile is provided. The approach provided for the first time in Ref. [32] has been shown to be very simple and accurate, among others. On the other hand, in Refs. [33,34,35], the refined zigzag theories are presented, where the zigzag functions are derived from the stiffness properties of the lamination scheme of the selected panel, leading to the so-called refined zigzag theories.
Among two-dimensional approaches, an efficient strategy to define the fundamental relations can be found in the well-known generalized formulation [36,37,38,39], where the through-the-thickness expansion of the configuration variable up to an arbitrary order is modeled regardless of the effective expression of the selected thickness function. In this way, several structural models with different kinematic assumptions can be derived. Furthermore, classical formulations like the First Order Shear Deformation Theory (FSDT) and the Third Order Shear Deformation Theory (TSDT) can be seen as particular cases of the generalized higher-order theory [40,41]. The adoption of a generalized configuration variable can be an efficient strategy when structures of different materials and various lamination schemes are studied with an arbitrary variation of mechanical properties, as happens in the case of Functionally Graded Materials (FGMs) [42,43,44], Carbon Nanotube (CNT) composites [45,46,47], honeycomb and lattice cells [48,49], as well as three-dimensional Variable Angle Tow (3D-VAT) composites [50,51,52]. Furthermore, the adoption of HSDTs allows one to study the effect of porosity because the presence of voids with an arbitrary distribution within the structure leads to a reduction in material stiffness [53,54].
For structures with complex shapes, materials, and boundary conditions, it is very difficult to derive a closed-form solution to the governing equations; where numerical approximations must be used to discretize the problem [55,56,57,58]. On the other hand, it is possible to derive a closed-form solution if some simplifications are taken into account [59,60,61,62,63]. In addition, some applications of practical interest can be examined with an infinite series expansion of each unknown variable [64,65,66,67]. It should be noted that exact solutions for simply supported layered structures are typically adopted to check the accuracy of the numerical solution. Among semi-analytical solutions for linear elasticity, a milestone is the research work by Pagano [68,69], where a three-dimensional closed-form solution is derived for laminated composite plates and sandwich panels. Then, two-dimensional models have been developed for laminated plates and curved sandwich shells. Starting from formulations based on the Classical Plate Theory (CPT), refined theories can be found in literature based on FSDT and TSDT. Recent developments in the field of smart materials have led to the development of new formulations regarding plates and shells with smart properties like piezoelectricity, magnetostriction, and heat transfer [70], and closed-form solutions have been derived for the validation of numerical simulations. Some preliminary works regarding laminated plates with piezoelectric [71,72] and piezomagnetic [73,74] properties must be cited.
It is interesting to note that closed-form solutions may not be used for practical applications because of the assumptions that are usually considered. However, some results of more practical interest can be obtained with semi-analytical formulations. In a semi-analytical theory, the solution is obtained with an expansion of degrees of freedom (DOFs) up to a sufficient order. The Navier method and the Levy procedure have been extensively adopted in several papers regarding linear elasticity problems for plates and shells [75,76,77,78]. More specifically, the Navier solution, based on the description of the field variable with a trigonometric series, can be adopted in the case of simply-supported laminated panels with cross-ply lamination schemes and antisymmetric angle-ply composites. In contrast, the Levy method [79,80] is suitable for panels with two simply supported edges and the other two subjected to arbitrary boundary conditions. An arbitrary load case can be analyzed with the Navier method if the external actions are expanded with a double Fourier trigonometric series. It has been shown in Ref. [81] that this methodology can also be adopted for closed panels of revolution. On the other hand, the Navier approach can be difficult to apply to structures of very complex shapes made of non-conventional materials because a significantly high number of terms of the harmonic expansion should be considered to obtain a sufficient level of accuracy. For this reason, a numerical model is more frequently adopted to derive an approximate solution under a limited number of hypotheses. Among numerical techniques, the Finite Element Method (FEM), which is extensively adopted in several applications and commercial codes, is based on a local interpolation of the unknown field variable in a discrete computational grid by means of the so-called shape functions [82,83,84]. In contrast, the class of spectral collocation methods [85,86,87] are based on a global interpolation of the solution by means of higher-order functions; therefore, a smoother solution can be derived with a reduced number of DOFs. Belonging to this class, the Generalized Differential Quadrature (GDQ) method [88,89,90,91,92] approximates the derivatives of an arbitrary function as a weighted sum of the function itself. It has been shown in several papers that the highest level of accuracy, computational stability, and efficiency is reached when the computational domain is discretized with non-uniform grids [93,94]. Moving from the GDQ numerical technique, the Generalized Integral Quadrature (GIQ) allows one to compute the integrals with a significantly reduced number of DOFs [95,96].
When a two-dimensional solution is derived, the reconstruction of the effective structural response of the panel can be difficult since it is not guaranteed, a priori, that the distribution of stress components satisfies the three-dimensional elasticity equations. This issue may lead to erroneous results, especially for moderately thick structures, where the stress components acting in the thickness direction cannot be neglected. For this reason, a correction of the stress and strain profiles should be conducted based on the three-dimensional elasticity equations. In Refs. [97,98,99], an effective stress and strain recovery procedure is provided for the evaluation of the static response of moderately thick doubly-curved shell structures made of laminated composite materials with arbitrary orientation and FGM, starting from a refined two-dimensional GDQ-based numerical solution. The recovery procedure has been demonstrated to be an effective tool for the reconstruction of the three-dimensional response of doubly-curved shell structures with advanced lamination schemes, starting from numerical solutions of refined formulations based on HSDTs. In some previous works [100,101], the recovery procedure has been applied to some two-dimensional GDQ-based numerical solutions for the prediction of the static response of doubly-curved shells with generally anisotropic materials. However, the effects of this procedure on two-dimensional semi-analytical solutions have not been checked.
In the present work, a two-dimensional model based on the ESL approach with HSDTs is presented to predict the linear static response of laminated curved panels. The geometry of the structure is described with the differential geometry basics and curvilinear principal coordinates. The generalized formulation is adopted for the description of the kinematic field, and a higher-order expression is provided along the thickness direction of the panel for each displacement field component. Furthermore, the zigzag effects are considered in the kinematic model. Following the ESL approach, the mechanical properties of each layer, modeled with a generally anisotropic constitutive relationship, are homogenized on the reference surface. The fundamental equations are derived from the Hamiltonian principle, accounting for an arbitrary distribution of the external surface loads, which are applied at the top and bottom of the laminated panel. Then, the semi-analytical Navier solution is found under some geometric and mechanical assumptions, taking into account a Fourier series expansion of the unknown field variable. Finally, a recovery procedure based on the three-dimensional elasticity equations for a doubly-curved solid is applied for the reconstruction of the three-dimensional response of the panel. Some examples of investigation are presented, where the accuracy of the semi-analytical formulation is checked for different curvatures, lamination schemes, and load cases. The numerical predictions have been performed with various kinematic field assumptions, and the results are compared to those coming from three-dimensional finite element models. Furthermore, some comparisons are conducted with the results of refined GDQ numerical models. It is shown that the recovery procedure determines when a semi-analytical procedure is adopted. In this way, it is possible to accurately predict the three-dimensional response of curved laminated structures with a reduced computational cost even though a two-dimensional formulation is adopted. On the other hand, if the stress and strain profiles are derived directly from the two-dimensional solution without any post-processing technique, some results may be obtained that are not consistent with the equilibrium of the structure under external loads.

2. Geometry of a Shell Structure

A doubly-curved shell is a three-dimensional solid within the Euclidean space whose position vector, denoted by R α 1 , α 2 , α 3 , is dependent on three parameters α i = α 1 , α 2 , α 3 . These parameters are defined in the closed interval α i α i 0 , α i 1 with i = 1 , 2 , 3 , where α i 0 < α i 1 are the extremes of the domain at issue. Following the ESL approach, the position vector R of a doubly-curved shell element is defined in terms of a reference surface, denoted by r α 1 , α 2 , located in the middle thickness h α 1 , α 2 of the panel:
R α 1 , α 2 , ζ = r α 1 , α 2 + h α 1 , α 2 2 z n α 1 , α 2
In the previous relation, z = 2 ζ / h is a dimensionless coordinate for the thickness direction, oriented alongside the outward normal unit vector n α 1 , α 2 . This vector can be calculated at each point of the reference surface in terms of the partial derivatives of r α 1 , α 2 with respect to the principal coordinates α 1 , α 2 , denoted by the symbol r , i = r / α i with i = 1 , 2 . One gets [16]:
n α 1 , α 2 = r , 1 × r , 2 r , 1 × r , 2
Following Equation (1), the geometric properties of the three-dimensional shell element can thus be defined from those of the reference surface r α 1 , α 2 . The principal radii of curvature R i α 1 , α 2 = R 1 , R 2 , referred to the α i = α 1 , α 2 principal direction, can thus be derived as follows:
R i α 1 , α 2 = r , i r , i r , i i n i = 1 , 2
Note that the symbol r , i j = 2 r / α i α j with i , j = 1 , 2 refers to the second order derivatives of the reference surface r with respect to α i , α j = α 1 , α 2 . For the sake of completeness, the principal curvature k n i = 1 / R i with i = 1 , 2 can be introduced along each principal direction. The Lamè parameters A i α 1 , α 2 = A 1 , A 2 can be computed at each point of the physical domain according to the following relation:
A i α 1 , α 2 = r , i r , i i = 1 , 2
The Lamè parameters A i = A 1 , A 2 are used for the computation of the infinitesimal arch lengths d s i = d s 1 , d s 2 along the principal directions α i = α 1 , α 2 of the reference surface. The arch length s i with i = 1 , 2 is defined so that s i 0 , L i , being L i the length of the parametric line. The following relation can thus be written:
d s i = A i d α i i = 1 , 2
On the other hand, the scaling parameters H i = H 1 , H 2 are defined in order to evaluate the scaling effects along the thickness direction that can be seen in a solid with one or more curvatures:
H i α 1 , α 2 , ζ = 1 + ζ R i i = 1 , 2
The three-dimensional shell solid is obtained from the superimposition of l different layers of thickness h k , therefore the total thickness of the structure is thus obtained as the sum of the width of each k -th lamina, with k = 1 , , l , which is located between its intrados and its extrados, denoted by ζ k and ζ k + 1 , respectively [16]:
h α 1 , α 2 = k = 1 l h k α 1 , α 2 = k = 1 l ζ k + 1 ζ k

3. Kinematic Relations

In the present section, a generalized ESL approach is presented for the expansion up to the N + 1 -th order of the three-dimensional displacement field vector U k α 1 , α 2 , ζ = U 1 k U 2 k U 3 k T along the thickness direction. To this end, three thickness functions F k τ α i ζ = F k τ α 1 , F k τ α 2 , F k τ α 3 are introduced for each τ -th kinematic expansion order with τ = 0 , , N + 1 . Remembering that the thickness coordinate is defined so that ζ ζ k , ζ k + 1 for k = 1 , , l , a generalized formulation can be derived, and the vector u τ α 1 , α 2 of the generalized displacement field components u 1 τ , u 2 τ , u 3 τ is introduced for each τ = 0 , , N + 1 . One gets the following relation [16]:
U k α 1 , α 2 , ζ = τ = 0 N + 1 F k τ ζ u τ α 1 , α 2     U 1 k U 2 k U 3 k = τ = 0 N + 1 F k τ α 1 0 0 0 F k τ α 2 0 0 0 F k τ α 3 u 1 τ u 2 τ u 3 τ
The ESL formulation of the previous relation allows one to derive a generalized two-dimensional theory for the mechanical elasticity problem. The selection of a particular expression of the thickness functions allows one to obtain several models for the static analysis of shell structures, including classical approaches like the CPT, FSDT, and TSDT. On the other hand, the thickness function related to τ = N + 1 simulates the zig-zag effect that occurs in the interlaminar region, which consists of an abrupt variation of the profile of each displacement field component. In the present work, power thickness functions are used for each τ = 0 , , N , together with Murakami’s zigzag function associated to τ = N + 1 :
F k τ α i ζ = ζ τ τ = 0 , , N 1 k z k = 1 k 2 ζ k + 1 ζ k ζ ζ k + 1 + ζ k ζ k + 1 ζ k τ = N + 1
When the thickness functions of Equation (9) are used, the notation EDZ- N can be used for the identification of the higher-order theory [16]. More specifically, “E” means that the two-dimensional theory is based on the ESL approach, whereas “D” means that an axiomatic expansion is adopted of the displacement field components, which are the configuration variables of the problem. When the zig-zag function is considered for τ = N + 1 , the letter “Z” is used. Finally, N denotes the maximum expansion order of the configuration variable.
At this point, the kinematic relations are derived for the ESL theory from those of the three-dimensional elasticity problem, as shown in Ref. [16]. If ε k α 1 , α 2 , ζ = ε 1 k ε 2 k γ 12 k γ 13 k γ 23 k ε 3 k T is the three-dimensional strain vector of the k -th layer, the following condensed relation can be taken into account:
ε k α 1 , α 2 , ζ = D U k = D ζ i = 1 3 D Ω α i U k
In the previous relation, D is the kinematic differential operator, which is split in those denoted by D ζ and D Ω α i = D Ω α 1 , D Ω α 2 , D Ω α 3 , which contain the derivatives with respect to the coordinate ζ and α 1 , α 2 , respectively. In an extended form, D ζ can be expressed as:
D ζ = 1 H 1 0 0 0 0 0 0 0 0 0 1 H 2 0 0 0 0 0 0 0 0 0 1 H 1 1 H 2 0 0 0 0 0 0 0 0 0 1 H 1 0 ζ 0 0 0 0 0 0 0 1 H 2 0 ζ 0 0 0 0 0 0 0 0 0 ζ
On the other hand, the kinematic operators D Ω α 1 , D Ω α 2 , D Ω α 3 of size 9 × 3 take the following aspect:
D Ω α 1 = D ¯ Ω α 1 0 0 D Ω α 2 = 0 D ¯ Ω α 2 0 D Ω α 3 = 0 0 D ¯ Ω α 3
The sub-operators D ¯ Ω α 1 , D ¯ Ω α 2 , D ¯ Ω α 3 are written in an extended form as:
D ¯ Ω α 1 = 1 A 1 α 1 1 A 1 A 2 A 2 α 1 1 A 1 A 2 A 1 α 2 1 A 2 α 2 1 R 1 0 1 0 0 , D ¯ Ω α 2 = 1 A 1 A 2 A 1 α 2 1 A 2 α 2 1 A 1 α 1 1 A 1 A 2 A 2 α 1 0 1 R 2 0 1 0 , D ¯ Ω α 3 = 1 R 1 1 R 2 0 0 1 A 1 α 1 1 A 2 α 2 0 0 1
Introducing Equation (8) in Equation (10), the generalized ESL kinematic relations are obtained, accounting for the effects of the curvature of the shell and the higher-order kinematic expansion of the displacement field variable [16]:
ε k = τ = 0 N + 1 i = 1 3 Z k τ α i ε τ α i
As can be seen, the three-dimensional strain vector ε k α 1 , α 2 , ζ has been expressed in terms of the generalized strain vector, for each τ = 0 , , N + 1 , denoted by ε τ α i α 1 , α 2 = ε 1 τ α i ε 2 τ α i γ 1 τ α i γ 2 τ α i γ 13 τ α i γ 23 τ α i ω 13 τ α i ω 23 τ α i ε 3 τ α i T , whose definition is reported in the following:
ε τ α i = D Ω α i u τ   for   τ = 0 , , N + 1 , i = 1 , 2 , 3
On the other hand, the ESL kinematic operator Z k τ α i is defined for each τ = 0 , , N + 1 in terms of the generalized thickness functions of Equation (8):
Z k τ α i = D ζ F k τ α i = F k τ α i H 1 0 0 0 0 0 0 0 0 0 F k τ α i H 2 0 0 0 0 0 0 0 0 0 F k τ α i H 1 F k τ α i H 2 0 0 0 0 0 0 0 0 0 F k τ α i H 1 0 F k τ α i ζ 0 0 0 0 0 0 0 F k τ α i H 2 0 F k τ α i ζ 0 0 0 0 0 0 0 0 0 F k τ α i ζ

4. Constitutive Relationship

In the present section, the constitutive behavior of the three-dimensional shell solid is described within the two-dimensional ESL model employing higher-order theories. Referring to an arbitrary k -th layer of the solid, a three-dimensional linear elastic constitutive relation is considered, defined as follows [16]:
σ k = E ¯ k ε k σ 1 k σ 2 k τ 12 k τ 13 k τ 23 k σ 3 k = E ¯ 11 k E ¯ 12 k E ¯ 16 k E ¯ 14 k E ¯ 15 k E ¯ 13 k E ¯ 12 k E ¯ 22 k E ¯ 26 k E ¯ 24 k E ¯ 25 k E ¯ 23 k E ¯ 16 k E ¯ 26 k E ¯ 66 k E ¯ 46 k E ¯ 56 k E ¯ 36 k E ¯ 14 k E ¯ 24 k E ¯ 46 k E ¯ 44 k E ¯ 45 k E ¯ 34 k E ¯ 15 k E ¯ 25 k E ¯ 56 k E ¯ 45 k E ¯ 55 k E ¯ 35 k E ¯ 13 k E ¯ 23 k E ¯ 36 k E ¯ 34 k E ¯ 35 k E ¯ 33 k ε 1 k ε 2 k γ 12 k γ 13 k γ 23 k ε 3 k           k = 1 , , l
In the previous relation, E ¯ k denotes the three-dimensional stiffness matrix, whose generic component E ¯ i j k with i , j = 1 , , 6 relates the i -th element of the stress vector, denoted by σ k α 1 , α 2 , ζ , to the j -th element of the strain vector ε k α 1 , α 2 , ζ . The constitutive behavior of each lamina is usually provided not in the geometric reference system, as happens in Equation (17), but in the material reference system O α ^ 1 α ^ 2 α ^ 3 , whose axes are oriented along the material symmetry directions. For this reason, the equation reported in the following should be considered, where σ ^ k and ε ^ k are the vectors of the three-dimensional stress and strain components, respectively, written with respect to O α ^ 1 α ^ 2 α ^ 3 material reference system, whereas E k is the corresponding stiffness matrix:
σ ^ k = E k ε ^ k σ ^ 1 k σ ^ 2 k τ ^ 12 k τ ^ 13 k τ ^ 23 k σ ^ 3 k = E 11 k E 12 k E 16 k E 14 k E 15 k E 13 k E 12 k E 22 k E 26 k E 24 k E 25 k E 23 k E 16 k E 26 k E 66 k E 46 k E 56 k E 36 k E 14 k E 24 k E 46 k E 44 k E 45 k E 34 k E 15 k E 25 k E 56 k E 45 k E 55 k E 35 k E 13 k E 23 k E 36 k E 34 k E 35 k E 33 k ε ^ 1 k ε ^ 2 k γ ^ 12 k γ ^ 13 k γ ^ 23 k ε ^ 3 k           k = 1 , , l
In the previous equation, E i j k with i , j = 1 , , 6 are the elements of the matrix E k . More specifically, they can be taken as equal to the three-dimensional stiffness components of the material of the -th layer, namely E i j k = C i j k . However, when the kinematic field assumption of Equation (8) neglects the stretching effect, E i j k turn out to be the reduced elastic coefficients Q i j k , which are calculated from the three-dimensional ones according to what is exerted in Ref. [16].
Finally, Equation (18) can revert to Equation (17) if the stiffness matrix E k is rotated by means of the transformation matrix T k , as defined in Ref. [16]. One gets:
E ¯ k = T k E k T k T
In the present study, the material axis α ^ 3 is taken along the thickness direction of the shell; therefore, the equivalence α ^ 3 = ζ is considered. As a consequence, the rotation matrix T k depends only on the angle ϑ k between α 1 and α ^ 1 reference axes such that T k = T k ϑ k . Introducing the generalized higher-order expansion of the kinematic relations of Equation (14) in the three-dimensional constitutive relationship of Equation (17), the following higher-order constitutive relationship comes out:
σ k = E ¯ k ε k = η = 0 N + 1 j = 1 3 E ¯ k Z k η α j ε η α j       for       k = 1 , , l
At this point, the previous relation can be used for the evaluation of the virtual variation δ Φ of the elastic strain energy of the shell:
δ Φ = k = 1 l α 1 α 2 ζ k ζ k + 1 δ ε k T σ k A 1 A 2 H 1 H 2 d α 1 d α 2 d ζ
The computation of the integrals of Equation (21) in the closed interval ζ k , ζ k + 1 allows one to introduce the generalized stress resultants, which are collected for each τ = 0 , , N + 1 and i = 1 , 2 , 3 in the generalized stress resultant vector S τ α i = 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 T [16]:
δ Φ = τ = 0 N + 1 i = 1 3 α 1 α 2 δ ε τ α i T S τ α i A 1 A 2 d α 1 d α 2
Finally, the higher-order ESL elastic constitutive relationship can be written in terms of S τ α i and ε τ α i by substituting Equation (20) in the final expression of the virtual variation of the elastic strain energy δ Φ , as reported in Equation (22). One gets:
S τ α i = η = 0 N + 1 j = 1 3 A τ η α i α j ε η α j for τ = 0 , , N + 1 , i = 1 , 2 , 3
In the previous equation, A τ η α i α j denotes the generalized higher-order constitutive operator, which is defined for each τ , η = 0 , , N + 1 and i , j = 1 , 2 , 3 as follows [16]:
A τ η α i α j = k = 1 l ζ k ζ k + 1 Z k τ α i T E ¯ k Z k η α j H 1 H 2 d ζ = A τ η 00 α i α j A τ η 01 α i α j A τ η 10 α i α j A τ η 11 α i α j
For the sake of clarity, the sub-matrices A τ η 00 α i α j , A τ η 01 α i α j , A τ η 10 α i α j , A τ η 11 α i α j are reported below in an extended form:
A τ η 00 α i α j = A 11 20 11 τ η 00 α i α j A 12 11 12 τ η 00 α i α j A 16 20 13 τ η 00 α i α j A 16 11 14 τ η 00 α i α j A 14 20 τ η 00 α i α j A 15 11 τ η 00 α i α j A 12 11 τ η 00 α i α j A 22 02 τ η 00 α i α j A 26 11 τ η 00 α i α j A 26 02 τ η 00 α i α j A 24 11 τ η 00 α i α j A 25 02 τ η 00 α i α j A 16 20 τ η 00 α i α j A 26 11 τ η 00 α i α j A 66 20 τ η 00 α i α j A 66 11 τ η 00 α i α j A 46 20 τ η 00 α i α j A 56 11 τ η 00 α i α j A 16 11 τ η 00 α i α j A 26 02 τ η 00 α i α j A 66 11 τ η 00 α i α j A 66 02 τ η 00 α i α j A 46 11 τ η 00 α i α j A 56 02 τ η 00 α i α j A 14 20 τ η 00 α i α j A 24 11 τ η 00 α i α j A 46 20 τ η 00 α i α j A 46 11 τ η 00 α i α j A 44 20 τ η 00 α i α j A 45 11 τ η 00 α i α j A 15 11 τ η 00 α i α j A 25 02 τ η 00 α i α j A 56 11 τ η 00 α i α j A 56 02 τ η 00 α i α j A 45 11 τ η 00 α i α j A 55 02 τ η 00 α i α j
A τ η 01 α i α j = A 14 10 τ η 01 α i α j A 15 10 τ η 01 α i α j A 13 10 τ η 01 α i α j A 24 01 τ η 01 α i α j A 25 01 τ η 01 α i α j A 23 01 τ η 01 α i α j A 46 10 τ η 01 α i α j A 56 10 τ η 01 α i α j A 36 10 τ η 01 α i α j A 46 01 τ η 01 α i α j A 56 01 τ η 01 α i α j A 36 01 τ η 01 α i α j A 44 10 τ η 01 α i α j A 45 10 τ η 01 α i α j A 34 10 τ η 01 α i α j A 45 01 τ η 01 α i α j A 55 01 τ η 01 α i α j A 35 01 τ η 01 α i α j
A τ η 10 α i α j = A 14 10 τ η 10 α i α j A 24 01 τ η 10 α i α j A 46 10 τ η 10 α i α j A 46 01 τ η 10 α i α j A 44 10 τ η 10 α i α j A 45 01 τ η 10 α i α j A 15 10 τ η 10 α i α j A 25 01 τ η 10 α i α j A 56 10 τ η 10 α i α j A 56 01 τ η 10 α i α j A 45 10 τ η 10 α i α j A 55 01 τ η 10 α i α j A 13 10 τ η 10 α i α j A 23 01 τ η 10 α i α j A 36 10 τ η 10 α i α j A 36 01 τ η 10 α i α j A 34 10 τ η 10 α i α j A 35 01 τ η 10 α i α j
A τ η 11 α i α j = A 44 00 τ η 11 α i α j A 45 00 τ η 11 α i α j A 34 00 τ η 11 α i α j A 45 00 τ η 11 α i α j A 55 00 τ η 11 α i α j A 35 00 τ η 11 α i α j A 34 00 τ η 11 α i α j A 35 00 τ η 11 α i α j A 33 00 τ η 11 α i α j
The generalized elastic coefficients of Equations (25)–(28) can be computed with the following condensed expression, setting the definitions 0 F k τ α i / ζ 0 = F k τ α i and 0 F k η α j / ζ 0 = F k η α j :
A n m p q τ η f g α i α j = k = 1 l ζ k ζ k + 1 B ¯ n m k f F k η α j ζ f g F k τ α i ζ g H 1 H 2 H 1 p H 2 q d ζ for τ , η = 0 , , N + 1 , n , m = 1 , , 6 , p , q = 0 , 1 , 2 , α i , α j = α 1 , α 2 , α 3 , f , g = 0 , 1
The quantities B ¯ n m k with n , m = 1 , , 6 occurring in the previous expression denote the three-dimensional elastic stiffness coefficients E ¯ n m k of Equation (17). As can be seen, their value is corrected with the shear correction factor, denoted by κ ζ , in order to consider the effects of shear stresses to the global deflection of the structure when lower order displacement field assumptions are considered:
B ¯ n m ( k ) = E ¯ n m ( k ) κ ζ E ¯ n m ( k ) for n , m = 1 , 2 , 3 , 6 , n , m = 4 , 5
Accordingly, when a linear distribution of the displacement field components is considered in Equation (8), the value κ ζ = 5 / 6 is assumed, whereas in the case of higher-order theories, a unitary value κ ζ = 1 is assigned to this quantity. The value κ ζ = 5 / 6 is here assumed, as in classical formulations for beams with rectangular cross-sections. The higher-order constitutive relation of Equation (23) can be written in terms of the kinematic field assumption of Equation (8) in order to provide a relationship between the stress resultant vector S τ α i and the generalized displacement field vector u τ . Substituting the definition of the vector ε τ α i of the generalized strain components of Equation (15) in the higher-order constitutive relation of Equation (23), the relation reported in the following is obtained:
S τ α i = η = 0 N + 1 j = 1 3 A τ η α i α j D Ω α j u η = η = 0 N + 1 j = 1 3 O τ η α i α j u η for τ = 0 , , N + 1 , i = 1 , 2 , 3
The matrix O τ η α i α j of size 9 × 3 , defined for each τ , η = 0 , , N + 1 and α i , α j = α 1 , α 2 , α 3 is reported in the following in an extended form:
O τ η α i α j = O 11 τ η α i α 1 O 12 τ η α i α 2 O 13 τ η α i α 3 O 21 τ η α i α 1 O 22 τ η α i α 2 O 23 τ η α i α 3 O 31 τ η α i α 1 O 32 τ η α i α 2 O 33 τ η α i α 3 O 41 τ η α i α 1 O 42 τ η α i α 2 O 43 τ η α i α 3 O 51 τ η α i α 1 O 52 τ η α i α 2 O 53 τ η α i α 3 O 61 τ η α i α 1 O 62 τ η α i α 2 O 63 τ η α i α 3 O 71 τ η α i α 1 O 72 τ η α i α 2 O 73 τ η α i α 3 O 81 τ η α i α 1 O 82 τ η α i α 2 O 83 τ η α i α 3 O 91 τ η α i α 1 O 92 τ η α i α 2 O 93 τ η α i α 3 for τ , η = 0 , , N + 1 , i , j = 1 , 2 , 3
The complete expression of the coefficients O g r τ η α i α j with g = 1 , , 9 and r = 1 , 2 , 3 can be found in Ref. [16].

5. Governing Equations

Once the kinematic and constitutive relations have been presented, the fundamental governing equations are derived for the linear static analysis of doubly-curved shell structures. Following an energetic procedure, the equilibrium configuration of the solid is derived from the minimum potential energy principle, taking into account the elastic deformation energy of the system, denoted by Φ , and the virtual work L e of external loads. If the virtual variation of each physical quantity is denoted by δ , the following relation is considered [16]:
δ Φ δ L e = 0
As shown in Equation (22), the virtual variation δ Φ of the elastic strain energy is written in terms of the virtual variation of the vector δ ε τ α i of the generalized strain components. Applying the integration by parts rule, an expression of the variation δ Φ is obtained in terms of the virtual variation of the generalized displacement field components δ u 1 τ , δ u 2 τ , δ u 3 τ :
δ Φ = τ = 0 N + 1 α 1 α 2 N 1 τ α 1 A 2 α 1 + N 21 τ α 1 A 1 α 2 + N 12 τ α 1 A 1 α 2 N 2 τ α 1 A 2 α 1 + T 1 τ α 1 R 1 P 1 τ α 1 A 1 A 2 δ u 1 τ + N 2 τ α 2 A 1 α 2 + N 12 τ α 2 A 2 α 1 + N 21 τ α 2 A 2 α 1 N 1 τ α 2 A 1 α 2 + T 2 τ α 2 R 2 P 2 τ α 2 A 1 A 2 δ u 2 τ + T 1 τ α 3 A 2 α 1 + T 2 τ α 3 A 1 α 2 N 1 τ α 3 R 1 + N 2 τ α 3 R 2 A 1 A 2 S 3 τ α 3 A 1 A 2 δ u 3 τ d α 1 d α 2 + + τ = 0 N + 1 α 1 N 21 τ α 1 δ u 1 τ + N 2 τ α 2 δ u 2 τ + T 2 τ α 3 δ u 3 τ A 1 d α 1 + + τ = 0 N + 1 α 2 N 1 τ α 1 δ u 1 τ + N 12 τ α 2 δ u 2 τ + T 1 τ α 3 δ u 3 τ A 2 d α 2
The virtual work of external actions, denoted by δ L e 3 D , is computed as the sum of the virtual work of the actions q i a + = q i a 1 and q i a = q i a 2 with i = 1 , 2 , 3 acting at the top j = 1 and the bottom j = 2 of the shell, respectively. Referring to a doubly-curved three-dimensional solid, one gets:
δ L e 3 D = α 1 α 2 j = 1 2 i = 1 3 q i a j δ U i j H 1 j H 2 j A 1 A 2 d α 1 d α 2
In the previous relation, U i + = U i 1 and U i = U i 2 with i = 1 , 2 , 3 are the virtual variations of the displacement field components at the top ζ = h / 2 and the bottom ζ = h / 2 surfaces of the three-dimensional solid, respectively. The application of the static equivalence principle introduces a set of generalized loads for each τ = 0 , , N + 1 kinematic expansion order, which are collected in the vector q τ α 1 , α 2 = q 1 τ q 2 τ q 3 τ T . They are associated to the virtual variation of the vector u τ α 1 , α 2 of the generalized displacement field components. The following relation is thus obtained:
δ L e = τ = 0 N + 1 α 1 α 2 δ u τ T q τ A 1 A 2 d α 1 d α 2 = τ = 0 N + 1 α 1 α 2 i = 1 3 q i τ δ u i τ A 1 A 2 d α 1 d α 2
According to the static equivalence principle, the virtual work δ L e 3 D of Equation (35) turns out to be equal to the virtual work δ L e of Equation (36):
δ L e 3 D = δ L e
Substituting in Equation (35) the kinematic assumptions of Equation (8) and introducing them in Equation (37), the following expression is derived for the generalized loads q i τ = q 1 τ , q 2 τ , q 3 τ [16]:
q i τ = j = 1 2 q i a j F τ α i j H 1 j H 2 j for     i = 1 , 2 , 3
In the previous equation, the quantity F τ α i j with j = 1 , 2 denotes the thickness function associated to U i + and U i , whereas H 1 j , H 2 j with j = 1 , 2 are the scaling parameters calculated at the top ζ = h / 2 and the bottom ζ = h / 2 surfaces of the shell.
Starting from Equation (35), it is possible to embed in the present formulation a general surface load. To this end, an arbitrary surface distribution q ¯ α 1 , α 2 is considered, whereas the magnitude of the external load at issue is denoted by q i j
q i a j α 1 , α 2 , ± h 2 = q i j ζ = ± h / 2 q ¯ α 1 , α 2 for i = 1 , 2 , 3 j = 1 , 2
In the case of a constant distribution of the external load, the distribution q ¯ α 1 , α 2 = 1 is considered.
Introducing in Equation (33) the expression of the virtual work of external actions and remembering Equation (34), the higher-order equilibrium equations are derived on the shell reference surface. Employing a compact notation, it gives:
i = 1 3 D Ω * α i S τ α i q τ = 0 for   τ = 0 , , N + 1
In the previous relation, D Ω * α i = D Ω * α 1 , D Ω * α 2 , D Ω * α 3 denote the equilibrium operators, which can be expressed with a matrix notation as follows:
D Ω * α 1 = D ¯ Ω * α 1 0 0 , D Ω * α 2 = 0 D ¯ Ω * α 2 0 , D Ω * α 3 = 0 0 D ¯ Ω * α 3
An extended version of the quantities D ¯ Ω * α i = D ¯ Ω * α 1 , D ¯ Ω * α 2 , D ¯ Ω * α 3 is reported below:
D ¯ Ω * α 1 = D ¯ Ω * α 1 1 D ¯ Ω * α 1 2 D ¯ Ω * α 1 3 D ¯ Ω * α 1 4 D ¯ Ω * α 1 5 D ¯ Ω * α 1 6 D ¯ Ω * α 1 7 D ¯ Ω * α 1 8 D ¯ Ω * α 1 9 D ¯ Ω * α 2 = D ¯ Ω * α 2 1 D ¯ Ω * α 2 2 D ¯ Ω * α 2 3 D ¯ Ω * α 2 4 D ¯ Ω * α 2 5 D ¯ Ω * α 2 6 D ¯ Ω * α 2 7 D ¯ Ω * α 2 8 D ¯ Ω * α 2 9 D ¯ Ω * α 3 = D ¯ Ω * α 3 1 D ¯ Ω * α 3 2 D ¯ Ω * α 3 3 D ¯ Ω * α 3 4 D ¯ Ω * α 3 5 D ¯ Ω * α 3 6 D ¯ Ω * α 3 7 D ¯ Ω * α 3 8 D ¯ Ω * α 3 9
where each term D ¯ Ω * α i with i = 1 , 2 , 3 reads as:
D ¯ Ω * α 1 1 = D ¯ Ω * α 2 3 = D ¯ Ω * α 3 5 = 1 A 1 α 1 + 1 A 1 A 2 A 2 α 1 , D ¯ Ω * α 1 4 = D ¯ Ω * α 2 2 = D ¯ Ω * α 3 6 = 1 A 2 α 2 + 1 A 1 A 2 A 1 α 2 , D ¯ Ω * α 1 3 = D ¯ Ω * α 2 1 = 1 A 1 A 2 A 1 α 2 , D ¯ Ω * α 1 2 = D ¯ Ω * α 2 4 = 1 A 1 A 2 A 2 α 1 , D ¯ Ω * α 1 5 = D ¯ Ω * α 3 1 = 1 R 1 , D ¯ Ω * α 2 6 = D ¯ Ω * α 3 2 = 1 R 2 , D ¯ Ω * α 1 7 = D ¯ Ω * α 2 8 = D ¯ Ω * α 3 9 = 1 , D ¯ Ω * α 1 6 = D ¯ Ω * α 1 8 = D ¯ Ω * α 1 9 = D ¯ Ω * α 2 5 = D ¯ Ω * α 2 7 = D ¯ Ω * α 2 9 = D ¯ Ω * α 3 3 = D ¯ Ω * α 3 4 = D ¯ Ω * α 3 7 = D ¯ Ω * α 3 8 = 0
Introducing in Equation (40) the definition of S τ α i in terms of u τ , as happens in Equation (31), the fundamental equations are derived in each point of the physical domain for the static analysis of doubly-curved shells with higher-order theories for each τ = 0 , , N + 1 [16]:
η = 0 N + 1 L τ η u η = q τ for τ = 0 , , N + 1
The fundamental matrix L τ η referred to an arbitrary τ , η = 0 , , N + 1 occurring in Equation (44) turns out to be of size 3 × 3 , reading as follows [16]:
L τ η = L 11 τ η α 1 α 1 L 12 τ η α 1 α 2 L 13 τ η α 1 α 3 L 21 τ η α 2 α 1 L 22 τ η α 2 α 2 L 23 τ η α 2 α 3 L 31 τ η α 3 α 1 L 32 τ η α 3 α 2 L 33 τ η α 3 α 3 = i = 1 3 j = 1 3 D Ω * α i A τ η α i α j D Ω α j     for τ , η = 0 , , N + 1
As shown in Equations (33) and (34), the application of the integration by parts rule with respect to the integration along α 1 , α 2 allows one to also derive the boundary conditions of the problem, which are applied at the edges of the rectangular physical domain α 1 0 , α 1 1 × α 2 0 , α 2 1 of extremes α i 0 , α i 1 with i = 1 , 2 . More specifically, the boundary conditions reported below are applied at the edges of the shell located at α 1 = α 1 0 or α 1 = α 1 1 :
N 1 τ α 1 = N ¯ 1 τ α 1 or u 1 τ = u ¯ 1 τ N 12 τ α 2 = N ¯ 12 τ α 2 or u 2 τ = u ¯ 2 τ T 1 τ α 3 = T ¯ 1 τ α 3 or u 3 τ = u ¯ 3 τ
where N ¯ 1 τ α 1 , N ¯ 12 τ α 2 , T ¯ 1 τ α 3 and u ¯ 1 τ , u ¯ 2 τ , u ¯ 3 τ are pre-determined values of the generalized stress resultants and the generalized displacement components, respectively, which can be assigned a-priori. In the same way, the following boundary conditions are derived for the physical domain edges with α 2 = α 2 0 or α 2 = α 2 1 :
N 21 τ α 1 = N ¯ 21 τ α 1 or u 1 τ = u ¯ 1 τ N 2 τ α 2 = N ¯ 2 τ α 2 or u 2 τ = u ¯ 2 τ T 2 τ α 3 = T ¯ 2 τ α 3 or u 3 τ = u ¯ 3 τ
In the previous relation, the fixed values of the generalized stress resultants and generalized displacement field components are denoted by N ¯ 21 τ α 1 , N ¯ 2 τ α 2 , T ¯ 2 τ α 3 and u ¯ 1 τ , u ¯ 2 τ , u ¯ 3 τ , respectively.
Starting from Equations (46) and (47), the boundary conditions of physical interests are derived because they assign a null value to the kinematic and static quantities along the shell sides. In particular, the Simply-supported (S) boundary conditions are defined in order to neglect in-plane displacements of the lateral surfaces of the doubly-curved shell solid:
N 1 τ α 1 = 0 , u 2 τ = u 3 τ = 0 at α 1 = α 1 0 or α 1 = α 1 1 N 2 τ α 2 = 0 , u 1 τ = u 3 τ = 0 at α 2 = α 2 0 or α 2 = α 2 1

6. Semi-Analytical Navier Solution

In the present section, a semi-analytical solution is found for the higher-order differential problem of Equation (44). To this end, the Navier method is adopted; therefore, some geometric assumptions are made. More specifically, it is assumed that the geometry of the shell is characterized by constant values of the Lamè parameters A 1 , A 2 and of the principal radii of curvature R 1 , R 2 . As a result, the relations reported in the following are considered:
A 1 = 1 n + m A 1 s 1 n s 2 m = 0 , A 2 = 1 n + m A 2 s 1 n s 2 m = 0 , R 1 = cos t n + m R 1 s 1 n s 2 m = 0 , R 2 = cos t n + m R 2 s 1 n s 2 m = 0
As a consequence, the lengths L 1 , L 2 of the curvilinear parametric lines can be calculated in terms of the radii R 1 , R 2 as follows:
L 1 = s 1 1 s 1 0 = α 1 1 α 1 0 R 1 L 2 = s 2 1 s 2 0 = α 2 1 α 2 0 R 2
In the case of a cylindrical surface with k n 1 = 0 and R 2 = R , the quantities L 1 , L 2 read as follows:
L 1 = s 1 1 s 1 0 = α 1 1 α 1 0 L 2 = s 2 1 s 2 0 = α 2 1 α 2 0 R
When a rectangular plate is studied, the principal curvatures are null, namely k n 1 = k n 2 = 0 , therefore, the lengths of the parametric lines are calculated from the following expression:
L 1 = s 1 1 s 1 0 = α 1 1 α 1 0 L 2 = s 2 1 s 2 0 = α 2 1 α 2 0
The semi-analytical Navier solution accounts for the harmonic distribution of the unknown field variable within the physical domain s 1 , s 2 0 , L 1 × 0 , L 2 . The displacement field components u 1 τ , u 2 τ and u 3 τ are thus expressed for each τ -th kinematic expansion order with τ = 0 , , N + 1 as follows:
u 1 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ U 1 n m τ cos n π L 1 s 1 sin m π L 2 s 2 u 2 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ U 2 n m τ sin n π L 1 s 1 cos m π L 2 s 2 u 3 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ U 3 n m τ sin n π L 1 s 1 sin m π L 2 s 2
In the previous equation, the quantities n and m denote the wave number of the solution along s 1 and s 2 , respectively, whereas the quantities U 1 n m τ , U 2 n m τ and U 3 n m τ are the wave amplitudes associated to each wave number n , m . The total number of waves along s 1 and s 2 is denoted by N ˜ and M ˜ , respectively. According to the Navier’s method, these quantities are assumed as N ˜ = M ˜ = . It can be seen that the harmonic expansion of Equation (53) respects the following boundary conditions along the edges of the physical domain:
u 1 τ s 1 = 0 , s 2 0 u 1 τ s 1 = L 1 , s 2 0 u 1 τ s 1 , s 2 = 0 = 0 u 1 τ s 1 , s 2 = L 2 = 0 u 2 τ s 1 = 0 , s 2 = 0 u 2 τ s 1 = L 1 , s 2 = 0 u 2 τ s 1 , s 2 = 0 0 u 2 τ s 1 , s 2 = L 2 0 u 3 τ s 1 = 0 , s 2 = 0 u 3 τ s 1 = L 1 , s 2 = 0 u 3 τ s 1 , s 2 = 0 = 0 u 3 τ s 1 , s 2 = L 2 = 0
As far as the external loads are concerned, the quantities q 1 ± , q 2 ± and q 3 ± of Equation (39), which are applied at the top and bottom surfaces of the shell, are also expanded in a harmonic form by means of the following expression:
q 1 ± s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 1 s n m ± cos n π L 1 s 1 sin m π L 2 s 2 q 2 ± s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 2 s n m ± sin n π L 1 s 1 cos m π L 2 s 2 q 3 ± s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 3 s n m ± sin n π L 1 s 1 sin m π L 2 s 2
with N ˜ = M ˜ = . The wave amplitudes of the external loads, defined for each wave number n , m , are denoted by Q 1 s n m ± , Q 2 s n m ± and Q 3 s n m ± . In Figure 1 and Figure 2 we report the expressions of the wave amplitudes of some load shapes of practical interest.
The generalized external actions q 1 τ , q 2 τ and q 3 τ introduced in Equation (38) which are on the shell reference surface, are expanded for each τ = 0 , , N + 1 with trigonometric series:
q 1 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 1 s n m τ cos n π L 1 s 1 sin m π L 2 s 2 q 2 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 2 s n m τ sin n π L 1 s 1 cos m π L 2 s 2 q 3 τ s 1 , s 2 = n = 1 N ˜ m = 1 M ˜ Q 3 s n m τ sin n π L 1 s 1 sin m π L 2 s 2
being Q 1 s n m τ , Q 2 s n m τ and Q 3 s n m τ the amplitudes of the generalized external actions q 1 τ , q 2 τ , q 3 τ associated to the τ -th kinematic expansion order, with τ = 0 , , N + 1 . Introducing Equations (55) and (56) in the static equivalence principle of Equation (37), the following definitions of the generalized amplitudes Q 1 s n m τ , Q 2 s n m τ , Q 3 s n m τ of the actions q 1 τ , q 2 τ , q 3 τ are obtained in terms of the wave amplitudes Q 1 s n m ± , Q 2 s n m ± , Q 3 s n m ± of the surface loads applied at the top and bottom surfaces of the shell:
Q 1 s n m τ = Q 1 s n m F τ 1 α 1 H 1 H 2 + Q 1 s n m + F τ l α 1 + H 1 + H 2 + Q 2 s n m τ = Q 2 s n m F τ 1 α 2 H 1 H 2 + Q 2 s n m + F τ l α 2 + H 1 + H 2 + Q 3 s n m τ = Q 3 s n m F τ 1 α 3 H 1 H 2 + Q 3 s n m + F τ l α 3 + H 1 + H 2 +
As it is well-known, the semi-analytical Navier solution can be found only for cross-ply lamination schemes. Therefore, the three-dimensional elastic constitutive relationship of Equation (17), written in the reference system of the problem for a generally anisotropic material, takes the following aspect:
σ 1 k σ 2 k τ 12 k τ 13 k τ 23 k σ 3 k = E ¯ 11 k E ¯ 12 k 0 0 0 E ¯ 13 k E ¯ 12 k E ¯ 22 k 0 0 0 E ¯ 23 k 0 0 E ¯ 66 k 0 0 0 0 0 0 E ¯ 44 k 0 0 0 0 0 0 E ¯ 55 k 0 E ¯ 13 k E ¯ 23 k 0 0 0 E ¯ 33 k ε 1 k ε 2 k γ 12 k γ 13 k γ 23 k ε 3 k           k = 1 , , l
In addition, the material orientation angle ϑ k occurring in Equation (19) is selected so that ϑ k = ± π / 2 or ϑ k = 0 . Introducing the harmonic expansion of the unknown field variables of Equation (53) and of the generalized external actions of Equation (56) in the fundamental relations of Equation (44), the vector U n m τ = U 1 n m τ U 2 n m τ U 3 n m τ T of the wave amplitude of the displacement field components is derived for each τ = 0 , , N + 1 from the following expression:
n = 1 N ˜ m = 1 M ˜ η = 0 N + 1 L ˜ n m τ η U n m η + Q s n m τ = 0
where the vector Q s n m τ = Q 1 s n m τ Q 2 s n m τ Q 3 s n m τ T collects the amplitudes of the generalized external actions of the τ -th expansion order. Furthermore, the quantity L ˜ n m τ η is the fundamental matrix related to the wave numbers n , m , defined for each τ , η = 0 , , N + 1 . In a more expanded form, the previous relation becomes:
n = 1 N ˜ m = 1 M ˜ η = 0 N + 1 L ˜ 11 n m τ η α 1 α 1 L ˜ 12 n m τ η α 1 α 2 L ˜ 13 n m τ η α 1 α 3 L ˜ 14 n m τ η α 1 α 4 L ˜ 15 n m τ η α 1 α 5 L ˜ 16 n m τ η α 1 α 6 L ˜ 17 n m τ η α 1 α 7 L ˜ 21 n m τ η α 2 α 1 L ˜ 22 n m τ η α 2 α 2 L ˜ 23 n m τ η α 2 α 3 L ˜ 24 n m τ η α 2 α 4 L ˜ 25 n m τ η α 2 α 5 L ˜ 26 n m τ η α 2 α 6 L ˜ 27 n m τ η α 2 α 7 L ˜ 31 n m τ η α 3 α 1 L ˜ 32 n m τ η α 3 α 2 L ˜ 33 n m τ η α 3 α 3 L ˜ 34 n m τ η α 3 α 4 L ˜ 35 n m τ η α 3 α 5 L ˜ 36 n m τ η α 3 α 6 L ˜ 37 n m τ η α 3 α 7 L ˜ 41 n m τ η α 4 α 1 L ˜ 42 n m τ η α 4 α 2 L ˜ 43 n m τ η α 4 α 3 L ˜ 44 n m τ η α 4 α 4 L ˜ 45 n m τ η α 4 α 5 L ˜ 46 n m τ η α 4 α 6 L ˜ 47 n m τ η α 4 α 7 L ˜ 51 n m τ η α 5 α 1 L ˜ 52 n m τ η α 5 α 2 L ˜ 53 n m τ η α 5 α 3 L ˜ 54 n m τ η α 5 α 4 L ˜ 55 n m τ η α 5 α 5 L ˜ 56 n m τ η α 5 α 6 L ˜ 57 n m τ η α 5 α 7 L ˜ 61 n m τ η α 6 α 1 L ˜ 62 n m τ η α 6 α 2 L ˜ 63 n m τ η α 6 α 3 L ˜ 64 n m τ η α 6 α 4 L ˜ 65 n m τ η α 6 α 5 L ˜ 66 n m τ η α 6 α 6 L ˜ 67 n m τ η α 6 α 7 L ˜ 71 n m τ η α 7 α 1 L ˜ 72 n m τ η α 7 α 2 L ˜ 73 n m τ η α 7 α 3 L ˜ 74 n m τ η α 7 α 4 L ˜ 75 n m τ η α 7 α 5 L ˜ 76 n m τ η α 7 α 6 L ˜ 77 n m τ η α 7 α 7 U 1 n m η U 2 n m η U 3 n m η Φ n m η Ψ n m η Ξ n m η Κ n m η + Q 1 s n m τ Q 2 s n m τ Q 3 s n m τ Q D s n m τ Q B s n m τ Q T s n m τ Q C s n m τ = 0 0 0 0 0 0 0
The complete expression of the coefficients L ˜ i j n m τ η α i α j of the fundamental matrix L ˜ n m τ η , with i , j = 1 , 2 , 3 , is reported below for each τ , η = 0 , , N + 1 :
L ˜ 11 n m τ η α 1 α 1 = A 11 20 τ η 00 α 1 α 1 n π L 1 2 A 66 02 τ η 00 α 1 α 1 m π L 2 2 A 44 20 τ η 00 α 1 α 1 R 1 2 + A 44 10 τ η 01 α 1 α 1 + A 44 10 τ η 10 α 1 α 1 R 1 A 44 00 τ η 11 α 1 α 1 L ˜ 12 n m τ η α 1 α 2 = A 12 11 τ η 00 α 1 α 2 + A 66 11 τ η 00 α 1 α 2 n π L 1 m π L 2 L ˜ 13 n m τ η α 1 α 3 = A 13 10 τ η 01 α 1 α 3 A 44 10 τ η 10 α 1 α 3 + A 11 20 τ η 00 α 1 α 3 + A 44 20 τ η 00 α 1 α 3 R 1 + A 12 11 τ η 00 α 1 α 3 R 2 n π L 1 L ˜ 21 n m τ η α 2 α 1 = A 12 11 τ η 00 α 2 α 1 + A 66 11 τ η 00 α 2 α 1 n π L 1 m π L 2 L ˜ 22 n m τ η α 2 α 2 = A 66 20 τ η 00 α 2 α 2 n π L 1 2 A 22 02 τ η 00 α 2 α 2 m π L 2 2 A 55 02 τ η 00 α 2 α 2 R 2 2 + A 55 01 τ η 01 α 2 α 2 + A 55 01 τ η 10 α 2 α 2 R 2 A 55 00 τ η 11 α 2 α 2 L ˜ 31 n m τ η α 3 α 1 = A 13 10 τ η 10 α 3 α 1 A 44 10 τ η 01 α 3 α 1 + A 11 20 τ η 00 α 3 α 1 + A 44 20 τ η 00 α 3 α 1 R 1 + A 12 11 τ η 00 α 3 α 1 R 2 n π L 1 L ˜ 32 n m τ η α 3 α 2 = A 23 01 τ η 10 α 3 α 2 A 55 01 τ η 01 α 3 α 2 + A 12 11 τ η 00 α 3 α 2 R 1 + A 22 02 τ η 00 α 3 α 2 + A 55 02 τ η 00 α 3 α 2 R 2 m π L 2 L ˜ 33 n m τ η α 3 α 3 = A 44 20 τ η 00 α 3 α 3 n π L 1 2 A 55 02 τ η 00 α 3 α 3 m π L 2 2 A 11 20 τ η 00 α 3 α 3 R 1 2 A 22 02 τ η 00 α 3 α 3 R 2 2 2 A 12 11 τ η 00 α 3 α 3 R 1 R 2 A 13 10 τ η 01 α 3 α 3 + A 13 10 τ η 10 α 3 α 3 R 1 A 23 01 τ η 01 α 3 α 3 + A 23 01 τ η 10 α 3 α 3 R 2 A 33 00 τ η 11 α 3 α 3
It should be noted that in the case of a cylindrical surface, the fundamental governing equations are modified, remembering the geometric relations of Equation (51). More specifically, when R 1 = , the generalized coefficients A n m p q τ η f g α i α j of Equation (29) are calculated for each τ , η = 0 , , N + 1 from the following relation:
A n m p q τ η f g α i α j = k = 1 l ζ k ζ k + 1 B ¯ n m k f F k η α j ζ f g F k τ α i ζ g H 2 H 2 q d ζ
In this case, the fundamental coefficients L ˜ i j n m τ η α i α j , which have been defined in Equation (61), are simplified because one principal direction is characterized by a null value of the principal curvature. The interested reader can find in Appendix A the extended expression of L ˜ i j n m τ η α i α j for the case of a cylindrical surface.
Finally, in the case of a rectangular plate with R 1 = R 2 = , the generalized coefficients A n m p q τ η f g α i α j are derived as follows:
A n m p q τ η f g α i α j = k = 1 l ζ k ζ k + 1 B ¯ n m k f F k η α j ζ f g F k τ α i ζ g d ζ
The expression of the coefficients L ˜ i j n m τ η α i α j of the fundamental matrix L ˜ n m τ η can be found for the case of a rectangular plate in Appendix B.

7. Generalized Differential Quadrature Method

In the present section, the main features of the GDQ method are presented for the derivation of the numerical solution of the three-dimensional equilibrium equations along the thickness direction, as shown previously.
A computational grid made of a discrete number of points is defined in the thickness direction, following a symmetric non-uniform distribution. In this context, the Chebyshev–Gauss–Lobatto (CGL) distribution [16] is adopted. Referring to the interval 1 , 1 , the location of the generic element x ¯ i of the grid at issue is derived as follows:
x ¯ i = cos i 1 I Q 1 π , i = 1 , , I Q , for x i 1 , 1
where I Q is the total number of discrete points. As stated previously, the GDQ method allows one to evaluate the derivative of a generic n -th order of a smooth function from a weighted sum of the values assumed by the function itself in its definition domain. Referring to an arbitrary univariate function f = f x defined in the closed interval a , b with a , b , its n -th order derivative in a point x i a , b with i = 1 , , I Q is thus calculated with the expression reported in the following [16]:
f n x i = n f x x n x = x i j = 1 I Q ς i j n f x j       i = 1 ,   2 , ,   I Q
In the previous relation, the quantity f x j with j = 1 , , I Q denotes the values assumed by the function in the discrete grid, whereas ς i j n are the weighting coefficients of the numerical method. As shown in other works, the present numerical approach provides a high level of accuracy for a sufficient number of discrete points, namely I Q > n . The GDQ weighting coefficients ς i j n are calculated with a recursive procedure [16] based on the adoption of the Lagrange polynomials, defined on the computational discrete grid, for the interpolation of the solution:
ς i j 1 = L 1 x i x i x j L 1 x j , ς i j ( n ) = n ς i j 1 ς i i n 1 ς i j n 1 x i x j i j ς i i n = j = 1   j i I N ς i j n i = j
In the previous relation, the quantities L 1 x i and L 1 x j denote the first order derivatives of the Lagrange polynomials at the points x i and x j , respectively. On the other hand, the definition ς i j 0 = δ i j with i , j = 1 , , I Q should be introduced, being δ i j the Kronecker delta operator.
The GDQ method can also be applied for the numerical computation of integrals within the GIQ numerical method. According to the GIQ, the integration, restricted to the closed interval x i , x j with x i , x j a , b , of a smooth function f = f x with x a , b can be evaluated as follows [16]:
x i x j f x d x = k = 1 I Q w ˜ k i j f x k = k = 1 I Q w j k w i k f x k
As can be seen, the definition interval of f is discretized with a grid of I Q points. The GIQ coefficients w i k , w j k with i , j , k = 1 , , I Q are collected in the matrix W of size I Q × I Q . The coefficients at issue are derived from the GDQ shifted coefficients of the first order derivative, denoted by ς ¯ i j 1 with i , j = 1 , I Q , which are defined as follows, setting ε = 1 × 10 10 :
ς ¯ i j 1 = r i ε r j ε ς ˜ i j 1 i j 1 r i ε + ς ˜ i j 1 i = j
The shifted coefficients of Equation (68) are collected in the matrix ς ¯ 1 , whose size is I Q × I Q . It can be shown that the matrix of the GIQ coefficients is the inverse of the matrix ς ¯ 1 :
W = ς ¯ 1 1
When the numerical integration is restricted to an arbitrary interval a , b , the domain 1 , 1 becomes a parent interval. For this reason, the coefficients w ˜ k 1 I Q of Equation (67) are transformed into those w k 1 I Q by means of the following GIQ mapping technique:
w k 1 I Q = b a 2 w ˜ k 1 I Q
In this way, the integral of f = f x restricted to the interval a , b are evaluated as follows [16]:
a b f ( x ) d x = k = 1 I Q w k 1 I Q f ( x k )

8. Stress and Strain Recovery Procedure

In the previous section, the two-dimensional Navier closed-form solution of the fundamental relations reported in Equation (44) was derived. Therefore, the actual response of the three-dimensional doubly-curved solid is now derived. The reconstruction of stress and strain profiles requires the adoption of three-dimensional equilibrium equations because only the adoption of the ESL kinematic and constitutive relations may lead to erroneous results. Referring to a doubly-curved shell solid with constant principal radii of curvature R 1 , R 2 along the physical domain and A 1 , A 2 = 1 , the three-dimensional equilibrium equations assume the following aspect [16]:
τ 13 k ζ + 2 R 1 + ζ + 1 R 2 + ζ τ 13 k = a k τ 23 k ζ + 1 R 1 + ζ + 2 R 2 + ζ τ 23 k = b k σ 3 k ζ + 1 R 1 + ζ + 1 R 2 + ζ σ 3 k = c k
where the parameters a k , b k and c k are written in an extended form as follows:
a k = 1 1 + ζ / R 1 σ 1 k s 1 1 1 + ζ / R 2 τ 12 k s 2 b k = 1 1 + ζ / R 1 τ 12 k s 1 1 1 + ζ / R 2 σ 2 k s 2 c k = σ 1 k R 1 + ζ + σ 2 k R 2 + ζ 1 1 + ζ / R 1 τ 13 k s 1 1 1 + ζ / R 2 τ 23 k s 2
At this point, a two-dimensional grid is extracted from the physical domain s 1 0 , s 1 1 × s 2 0 , s 2 1 made of I N × I M nodes, starting from the non-uniform one-dimensional CGL distribution of Equation (64).
As far as the thickness direction is concerned, a discrete grid of I T points is defined in each interval ζ k , ζ k + 1 of length h k = ζ k + 1 ζ k related to an arbitrary k -th layer of the stacking sequence, setting k = 1 , , l . The generic element of this grid is denoted by ζ m ˜ k for m ˜ = 1 , , I T , with ζ m ˜ k ζ k , ζ k + 1 . The quantity ζ m ˜ k is defined from an arbitrary distribution within the dimensionless interval ξ m ˜ 1 , 1 as follows:
ζ m ˜ k = h k 2 ξ m ˜ + ζ k + 1 + ζ k 2
Finally, the elements ζ m ˜ k are arranged in the vector ζ k = ζ 1 k ζ m ˜ k ζ I T k T of size I T × 1 , defined in each k -th lamina. At this point, a new vector of size l I T × 1 is introduced, which contains all the discrete points ζ m belonging to the interval h / 2 , h / 2 that are located in the thickness direction. To this end, the index m = 1 , , l I T is introduced, defined as m = k 1 I T + m ˜ . As a consequence, the vector ζ k is arranged in the global vector ζ 1 ζ m ζ l I T T = ζ 1 T ζ k T ζ l T of size l I T × 1 .
The through-the-thickness displacement field profile can be evaluated for each point s 1 i , s 2 j of the reference surface, remembering the kinematic assumption of Equation (8). The relation reported in the following is thus considered for each m = 1 , , l I T :
U i j m k = τ = 0 N + 1 F τ i j m k u τ s 1 i , s 2 j
In the same way, from Equation (14) the profiles are derived of the three-dimensional strain components of the vector ε i j m k = ε 1 i j m k ε 2 i j m k γ 12 i j m k γ 13 i j m k γ 23 i j m k ε 3 i j m k T for each point s 1 i , s 2 j of the reference surface of the shell:
ε i j m k = τ = 0 N + 1 i = 1 3 Z i j m k τ α i ε i j τ α i
The quantity ε i j τ α i = ε τ α i s 1 i , s 2 j with i = 1 , 2 , 3 denotes the generalized strain vector, defined for each τ -th kinematic expansion order, evaluated in each point of the reference surface. Starting from the three-dimensional constitutive relation of Equation (17) with E ¯ i j k = C ¯ i j k for i , j = 1 , , 6 , the distribution of the membrane stresses σ 1 k , σ 2 k and τ 12 k is derived introducing in each through-the-thickness discrete point of the computational grid the three-dimensional strain components of Equation (76), remembering the hypotheses made for the derivation of the Navier closed-form solution:
σ 1 i j m k σ 2 i j m k τ 12 i j m k = C ¯ 11 m k C ¯ 12 m k 0 0 0 C ¯ 13 m k C ¯ 12 m k C ¯ 22 m k 0 0 0 C ¯ 23 m k 0 0 C ¯ 66 m k 0 0 0 ε 1 i j m k ε 2 i j m k γ 12 i j m k γ 13 i j m k γ 23 i j m k ε 3 i j m k
Once the discrete distribution of membrane stresses is derived along the shell thickness according to Equation (77), it is useful to compute their first order derivative with respect to s 1 , s 2 parametric lines. These derivations are performed numerically by means of the GDQ method of Equation (65). One gets:
σ 1 , 1 k = σ 1 k s 1 i j m r = 1 I N ς i r s 1 1 σ 1 r j m k ,     σ 2 , 2 k = σ 2 k s 2 i j m r = 1 I M ς j r s 2 1 σ 2 i r m k , τ 12 , 1 k = τ 12 k s 1 i j m r = 1 I N ς i r s 1 1 τ 12 r j m k ,     τ 12 , 2 k = τ 12 k s 2 i j m r = 1 I M ς j r s 2 1 τ 12 i r m k
where the notations ς i r 1 = ς i r s 1 1 and ς j r 1 = ς j r s 2 1 mean that the GDQ rule is applied along s 1 and s 2 parametric line, respectively. Note that the partial derivatives of the stresses reported in Equation (78) are not solved according to the semi-analytical procedure, because in this case they should be evaluated for each n , m of the trigonometric expansion of Equation (53), and the post-processing recovery procedure should be applied many times. In contrast, the adoption of the GDQ numerical method allows one to apply directly the procedure to the expanded solution. The out-of-plane stress components τ 13 k and τ 23 k are evaluated from the first two three-dimensional equilibrium Equation (72) of a doubly-curved shell, which are written in each k -th layer of the shell in a discrete form as follows:
r = 1 I T ς m ˜ r ζ 1 τ ¯ 13 i j k 1 I T + r k + 2 R 1 + ζ m + 1 R 2 + ζ m τ ¯ 13 i j m k = a i j m k r = 1 I T ς m ˜ r ζ 1 τ ¯ 23 i j k 1 I T + r k + 1 R 1 + ζ m + 2 R 2 + ζ m τ ¯ 23 i j m k = b i j m k
The parameters a k and b k , dependent on the membrane stresses σ 1 k , σ 2 k and τ 12 k , are written in a discrete form as follows:
a i j m k = 1 1 + ζ m / R 1 σ 1 k s 1 i j m 1 1 + ζ m / R 2 τ 12 k s 2 i j m b i j m k = 1 1 + ζ m / R 1 τ 12 k s 1 i j m 1 1 + ζ m / R 2 σ 2 k s 2 i j m
The solution of Equation (79) is derived numerically by means of the GDQ method. On the other hand, the boundary conditions are enforced at the bottom surface of the shell, remembering the reciprocity principle of stress components. In other words, the out-of-plane shear stresses must be at the bottom of the shell, equal to the in-plane loads q 1 and q 2 . In the same way, at the top surface of the shell, the shear stress profile must guarantee equilibrium with the external loads q 1 + and q 2 + :
τ 13 1 ζ = h / 2 = q 1 or τ 13 l ζ = h / 2 = q 1 + τ 23 1 ζ = h / 2 = q 2 or τ 23 l ζ = h / 2 = q 2 +
The boundary conditions introduced in the previous equation are written below in discrete form:
τ ¯ 13 i j 1 1 = q 1 i j , τ ¯ 23 i j 1 1 = q 2 i j
τ 13 i j l I T l = q 1 i j + , τ 23 i j l I T l = q 2 i j +
At this point, a numerical solution is found for k = 1 , taking into account the external constraints of Equation (82). Once the equilibrium Equation (79) is solved in the first layer, the boundary conditions for the generic k -th layer with k = 2 , , l are assessed starting from the results obtained for k 1 . Note that the discrete points associated to the indexes m = k 1 I T and m = k 1 I T + 1 are located at the same height in the interface region along the thickness direction, therefore, the relations reported below can be enforced at the interface between two adjacent layers:
τ ¯ 13 i j k 1 I T + 1 k = τ ¯ 13 i j k 1 I T k 1 τ ¯ 23 i j k 1 I T + 1 k = τ ¯ 23 i j k 1 I T k 1
Finally, the boundary conditions of Equation (83) are enforced at the top surface of the solid, remembering that the solution of the differential system of Equation (79) is defined as less than a linear transformation. For this reason, the through-the-thickness profiles of the stresses τ ¯ 13 k and τ ¯ 23 k , derived numerically, are corrected by adding a linear term dependent on the thickness coordinate ζ m [16], as shown in Figure 3:
τ 13 i j m k = τ ¯ 13 i j m k + q 1 i j + τ ¯ 13 i j l I T l h ζ m + h 2 τ 23 i j m k = τ ¯ 23 i j m k + q 2 i j + τ ¯ 23 i j l I T l h ζ m + h 2 m = 1 , , l I T
At this point, the derivative of τ 13 k and τ 23 k shear stresses with respect to s 1 , s 2 can be evaluated with the GDQ method as follows:
τ 13 , 1 k = τ 13 k s 1 i j m r = 1 I N ς i r s 1 1 τ 13 r j m k τ 23 , 2 k = τ 23 k s 2 i j m r = 1 I M ς j r s 2 1 τ 23 i r m k m = 1 , , l I T
The third equilibrium equation of Equation (72), reported below in discrete form, is adopted for the derivation of the actual profile of the normal stress σ 3 k :
r = 1 I   T ς m ˜ r ζ 1 σ ¯ 3 i j k 1 I T + r k + 1 R 1 + ζ m + 1 R 2 + ζ m σ ¯ 3 i j m k = c i j m k
where the parameter c i j m k accounts for the following expression:
c i j m k = σ 1 i j m k R 1 + ζ m + σ 2 i j m k R 2 + ζ m 1 1 + ζ m / R 1 τ 13 k s 1 i j m 1 1 + ζ m / R 2 τ 23 k s 2 i j m
Two possible sets of boundary conditions are considered for the derivation of the numerical solution of Equation (87), setting q 3 + and q 3 the value of the external actions oriented along the normal direction applied at the top and bottom surfaces of the shell, respectively:
σ ¯ 3 1 ζ = h / 2 = q 3 σ ¯ 3 l ζ = h / 2 = q 3 +
In discrete form, the last two relations become:
σ ¯ 3 i j 1 1 = q 3 i j σ ¯ 3 i j l I T l = q 3 i j +
Following the same approach as Equation (84), once the boundary condition in Equation (89) is employed for the derivation of the numerical solution of Equation (87) of the first layer k = 1 , the following boundary condition is considered for the numerical assessment of the equilibrium Equation (87) in the remaining laminae, namely for k = 2 , , l :
σ ¯ 3 i j k 1 I T + 1 k = σ ¯ 3 i j k 1 I T k 1
The second boundary condition of Equation (90) is applied by means of a linear correction [16] of the profile of the normal stress σ ¯ 3 k derived from Equation (87):
σ 3 i j m k = σ ¯ 3 i j m k + q 3 i j + σ ¯ 3 i j l I T l h ζ m + h 2
The correction of the through-the-thickness profile of the normal stress has been represented in Figure 3. Once the three-dimensional stresses are derived from the present recovery procedure, the constitutive relationship of Equation (17) is used for the derivation of the updated profile of the out-of-plane strain components γ 13 k , γ 23 k and ε 3 k . To this end, the following relation is considered at each point s 1 i , s 2 j , ζ m of the three-dimensional solid:
C ¯ 44 m k γ 13 i j m k = τ 13 i j m k C ¯ 55 m k γ 23 i j m k = τ 23 i j m k C ¯ 33 m k ε 3 i j m k = σ 3 i j m k C ¯ 13 m k ε 1 i j m k C ¯ 23 m k ε 2 i j m k
Starting from the previous equation, the expression of the quantities γ 13 i j m k , γ 23 i j m k and ε 3 i j m k is easily derived.
Finally, the out-of-plane strain components are introduced in the three-dimensional constitutive relation of Equation (77). In this way, the membrane stresses σ 1 i j m k , σ 2 i j m k and τ 12 i j m k are corrected because, in the previous step, they were calculated without considering the three-dimensional equilibrium equations.

9. Application and Results

In the present section, some examples of investigation are shown, where the present semi-analytical theory is adopted for the derivation of the static response of structures of different curvatures and lamination schemes for different load cases. The accuracy of the solution is validated with success with respect to highly computationally demanding three-dimensional finite element models, shown in Figure 4, as implemented in the ABAQUS code. In all examples, it is shown that the present semi-analytical solution, together with the recovery of stresses and strains, can accurately predict the three-dimensional response of curved and laminated structures with a reduced computational cost. Furthermore, the convergence of the method is studied when load shapes of practical interest are considered, like uniform, concentrated, line, and hydrostatic loads. Finally, some examples are presented where curved and layered structures are subjected to generally-shaped pressures. In all the examples, lamination schemes are made starting from layers of graphite-epoxy ρ k = 1450 kg / m 3 , whose engineering constants, obtained from Ref. [101], are reported in the following:
E 1 k = 137.90 GPa G 23 k = 6.21 GPa ν 23 k = 0.49 E 2 k = E 3 k = 8.96 GPa G 12 k = G 13 k = 7.10 GPa ν 12 k = ν 13 k = 0.30
In the previous relation, the quantities ν 12 k , ν 13 k , ν 23 k denote the Poisson’s coefficients of the orthotropic materials, whereas E 1 k , E 2 k , E 3 k and G 12 k , G 13 k , G 23 k are the elastic moduli and the shear moduli, respectively. They are related to the three-dimensional stiffness constants C i j k = E i j k with i , j = 1 , , 6 of Equation (18) according to the procedure extensively detailed in Ref. [16]. The stiffness constants C i j k of soft orthotropic layers, denoted by graphite-epoxy soft5 and graphite-epoxy soft10, are five and ten times softer, respectively, than those of the graphite-epoxy of Equation (94). As a consequence, the engineering constants of the graphite-epoxy soft5, whose density is thus ρ k = 290 kg / m 3 , the following aspect:
E 1 k = 27.58 GPa G 23 k = 1.242 GPa ν 23 k = 0.49 E 2 k = E 3 k = 1.792 GPa G 12 k = G 13 k = 1.42 GPa ν 12 k = ν 13 k = 0.30
In the same way, the mechanical properties of the graphite-epoxy soft10 ρ k = 145 kg / m 3 can be written as:
E 1 k = 13.790 GPa G 23 k = 0.621 GPa ν 23 k = 0.49 E 2 k = E 3 k = 0.896 GPa G 12 k = G 13 k = 0.710 GPa ν 12 k = ν 13 k = 0.30
In each simulation, the lamination scheme consists of three layers of thicknesses h 1 = h 3 = 0.03 m and h 2 = 0.04 m with material orientation angles equal to ϑ 1 = ϑ 3 = 0 and ϑ 2 = π / 2 . Finally, the recovery procedure has been applied setting I T = 31 discrete points in the thickness direction for each k -th layer. Note that the computational cost of the present semi-analytical formulation is here intended as the number of terms, denoted by N ˜ and M ˜ , occurring in the harmonic expansion of the solution according to Equation (53). No details are given on the computational time, whose aspect depends on the machine properties and on the numerical implementation of the linear system (59), which is not the main focus of this work.
The first example consists of a simply-supported rectangular plate of dimensions L 1 = 2 m and L 2 = 1 m made of three layers of graphite-epoxy with properties as in Equation (94) of thicknesses h 1 = h 3 = 0.03 m and h 2 = 0.04 m subjected to two different patch loads, one at the top surface and the other at the bottom surface with magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. The position and shape parameters of the load distribution in hand have been selected so that the external pressure is applied in a quarter of the physical domain. On the other hand, the through-the-thickness distribution of the three-dimensional kinematic and mechanical quantities have been evaluated in a region where the external load is not applied, namely 0.25 L 1 , 0.25 L 2 . In Figure 5 the reader can find some information regarding the convergence rate of the present semi-analytical solution, which is computed by the following percentage error e % :
e % = w N ˜ w FEM w FEM × 100 %
where w N ˜ denotes the vertical deflection of the central point of the three-dimensional solid derived from Equation (75), while w FEM = 0.00150236 m is the corresponding value derived from the finite element simulation.
As visible in Table 1, a very rapid convergence rate is found for a limited number of terms within the harmonic expansion (53). In this way, the results of the simulation are shown to be very stable for the selected wave numbers. Note that the selected value of N ˜ = M ˜ not only depends on the convergence of the vertical deflection to the reference value, but also on the fulfillment of the loading conditions at the top and bottom surfaces of the plate under consideration.
The distributions of the displacement field components, the three-dimensional strain, and the stresses have been reported in Figure 5, Figure 6 and Figure 7, respectively. For each quantity, a 3D FEM model of parabolic C3D20 elements, consisting of 582327 DOFs, has been adopted for the derivation of a reference solution. Furthermore, a two-dimensional numerical solution is derived with the GDQ method, accounting for the FSDT and TSDT theories as well as the EDZ4. The physical domain is discretized with a two-dimensional grid made of I N = I M = 31 discrete points, following the CGL distribution. The present semi-analytical theory has been used with the EDZ4 displacement field assumption, setting N ˜ = M ˜ = 150 . As can be seen from Figure 5, very accurate results are provided in terms of the displacement field components U 1 , U 2 , U 3 . On the other hand, in Figure 6, it is shown that the recovery procedure is determining the correct prediction of the out-of-plane strain components γ 13 , γ 23 , ε 3 , especially in the central layer. Furthermore, very accurate results are obtained for in-plane quantities ε 1 , ε 2 , γ 12 . As far as the three-dimensional stress components are concerned, in Figure 7, it is shown that the reconstruction of out-of-plane stresses from the semi-analytical Navier solution by means of the three-dimensional constitutive relationship of Equation (17) leads to erroneous results, especially for the σ 3 normal stress. On the other hand, when the equilibrium-based recovery procedure is applied, the results coming from the two-dimensional semi-analytical solution, denoted by (E), perfectly match those coming from the 3D FEM simulation.
The next example takes into account the same rectangular panel subjected to hydrostatic loads. More specifically, the first load of magnitude q ¯ 3 + = 7 × 10 5 N / m 2 and directed along α 1 principal direction is applied at the top surface, whereas the second load of magnitude q ¯ 3 = 3 × 10 5 N / m 2 , applied at the bottom surface, is directed along α 2 principal direction. Three different load cases have been considered, denoted by H1, H2, and H3. In the first one, only the top surface is loaded, whereas in the second one, the external load is applied only to the bottom surface. Finally, in the H3 load case, the two hydrostatic loads have been considered together. Thickness plots calculated at the point located at 0.25 L 1 , 0.25 L 2 have been collected in Figure 8, Figure 9 and Figure 10.
For each load configuration, the semi-analytical solution has been calculated with N ˜ = M ˜ = 150 . The results are compared to those of a 3D FEM simulation and to those coming from a GDQ numerical model with classical and higher-order theories, showing a very good agreement between different approaches. When lower order theories like FSDT and TSDT are employed, a slight discrepancy in results can be seen in the thickness plot of U 3 , reported in Figure 8, whereas the EDZ4 theory perfectly matches the three-dimensional predictions of all quantities in each load case. With particular reference to both in-plane and out-of-plane strain and stress profiles in Figure 9 and Figure 10, it should be noted that for each case, the loading conditions are perfectly respected. Furthermore, the solution coming from the H3 load case can be seen as the sum of H1 and H2, due to the additivity property of the linear solutions calculated previously.
At this point, a laminated cylindrical panel of radius R 2 = 1.2 m is considered with L 1 = 2 m and α 2 0 , α 2 1 = π / 3 , π / 3 . Two externally concentrated loads of magnitude q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N are applied at the top and bottom surfaces, respectively. The position parameters are s 10 = 0.5 L 1 and s 20 = 0.5 L 2 , selected so that the structure is loaded in the center of the physical domain. The lamination scheme consists of two external layers of graphite-epoxy, as defined in Equation (94), whereas the central core is made of graphite-epoxy soft10, as defined in Equation (96).
Thickness plots are provided at the point located at 0.25 L 1 , 0.25 L 2 , and they are collected in Figure 11, Figure 12 and Figure 13. A 3D FEM model with 741975 DOFs made of parabolic C3D20 brick elements has been developed, and a three-dimensional reference solution has been provided.
In addition, a numerical solution based on the GDQ numerical technique has been derived, taking into account both the FSDT and the TSDT displacement field assumptions. The zigzag effects are clearly visible in the deflection of the panel because an abrupt variation of the material stiffness occurs between two adjacent laminae. For this reason, a parametric investigation has been conducted with the semi-analytical approach with N ˜ = M ˜ = 500 , and the static response of the cylindrical panel has been evaluated employing various higher-order theories. As can be seen from Figure 11, the exact solutions accounting for the EDZ3 and the EDZ4 higher-order theories perfectly match the predictions of the three-dimensional FEM reference model. Similar considerations can be made for the plots of the three-dimensional strain components in Figure 12. The adoption of a higher-order displacement field is key for the prediction of both in-plane and out-of-plane kinematic quantities. In the same way, the three-dimensional stress components of Figure 13 are well described by the EDZ4 model, but with a significantly reduced computational cost if compared to the 3D-FEM simulation.
The EDZ4 model has been taken as a reference also in the next example, where the same simply-supported cylindrical panel, made of three layers of graphite-epoxy, as defined in Equation (94), has been loaded with a line load of magnitude q ¯ 3 + = 4.58 × 10 3 N / m distributed along α 1 principal direction, located at s 20 = 0.75 L 2 . The thickness plots are evaluated at 0.25 L 1 , 0.25 L 2 . The GDQ numerical solution has been calculated with I N = I M = 31 , accounting for FSDT, TSDT, and EDZ4 theories, whereas the semi-analytical Navier solution is derived with the EDZ4 kinematic assumption setting N ˜ = M ˜ = 500 . Thickness plots of displacement field, strain, and stress components are reported in Figure 14, Figure 15 and Figure 16, respectively. As shown in Figure 14, classical approaches like the FSDT and the TSDT provide a uniform value of U 3 , and the stretching effect predicted by the 3D-FEM model cannot be evaluated. On the other hand, when the EDZ4 theory is adopted, the parabolic distribution of the displacement field components is predicted with a sufficient level of accuracy. The distribution of the three-dimensional strain components, reported in Figure 15, provides refined results with respect to 3D FEM if a higher-order displacement field is considered. It is important to underline that the recovery procedure provides very accurate results because of the external load cases, which are the boundary conditions of the problem; therefore, an accurate distribution of the stress components is derived. Finally, from the results reported in Figure 16, it can be said that the three-dimensional in-plane stress profiles derived from the 3D FEM model are in line with those provided by both the GDQ numerical solution and the present approach.
Next example points out the high convergence rate of the present semi-analytical method in the case of general loads. Let us consider a cylindrical panel of graphite epoxy with an internal region made of graphite-epoxy soft10 subjected to different external pressures applied at the extrados of the shell. More specifically, two hydrostatic loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 4 × 10 5 N / m 2 are considered, which are denoted by H1 and H2, respectively. In addition, a uniform pressure (U) of magnitude q ¯ 3 + = 2 × 10 5 N / m 2 acting along the thickness direction is applied. Five different load cases are considered, accounting for different combinations of the external pressures introduced previously. They are summarized as following:
Case   01 ( C 1 ) H 1 Case   02 ( C 2 ) H 2 Case   03 ( C 3 ) H 1 + H 2 Case   04 ( C 4 ) U Case   05 ( C 5 ) H 1 + H 2 + U
For each load case, the thickness plots are evaluated 0.25 L 1 , 0.25 L 2 and collected in Figure 17, Figure 18 and Figure 19, and a 3D finite element reference solution has been derived with commercial software. In addition, a GDQ-based numerical solution has been evaluated with the EDZ4 theory, which matches the 3D FEM predictions. The semi-analytical solution is evaluated with the EDZ4 displacement field assumption, set N ˜ = M ˜ = 150 in each load case. The results are in line with the reference solution in terms of displacement field components, as can be seen in Figure 17. It should be noted that the selected lamination scheme presents some zigzag effects, especially for case C5. The adoption of Murakami’s zigzag function in the present model allows one to consider the piecewise inclination of the profile of the displacement field components. Similar considerations are made for the three-dimensional strain components of Figure 18, where the 3D FEM reference solution is perfectly predicted in a reduced amount of time by the higher-order semi-analytical model, and a good level of accuracy is also reached in the central region of the structure. In order to reduce the computational effort of the simulation, the results of more complicated load cases like C3 are obtained from the algebraic sum of C1 and C2 simulations. In the same way, load case C5 is derived from the sum of C3 and C4. As a consequence, the results of the simulations referred to in C3 and C5 can be efficiently obtained as an algebraic sum of the profiles of all kinematic and mechanical quantities recovered previously in the corresponding two-dimensional semi-analytical solutions, as shown in Figure 19 in the case of the three-dimensional stress components. The recovery procedure is not applied in C3 and C5 because the equilibrium equations in the thickness direction have already been solved independently in C1, C2, and C4. The 3D FEM numerical predictions are perfectly matched by the semi-analytical model, with a significantly reduced computational cost and time.
The next example presents a shallow spherical panel of radius R = 3 m , whose physical domain is defined so that α 1 0 , α 1 1 = 5 π / 12 , 7 π / 12 and α 2 0 , α 2 1 = π / 9 , π / 9 . Three different lamination schemes are considered made of two external layers of graphite-epoxy (94), whereas the central core consists of graphite-epoxy (94), graphite-epoxy soft5 (95) and graphite-epoxy soft10 (96) for cases C1, C2, and C3, respectively. In the first load case, the panel under consideration is subjected to sinusoidal loads of magnitude q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with N ˜ = M ˜ = 1 . Two reference solutions are provided, developed with 3D finite element solution and a two-dimensional GDQ-based formulation, accounting for the EDZ4 displacement field theory.
The 3D FEM model is made of parabolic C3D20 brick elements with a total number of 293847 DOFs, whereas the 2D-GDQ model is built starting from a two-dimensional CGL grid with I N = I M = 31 . The solution obtained from the semi-analytical simulation is based on the EDZ3 and EDZ4 theories. The thickness plots, reported in Figure 20, Figure 21 and Figure 22, are provided for the point located at 0.25 L 1 , 0.25 L 2 within the physical domain, where the quantities L 1 and L 2 have been defined in Equation (50).
As shown in Figure 20, the reduction of the stiffness of the central core of the panel leads to a typical zigzag profile of the in-plane displacement field components U 1 and U 2 . Furthermore, when the central layer stiffness is reduced, the vertical deflection U 3 increases. Similar considerations can be made for the strain components of Figure 21, where ε 1 , ε 2 and ε 3 assume in the central lamina a non-linear profile in the case of C2 and C3, whereas the value of γ 12 , γ 13 and γ 23 is increased. Furthermore, for all strain components, a good agreement can be seen between the predictions of various numerical approaches and the semi-analytical results.
Referring to the results of Figure 22, the values of both in-plane and out-of-plane stress components derived from both the 3D FEM and the GDQ are predicted with success by the present semi-analytical formulation. Furthermore, the boundary conditions are respected at the top and bottom surfaces.
Once the semi-analytical model of the spherical panel has been validated for the case of sinusoidal loads ( N ˜ = M ˜ = 1 ), the linear static response of the same structure is derived for the case of uniform transverse loads q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 applied at the top and bottom surfaces, respectively. In this case, the results obtained with the present semi-analytical theory have been derived setting N ˜ = M ˜ = 150 , taking into account the EDZ4 higher-order theory. The thickness plots are provided for the point 0.25 L 1 , 0.25 L 2 and collected in Figure 23, Figure 24 and Figure 25.
The reference solution has been calculated with a 3D-FEM model and some GDQ numerical simulations, based on the FSDT and the TSDT kinematic field assumptions. The profiles of the displacement field components in Figure 23 show that the predictions of the reference models can be obtained only when a higher-order displacement field is considered in the semi-analytical model. In fact, in the case of softcore lamination schemes, classical approaches like the FSDT and the TSDT are not consistent, whereas the results provided with the EDZ4 theory match the 3D-FEM predictions. In Figure 24, the through-the-thickness profiles of the three-dimensional strain components are provided. As can be seen, for both hardcore and softcore configurations of the stacking sequence, the present higher-order semi-analytical solution predicts with success the strain profiles provided by the three-dimensional model, even in the central softcore lamina. Furthermore, the presence of the zigzag function allows one to see what happens in the interface region between two adjacent laminae.
As far as the three-dimensional stress components are concerned, Figure 25 shows that for each lamination scheme that has been investigated, the results of the three-dimensional model are well predicted when higher-order thickness functions are considered in the two-dimensional semi-analytical solution.

10. Conclusions

In the present manuscript, an efficient two-dimensional semi-analytical model has been presented for the evaluation of the static response of curved laminated panels subjected to arbitrary loads. The fundamental governing equations have been written according to the ESL approach with a generalized description of the kinematic field variable along the thickness direction. The solution to the semi-analytical problem has been provided with the Navier method, coupled with a post-processing recovery procedure based on the GDQ numerical technique. The model has been employed in some examples of investigation, where the three-dimensional linear static response of panels with different curvatures, lamination schemes, and load shapes has been derived and successfully compared to the numerical predictions coming from refined 3D FEM simulations. It has been shown that when the semi-analytical approach is used for the derivation of the solution, the GDQ-based recovery procedure allows the model to perfectly fulfil the load conditions in each point of the panel. In addition, the adoption of higher-order theories, together with stress and strain recovery, allows for a good level of accuracy in evaluating the three-dimensional behavior of structures with more cross-ply lamination schemes. Finally, a high level of accuracy is also seen in the case of softcore layers when a higher-order two-dimensional theory is considered, thus reducing the computational effort of each simulation. The numerical examples show that the EDZ4 theory is a valid tool for many lamination schemes, and the semi-analytical solution perfectly matches the 3D finite element predictions, especially in the case of rectangular plates and cylindrical panels. It has been shown that when uniformly-distributed patch and hydrostatic loads are applied to the panel, the convergence of results is seen for N ˜ = M ˜ = 150 , while further terms are required in the case of concentrated and line loads, namely N ˜ = M ˜ = 500 . It is shown that this issue does not depend on the geometry or lamination scheme, but only on the applied load shape. The present semi-analytical solution can be a valid alternative to well-established finite element simulations. In addition, among two-dimensional theories, it allows for a rapid and accurate solution of the problem for structures of constant curvatures with cross-ply whose lamination schemes are made of orthotropic materials. For this reason, it can be applied to fiber-reinforced composite materials as well as lattice and honeycomb panels and structures reinforced with dispersed short fibers without significant computational effort if compared to trustworthy numerical models.

Author Contributions

Conceptualization, F.T. and R.D.; Methodology, F.T., M.V. and R.D.; Software, F.T.; Validation, F.T., M.V. and R.D.; Formal analysis, F.T., M.V. and R.D.; Investigation, F.T., M.V. and R.D.; Data curation, R.D.; Writing—original draft, F.T. and M.V.; Writing—review & editing, F.T. and R.D.; Supervision, F.T. and R.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In this Appendix, we provide the complete definition of the fundamental coefficients L ˜ i j n m τ η α i α j with i , j = 1 , 2 , 3 of the matrix L ˜ n m τ η occurring in Equation (60) for the case of a cylindrical panel characterized by R 1 = and R 2 = R . As can be seen, they depend on the indexes, τ , η = 0 , , N + 1 which are the kinematic expansion order of Equation (8), and on the wave numbers n , m introduced in Equation (53) for the derivation of the semi-analytical Navier solution.
L ˜ 11 n m τ η α 1 α 1 = A 11 20 τ η 00 α 1 α 1 n π L 1 2 A 66 02 τ η 00 α 1 α 1 m π L 2 2 A 44 00 τ η 11 α 1 α 1
L ˜ 12 n m τ η α 1 α 2 = A 12 11 τ η 00 α 1 α 2 + A 66 11 τ η 00 α 1 α 2 n π L 1 m π L 2
L ˜ 13 n m τ η α 1 α 3 = A 13 10 τ η 01 α 1 α 3 A 44 10 τ η 10 α 1 α 3 + A 12 11 τ η 00 α 1 α 3 R n π L 1
L ˜ 21 n m τ η α 2 α 1 = A 12 11 τ η 00 α 2 α 1 + A 66 11 τ η 00 α 2 α 1 n π L 1 m π L 2
L ˜ 22 n m τ η α 2 α 2 = A 66 20 τ η 00 α 2 α 2 n π L 1 2 A 22 02 τ η 00 α 2 α 2 m π L 2 2 A 55 02 τ η 00 α 2 α 2 R 2 + A 55 01 τ η 01 α 2 α 2 + A 55 01 τ η 10 α 2 α 2 R A 55 00 τ η 11 α 2 α 2
L ˜ 23 n m τ η α 2 α 3 = A 23 01 τ η 01 α 2 α 3 A 55 01 τ η 10 α 2 α 3 + A 22 02 τ η 00 α 2 α 3 + A 55 02 τ η 00 α 2 α 3 R m π L 2
L ˜ 31 n m τ η α 3 α 1 = A 13 10 τ η 10 α 3 α 1 A 44 10 τ η 01 α 3 α 1 + A 12 11 τ η 00 α 3 α 1 R n π L 1
L ˜ 32 n m τ η α 3 α 2 = A 23 01 τ η 10 α 3 α 2 A 55 01 τ η 01 α 3 α 2 + A 22 02 τ η 00 α 3 α 2 + A 55 02 τ η 00 α 3 α 2 R m π L 2
L ˜ 33 n m τ η α 3 α 3 = A 44 20 τ η 00 α 3 α 3 n π L 1 2 A 55 02 τ η 00 α 3 α 3 m π L 2 2 A 22 02 τ η 00 α 3 α 3 R 2 A 23 01 τ η 01 α 3 α 3 + A 23 01 τ η 10 α 3 α 3 R A 33 00 τ η 11 α 3 α 3

Appendix B

In the following, we report the complete expression of the fundamental coefficients L ˜ i j n m τ η α i α j , for the case of a rectangular plate of dimensions L 1 , L 2 . They can be seen as a particular case of Equation (61), setting R 1 = R 2 = .
L ˜ 11 n m τ η α 1 α 1 = A 11 20 τ η 00 α 1 α 1 n π L 1 2 A 66 02 τ η 00 α 1 α 1 m π L 2 2 A 44 00 τ η 11 α 1 α 1
L ˜ 12 n m τ η α 1 α 2 = A 12 11 τ η 00 α 1 α 2 + A 66 11 τ η 00 α 1 α 2 n π L 1 m π L 2
L ˜ 13 n m τ η α 1 α 3 = A 13 10 τ η 01 α 1 α 3 A 44 10 τ η 10 α 1 α 3 n π L 1
L ˜ 21 n m τ η α 2 α 1 = A 12 11 τ η 00 α 2 α 1 + A 66 11 τ η 00 α 2 α 1 n π L 1 m π L 2
L ˜ 22 n m τ η α 2 α 2 = A 66 20 τ η 00 α 2 α 2 n π L 1 2 A 22 02 τ η 00 α 2 α 2 m π L 2 2 A 55 00 τ η 11 α 2 α 2
L ˜ 23 n m τ η α 2 α 3 = A 23 01 τ η 01 α 2 α 3 A 55 01 τ η 10 α 2 α 3 m π L 2
L ˜ 31 n m τ η α 3 α 1 = A 13 10 τ η 10 α 3 α 1 + A 44 10 τ η 01 α 3 α 1 n π L 1
L ˜ 32 n m τ η α 3 α 2 = A 23 01 τ η 10 α 3 α 2 + A 55 01 τ η 01 α 3 α 2 m π L 2
L ˜ 33 n m τ η α 3 α 3 = A 44 20 τ η 00 α 3 α 3 n π L 1 2 A 55 02 τ η 00 α 3 α 3 m π L 2 2 A 33 00 τ η 11 α 3 α 3

References

  1. Jones, R.M. Mechanics of Composite Materials; CRC Press: Boca Raton, FL, USA, 2018. [Google Scholar]
  2. Vannucci, P.; Verchery, G. Stiffness design of laminates using the polar method. Int. J. Solids Struct. 2001, 38, 9281–9294. [Google Scholar] [CrossRef]
  3. Reddy, J.N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis; CRC Press: Boca Raton, FL, USA, 2003. [Google Scholar]
  4. Tsai, S. Introduction to Composite Materials; Technomic Publishing Company: Lancaster, UK, 1980. [Google Scholar]
  5. Baker, A.A. Composite Materials for Aircraft Structures; AIAA: Blacksburg, VI, USA, 2004. [Google Scholar]
  6. Cottrell, J.A.; Reali, A.; Bazilevs, Y.; Hughes, T.J. Isogeometric analysis of structural vibrations. Comput. Methods Appl. Mech. Eng. 2006, 195, 5257–5296. [Google Scholar] [CrossRef]
  7. Piegl, L.; Tiller, W. The NURBS Book; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1996. [Google Scholar]
  8. Hughes, T.J.; Liu, W.K. Nonlinear finite element analysis of shells: Part I. Three-dimensional shells. Comput. Methods Appl. Mech. Eng. 1981, 26, 331–362. [Google Scholar] [CrossRef]
  9. Hii, A.K.W.; Minera, S.; Groh, R.M.J.; Pirrera, A.; Kawashita, L.F. Three-dimensional stress analyses of complex laminated shells with a variable-kinematics continuum shell element. Compos. Struct. 2019, 229, 111405. [Google Scholar] [CrossRef]
  10. Punera, D.; Kant, T. Free vibration of functionally graded open cylindrical shells based on several refined higher order displacement models. Thin-Walled Struct. 2017, 119, 707–726. [Google Scholar] [CrossRef]
  11. Patel, B.P.; Gupta, S.S.; Loknath, M.S.; Kadu, C.P. Free vibration analysis of functionally graded elliptical cylindrical shells using higher-order theory. Compos. Struct. 2005, 69, 259–270. [Google Scholar] [CrossRef]
  12. Guo, J.; Qin, Z.; Zhang, Y. Nonlinear aerodynamic analysis of functional graded plates using an HSDT-based isogeometric approach. Thin-Walled Struct. 2023, 186, 110658. [Google Scholar] [CrossRef]
  13. Farzam, A.; Hassani, B. A new efficient shear deformation theory for FG plates with in-plane and through-thickness stiffness variations using isogeometric approach. Mech. Adv. Mater. Struct. 2019, 26, 512–525. [Google Scholar] [CrossRef]
  14. Kraus, H. Thin Elastic Shells; Wiley: New York, NY, USA, 1967. [Google Scholar]
  15. Calladine, C.R. Theory of Shell Structures; Cambridge University Press: New York, NY, USA, 1983. [Google Scholar]
  16. Tornabene, F. Hygro-Thermo-Magneto-Electro-Elastic Theory of Anisotropic Doubly-Curved Shells; Esculapio: Bologna, Italy, 2023. [Google Scholar]
  17. Liew, K.M.; Pan, Z.Z.; Zhang, L.W. An overview of layerwise theories for composite laminates and structures: Development, numerical implementation and application. Compos. Struct. 2019, 216, 240–259. [Google Scholar] [CrossRef]
  18. Xavier, P.B.; Chew, C.H.; Lee, K.H. Buckling and vibration of multilayer orthotropic composite shells using a simple higher-order layerwise theory. Int. J. Solids Struct. 1995, 32, 3479–3497. [Google Scholar] [CrossRef]
  19. Reddy, J.N. An evaluation of equivalent-single-layer and layerwise theories of composite laminates. Compos. Struct. 1993, 25, 21–35. [Google Scholar] [CrossRef]
  20. Reddy, J.N.; Robbins, D.H., Jr. Theories and computational models for composite laminates. Appl. Mech. Rev. 1994, 47, 147–169. [Google Scholar] [CrossRef]
  21. Tornabene, F.; Viscoti, M.; Dimitri, R. Generalized higher order layerwise theory for the dynamic study of anisotropic doubly-curved shells with a mapped geometry. Eng. Anal. Bound. Elements 2022, 134, 147–183. [Google Scholar] [CrossRef]
  22. Tornabene, F.; Viscoti, M.; Dimitri, R.; Reddy, J.N. Higher order theories for the vibration study of doubly-curved anisotropic shells with a variable thickness and isogeometric mapped geometry. Compos. Struct. 2021, 267, 113829. [Google Scholar] [CrossRef]
  23. Tornabene, F.; Viscoti, M.; Dimitri, R. General boundary conditions implementation for the static analysis of anisotropic doubly-curved shells resting on a Winkler foundation. Compos. Struct. 2023, 322, 117198. [Google Scholar] [CrossRef]
  24. Aminipour, H.; Janghorban, M.; Civalek, O. Analysis of functionally graded doubly-curved shells with different materials via higher order shear deformation theory. Compos. Struct. 2020, 251, 112645. [Google Scholar] [CrossRef]
  25. Mantari, J.L.; Oktem, A.S.; Soares, C.G. A new trigonometric shear deformation theory for isotropic, laminated composite and sandwich plates. Int. J. Solids Struct. 2012, 49, 43–53. [Google Scholar] [CrossRef]
  26. Ferreira, A.J.M.; Roque, C.M.C.; Jorge, R.M.N. Analysis of composite plates by trigonometric shear deformation theory and multiquadrics. Comput. Struct. 2005, 83, 2225–2237. [Google Scholar] [CrossRef]
  27. Al-Osta, M.A. An exponential-trigonometric quasi-3D HSDT for wave propagation in an exponentially graded plate with microstructural defects. Compos. Struct. 2022, 297, 115984. [Google Scholar] [CrossRef]
  28. Mantari, J.L.; Bonilla, E.M.; Soares, C.G. A new tangential-exponential higher order shear deformation theory for advanced composite plates. Compos. Part B Eng. 2014, 60, 319–328. [Google Scholar] [CrossRef]
  29. Sayyad, A.S.; Ghugal, Y.M. Bending, buckling and free vibration of laminated composite and sandwich beams: A critical review of literature. Compos. Struct. 2017, 171, 486–504. [Google Scholar] [CrossRef]
  30. Gherlone, M. On the use of zigzag functions in equivalent single layer theories for laminated composite and sandwich beams: A comparative study and some observations on external weak layers. J. Appl. Mech. 2013, 80, 061004. [Google Scholar] [CrossRef]
  31. Icardi, U.; Sola, F. Assessment of recent zig-zag theories for laminated and sandwich structures. Compos. Part B Eng. 2016, 97, 26–52. [Google Scholar] [CrossRef]
  32. Murakami, H. Laminated Composite Plate Theory with Improved In-Plane Responses. ASME J. Appl. Mech. 1986, 53, 661–666. [Google Scholar] [CrossRef]
  33. Tessler, A.; di Sciuva, M.; Gherlone, M. A refined zigzag beam theory for composite and sandwich beams. J. Compos. Mater. 2009, 43, 1051–1081. [Google Scholar] [CrossRef]
  34. Iurlaro, L.; Gherlone, M.D.; di Sciuva, M.; Tessler, A. Assessment of the refined zigzag theory for bending, vibration, and buckling of sandwich plates: A comparative study of different theories. Compos. Struct. 2013, 106, 777–792. [Google Scholar] [CrossRef]
  35. Di Sciuva, M.; Sorrenti, M. Bending and free vibration analysis of functionally graded sandwich plates: An assessment of the Refined Zigzag Theory. J. Sandw. Struct. Mater. 2021, 23, 760–802. [Google Scholar] [CrossRef]
  36. Rezaei, A.S.; Saidi, A.R. Application of Carrera Unified Formulation to study the effect of porosity on natural frequencies of thick porous–cellular plates. Compos. Part B Eng. 2016, 91, 361–370. [Google Scholar] [CrossRef]
  37. Reddy, J.N. A generalization of two-dimensional theories of laminated composite plates. Commun. Appl. Numer. Methods 1987, 3, 173–180. [Google Scholar] [CrossRef]
  38. Barbero, E.J.; Reddy, J.N.; Teply, J. An accurate determination of stresses in thick laminates using a generalized plate theory. Int. J. Numer. Methods Eng. 1990, 29, 1–14. [Google Scholar] [CrossRef]
  39. Washizu, K. Variational Methods in Elasticity & Plasticity; Pergamon Press: Oxford, UK, 1975. [Google Scholar]
  40. Wang, C.M.; Reddy, J.N.; Lee, K.H. (Eds.) . Shear Deformable Beams and Plates: Relationships with Classical Solutions; Elsevier: Oxford, UK, 2000. [Google Scholar]
  41. Reddy, J.N.; Wang, C.M. An overview of the relationships between solutions of the classical and shear deformation plate theories. Compos. Sci. Technol. 2000, 60, 2327–2335. [Google Scholar] [CrossRef]
  42. Tornabene, F. Free vibration analysis of functionally graded conical, cylindrical shell and annular plate structures with a four-parameter power-law distribution. Comput. Methods Appl. Mech. Eng. 2009, 198, 2911–2935. [Google Scholar] [CrossRef]
  43. Tornabene, F.; Viscoti, M.; Dimitri, R. Higher order theories for the modal analysis of anisotropic doubly-curved shells with a three-dimensional variation of the material properties. Eng. Anal. Bound. Elem. 2024, 158, 486–519. [Google Scholar] [CrossRef]
  44. Tornabene, F. On the critical speed evaluation of arbitrarily oriented rotating doubly-curved shells made of functionally graded materials. Thin-Walled Struct. 2019, 140, 85–98. [Google Scholar] [CrossRef]
  45. Malekzadeh, P.; Zarei, A.R. Free vibration of quadrilateral laminated plates with carbon nanotube reinforced composite layers. Thin-Walled Struct. 2014, 82, 221–232. [Google Scholar] [CrossRef]
  46. Safaei, B.; Ahmed, N.A.; Fattahi, A.M. Free vibration analysis of polyethylene/CNT plates. Eur. Phys. J. Plus 2019, 134, 271. [Google Scholar] [CrossRef]
  47. Rezaiee-Pajand, M.; Sobhani, E.; Masoodi, A.R. Free vibration analysis of functionally graded hybrid matrix/fiber nanocomposite conical shells using multiscale method. Aerosp. Sci. Technol. 2020, 105, 105998. [Google Scholar] [CrossRef]
  48. Tornabene, F.; Viscoti, M.; Dimitri, R.; Aiello, M.A. Higher order formulations for doubly-curved shell structures with a honeycomb core. Thin-Walled Struct. 2021, 164, 107789. [Google Scholar] [CrossRef]
  49. Tornabene, F.; Viscoti, M.; Dimitri, R.; Aiello, M.A. Higher-order modeling of anisogrid composite lattice structures with complex geometries. Eng. Struct. 2021, 244, 112686. [Google Scholar] [CrossRef]
  50. Demasi, L.; Biagini, G.; Vannucci, F.; Santarpia, E.; Cavallaro, R. Equivalent Single Layer, Zig-Zag, and Layer Wise theories for variable angle tow composites based on the Generalized Unified Formulation. Compos. Struct. 2017, 177, 54–79. [Google Scholar] [CrossRef]
  51. Tornabene, F.; Viscoti, M.; Dimitri, R. Free vibration analysis of laminated doubly-curved shells with arbitrary material orientation distribution employing higher order theories and differential quadrature method. Eng. Anal. Bound. Elements 2023, 152, 397–445. [Google Scholar] [CrossRef]
  52. Montemurro, M.; Catapano, A. A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates. Compos. Struct. 2018, 209, 561–578. [Google Scholar] [CrossRef]
  53. Addou, F.Y.; Meradjah, M.; Bousahla, A.A.; Benachour, A.; Bourada, F.; Tounsi, A.; Mahmoud, S.R. Influences of porosity on dynamic response of FG plates resting on Winkler/Pasternak/Kerr foundation using quasi 3D HSDT. Comput. Concr. 2019, 24, 347–367. [Google Scholar]
  54. Vinyas, M. On frequency response of porous functionally graded magneto-electro-elastic circular and annular plates with different electro-magnetic conditions using HSDT. Compos. Struct. 2020, 240, 112044. [Google Scholar] [CrossRef]
  55. Shojaei, M.F.; Ansari, R. Variational differential quadrature: A technique to simplify numerical analysis of structures. Appl. Math. Model. 2017, 49, 705–738. [Google Scholar] [CrossRef]
  56. Tornabene, F.; Viscoti, M.; Dimitri, R. Static analysis of anisotropic doubly-curved shells with arbitrary geometry and variable thickness resting on a Winkler-Pasternak support and subjected to general loads. Eng. Anal. Bound. Elem. 2022, 140, 618–673. [Google Scholar] [CrossRef]
  57. Tornabene, F.; Viscoti, M.; Dimitri, R.; Rosati, L. Dynamic analysis of anisotropic doubly-curved shells with general boundary conditions, variable thickness and arbitrary shape. Compos. Struct. 2023, 309, 116542. [Google Scholar] [CrossRef]
  58. Liew, K.M.; Teo, T.M.; Han, J.B. Three-dimensional static solutions of rectangular plates by variant differential quadrature method. Int. J. Mech. Sci. 2001, 43, 1611–1628. [Google Scholar] [CrossRef]
  59. Brischetto, S. An exact 3D solution for free vibrations of multilayered cross-ply composite and sandwich plates and shells. Int. J. Appl. Mech. 2014, 6, 1450076. [Google Scholar] [CrossRef]
  60. Brischetto, S.; Torre, R. Exact 3D solutions and finite element 2D models for free vibration analysis of plates and cylinders. Curved Layer. Struct. 2014, 1, 59–92. [Google Scholar] [CrossRef]
  61. Duan, W.H.; Wang, C.M. Exact solutions for axisymmetric bending of micro/nanoscale circular plates based on nonlocal plate theory. Nanotechnology 2007, 18, 385704. [Google Scholar] [CrossRef]
  62. Vel, S.S.; Batra, R.C. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J. Sound Vib. 2004, 272, 703–730. [Google Scholar] [CrossRef]
  63. Pan, E. Exact solution for simply supported and multilayered magneto-electro-elastic plates. J. Appl. Mech. 2001, 68, 608–618. [Google Scholar] [CrossRef]
  64. Mama, B.O.; Nwoji, C.U.; Ike, C.C.; Onah, H.N. Analysis of simply supported rectangular Kirchhoff plates by the finite Fourier sine transform method. Int. J. Adv. Eng. Res. Sci. 2017, 4, 285–291. [Google Scholar]
  65. Reddy, J.N. Theory and Analysis of Elastic Plates and Shells; Taylor and Francis: Philadelphia, PA, USA, 1999. [Google Scholar]
  66. Monge, J.C.; Mantari, J.L.; Arciniega, R.A. Computational semi-analytical method for the 3D elasticity bending solution of laminated composite and sandwich doubly-curved shells. Eng. Struct. 2020, 221, 110938. [Google Scholar] [CrossRef]
  67. Ye, W.; Li, Z.; Liu, J.; Zang, Q.; Lin, G. Higher order semi-analytical solution for bending of angle-ply composite laminated cylindrical shells based on three-dimensional theory of elasticity. Thin-Walled Struct. 2019, 145, 106392. [Google Scholar] [CrossRef]
  68. Pagano, N.J. Exact solutions for composite laminates in cylindrical bending. J. Compos. Mater. 1969, 3, 398–411. [Google Scholar] [CrossRef]
  69. Pagano, N.J. Exact solutions for rectangular bidirectional composites and sandwich plates. J. Compos. Mater. 1970, 4, 20–34. [Google Scholar] [CrossRef]
  70. Pan, E.; Han, F. Exact solution for functionally graded and layered magneto-electro-elastic plates. Int. J. Eng. Sci. 2005, 43, 321–339. [Google Scholar] [CrossRef]
  71. Wu, C.P.; Liu, Y.C. A review of semi-analytical numerical methods for laminated composite and multilayered functionally graded elastic/piezoelectric plates and shells. Compos. Struct. 2016, 147, 1–15. [Google Scholar] [CrossRef]
  72. Qing, G.; Qiu, J.; Liu, Y. A semi-analytical solution for static and dynamic analysis of plates with piezoelectric patches. Int. J. Solids Struct. 2006, 43, 1388–1403. [Google Scholar] [CrossRef]
  73. Zhang, P.; Qi, C.; Fang, H.; Ma, C.; Huang, Y. Semi-analytical analysis of static and dynamic responses for laminated magneto-electro-elastic plates. Compos. Struct. 2019, 222, 110933. [Google Scholar] [CrossRef]
  74. Zhang, P.; Qi, C.; Fang, H.; Sun, X. A semi-analytical approach for the flexural analysis of in-plane functionally graded magneto-electro-elastic plates. Compos. Struct. 2020, 250, 112590. [Google Scholar] [CrossRef]
  75. Vel, S.S.; Batra, R.C. Three-dimensional analytical solution for hybrid multilayered piezoelectric plates. J. Appl. Mech. 2000, 67, 558–567. [Google Scholar] [CrossRef]
  76. Li, H.; Pang, F.; Miao, X.; Du, Y.; Tian, H. A semi-analytical method for vibration analysis of stepped doubly-curved shells of revolution with arbitrary boundary conditions. Thin-Walled Struct. 2018, 129, 125–144. [Google Scholar] [CrossRef]
  77. Pang, F.; Li, H.; Wang, X.; Miao, X.; Li, S. A semi analytical method for the free vibration of doubly-curved shells of revolution. Comput. Math. Appl. 2018, 75, 3249–3268. [Google Scholar] [CrossRef]
  78. Xavier, P.B.; Lee, K.H.; Chew, C.H. An improved zig-zag model for the bending of laminated composite shells. Compos. Struct. 1993, 26, 123–138. [Google Scholar] [CrossRef]
  79. Thai, H.T.; Kim, S.E. Analytical solution of a two variable refined plate theory for bending analysis of orthotropic Levy-type plates. Int. J. Mech. Sci. 2012, 54, 269–276. [Google Scholar] [CrossRef]
  80. Cooke, D.W.; Levinson, M. Thick rectangular plates—II: The generalized Lévy solution. Int. J. Mech. Sci. 1983, 25, 207–215. [Google Scholar] [CrossRef]
  81. Brischetto, S. Exact elasticity solution for natural frequencies of functionally graded simply-supported structures. Comput. Model. Eng. Sci. 2013, 95, 391–430. [Google Scholar]
  82. Reddy, J.N. Introduction to the Finite Element Method; McGraw-Hill Education: New York, NY, USA, 2019. [Google Scholar]
  83. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method: Solid Mechanics; McGraw-Hill: New York, NY, USA, 1967; Volume 2. [Google Scholar]
  84. Oden, J.T.; Reddy, J.N. An Introduction to the Mathematical Theory of Finite Elements; Dover Publications: New York, NY, USA, 2012. [Google Scholar]
  85. Malik, M.R.; Zang, T.A.; Hussaini, M.Y. A spectral collocation method for the Navier-Stokes equations. J. Comput. Phys. 1985, 61, 64–88. [Google Scholar] [CrossRef]
  86. Ma, X.; Huang, C. Spectral collocation method for linear fractional integro-differential equations. Appl. Math. Model. 2014, 38, 1434–1448. [Google Scholar] [CrossRef]
  87. Orszag, S.A. Numerical methods for the simulation of turbulence. Phys. Fluids 1969, 12, II-250. [Google Scholar] [CrossRef]
  88. Shu, C.; Chen, W.; Du, H. Free vibration analysis of curvilinear quadrilateral plates by the differential quadrature method. J. Comput. Phys. 2000, 163, 452–466. [Google Scholar] [CrossRef]
  89. Shu, C. Differential Quadrature and Its Application in Engineering; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  90. Shu, C.; Richards, B.E. Application of generalized differential quadrature to solve two-dimensional incompressible Navier-Stokes equations. Int. J. Numer. Methods Fluids 1992, 15, 791–798. [Google Scholar] [CrossRef]
  91. Wu, T.Y.; Liu, G. A differential quadrature as a numerical method to solve differential equations. Comput. Mech. 1999, 24, 197–205. [Google Scholar] [CrossRef]
  92. Tornabene, F. Generalized Differential and Integral Quadrature; Esculapio: Bologna, Italy, 2023. [Google Scholar]
  93. Shu, C.; Chen, W.; Xue, H.; Du, H. Numerical study of grid distribution effect on accuracy of DQ analysis of beams and plates by error estimation of derivative approximation. Int. J. Numer. Methods Eng. 2001, 51, 159–179. [Google Scholar] [CrossRef]
  94. Shu, C.; Chen, W. On optimal selection of interior points for applying discretized boundary conditions in DQ vibration analysis of beams and plates. J. Sound Vib. 1999, 222, 239–257. [Google Scholar] [CrossRef]
  95. Tornabene, F.; Viscoti, M.; Dimitri, R. Equivalent single layer higher order theory based on a weak formulation for the dynamic analysis of anisotropic doubly-curved shells with arbitrary geometry and variable thickness. Thin-Walled Struct. 2022, 174, 109119. [Google Scholar] [CrossRef]
  96. Shu, C.; Chew, Y.T.; Richards, B.E. Generalized differential and integral quadrature and their application to solve boundary layer equations. Int. J. Numer. Methods Fluids 1995, 21, 723–733. [Google Scholar] [CrossRef]
  97. Patton, A.; Antolin, P.; Dufour, J.E.; Kiendl, J.; Reali, A. Accurate equilibrium-based interlaminar stress recovery for isogeometric laminated composite Kirchhoff plates. Compos. Struct. 2021, 256, 112976. [Google Scholar] [CrossRef]
  98. Patton, A.; Antolin, P.; Kiendl, J.; Reali, A. Efficient equilibrium-based stress recovery for isogeometric laminated curved structures. Compos. Struct. 2021, 272, 113975. [Google Scholar] [CrossRef]
  99. Daniel, P.M.; Främby, J.; Fagerström, M.; Maimí, P. Complete transverse stress recovery model for linear shell elements in arbitrarily curved laminates. Compos. Struct. 2020, 252, 112675. [Google Scholar] [CrossRef]
  100. Tornabene, F.; Viscoti, M.; Dimitri, R. Static analysis of anisotropic doubly-curved shell subjected to concentrated loads employing higher order layer-wise theories. Comput. Model. Eng. Sci. 2023, 134, 1393–1468. [Google Scholar] [CrossRef]
  101. Tornabene, F.; Viscoti, M.; Dimitri, R. Static analysis of doubly-curved shell structures of smart materials and arbitrary shape subjected to general loads employing higher order theories and generalized differential quadrature method. Comput. Model. Eng. Sci. 2022, 133, 719–798. [Google Scholar] [CrossRef]
Figure 1. Coefficients for the harmonic expansion of the wave amplitudes Q s n m ± = Q 1 n m ± , Q 2 n m ± , Q 3 n m ± of external surface loads q s ± = q 1 ± , q 2 ± , q 3 ± of reference magnitudes q ¯ s ± = q ¯ 1 ± , q ¯ 2 ± , q ¯ 3 ± applied on a doubly-curved shell panel.
Figure 1. Coefficients for the harmonic expansion of the wave amplitudes Q s n m ± = Q 1 n m ± , Q 2 n m ± , Q 3 n m ± of external surface loads q s ± = q 1 ± , q 2 ± , q 3 ± of reference magnitudes q ¯ s ± = q ¯ 1 ± , q ¯ 2 ± , q ¯ 3 ± applied on a doubly-curved shell panel.
Materials 17 00588 g001
Figure 2. Further coefficients for the harmonic expansion of the wave amplitudes Q s n m ± = Q 1 n m ± , Q 2 n m ± , Q 3 n m ± of external surface loads q s ± = q 1 ± , q 2 ± , q 3 ± of reference magnitudes q ¯ s ± = q ¯ 1 ± , q ¯ 2 ± , q ¯ 3 ± applied on a doubly-curved shell panel.
Figure 2. Further coefficients for the harmonic expansion of the wave amplitudes Q s n m ± = Q 1 n m ± , Q 2 n m ± , Q 3 n m ± of external surface loads q s ± = q 1 ± , q 2 ± , q 3 ± of reference magnitudes q ¯ s ± = q ¯ 1 ± , q ¯ 2 ± , q ¯ 3 ± applied on a doubly-curved shell panel.
Materials 17 00588 g002
Figure 3. Linear correction σ i j = τ 13 , τ 23 , σ 3 of the profile of the out-of-plane stress components by means of the external load q s + = q 1 + , q 2 + , q 3 + acting at the top surface of the shell. The external pressure at the bottom surface, denoted by q s = q 1 , q 2 , q 3 , has been previously modeled as the boundary condition of the three-dimensional elasticity equation for the derivation of σ ˜ i j = τ ¯ 13 , τ ¯ 23 , σ ¯ 3 .
Figure 3. Linear correction σ i j = τ 13 , τ 23 , σ 3 of the profile of the out-of-plane stress components by means of the external load q s + = q 1 + , q 2 + , q 3 + acting at the top surface of the shell. The external pressure at the bottom surface, denoted by q s = q 1 , q 2 , q 3 , has been previously modeled as the boundary condition of the three-dimensional elasticity equation for the derivation of σ ˜ i j = τ ¯ 13 , τ ¯ 23 , σ ¯ 3 .
Materials 17 00588 g003
Figure 4. Representation of 3D FEM models of structures of various curvatures, and 2D GDQ-based models for the present semi-analytical formulation.
Figure 4. Representation of 3D FEM models of structures of various curvatures, and 2D GDQ-based models for the present semi-analytical formulation.
Materials 17 00588 g004
Figure 5. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 5. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g005
Figure 6. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 6. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g006
Figure 7. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 7. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to patch loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with c 10 = 0.25 L 1 and c 20 = 0.25 L 2 applied at s 10 , s 20 = 0.75 L 1 , 0.75 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g007
Figure 8. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 8. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g008
Figure 9. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 9. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g009
Figure 10. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 10. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g010
Figure 11. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 11. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g011
Figure 12. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 12. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g012
Figure 13. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 13. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to concentrated loads of magnitudes q ¯ 3 + = 2000 N and q ¯ 3 = 2000 N with s 10 = 0.5 L 1 and s 20 = 0.5 L 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g013
Figure 14. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 14. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g014
Figure 15. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 15. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g015
Figure 16. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 16. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a line load distributed along α 1 principal direction of magnitude q ¯ 3 + = 4.58 × 10 3 N / m . The Navier solution has been calculated setting n = m = 500 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g016
Figure 17. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 17. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated cylindrical panel subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g017
Figure 18. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 18. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g018
Figure 19. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 19. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated rectangular plate subjected to a combination of uniform and hydrostatic loads along α 1 and α 2 principal directions of magnitudes q ¯ 3 + = 2 × 10 5 N / m 2 , q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 + = 4 × 10 5 N / m 2 , respectively. Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g019
Figure 20. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with n = m = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 20. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with n = m = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g020
Figure 21. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with N ˜ = M ˜ = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 21. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with N ˜ = M ˜ = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g021
Figure 22. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with N ˜ = M ˜ = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 22. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to sinusoidal loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 with N ˜ = M ˜ = 1 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g022
Figure 23. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 23. Reconstruction of the profiles of the components of the three-dimensional displacement field vector U α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g023
Figure 24. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 24. Reconstruction of the profiles of the components of the three-dimensional strain vector ε α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g024
Figure 25. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Figure 25. Reconstruction of the profiles of the components of the three-dimensional stress vector σ α 1 , α 2 , ζ along the thickness direction of a laminated spherical panel subjected to uniform loads of magnitudes q ¯ 3 + = 7 × 10 5 N / m 2 and q ¯ 3 = 3 × 10 5 N / m 2 . Thickness plots have been provided for the point located at 0.25 L 1 , 0.25 L 2 .
Materials 17 00588 g025
Table 1. Static analysis of a rectangular plate subjected to a patch load. The percentage error in Equation (97) defines the convergence rate of the semi-analytical solution with respect to the reference simulation derived from a three-dimensional finite element simulation. The quantity w N ˜ is expressed in 10 4 m .
Table 1. Static analysis of a rectangular plate subjected to a patch load. The percentage error in Equation (97) defines the convergence rate of the semi-analytical solution with respect to the reference simulation derived from a three-dimensional finite element simulation. The quantity w N ˜ is expressed in 10 4 m .
N ˜   =   M ˜ w N ˜ e %
515.04760.16%
1015.04200.12%
1215.00770.11%
2015.00840.10%
5015.00850.10%
10015.00850.10%
12015.00850.10%
15015.00850.10%
20015.00850.10%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Tornabene, F.; Viscoti, M.; Dimitri, R. On the Importance of the Recovery Procedure in the Semi-Analytical Solution for the Static Analysis of Curved Laminated Panels: Comparison with 3D Finite Elements. Materials 2024, 17, 588. https://doi.org/10.3390/ma17030588

AMA Style

Tornabene F, Viscoti M, Dimitri R. On the Importance of the Recovery Procedure in the Semi-Analytical Solution for the Static Analysis of Curved Laminated Panels: Comparison with 3D Finite Elements. Materials. 2024; 17(3):588. https://doi.org/10.3390/ma17030588

Chicago/Turabian Style

Tornabene, Francesco, Matteo Viscoti, and Rossana Dimitri. 2024. "On the Importance of the Recovery Procedure in the Semi-Analytical Solution for the Static Analysis of Curved Laminated Panels: Comparison with 3D Finite Elements" Materials 17, no. 3: 588. https://doi.org/10.3390/ma17030588

APA Style

Tornabene, F., Viscoti, M., & Dimitri, R. (2024). On the Importance of the Recovery Procedure in the Semi-Analytical Solution for the Static Analysis of Curved Laminated Panels: Comparison with 3D Finite Elements. Materials, 17(3), 588. https://doi.org/10.3390/ma17030588

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop