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

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.


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 semianalytical 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 twodimensional 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.

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 , where α 0 i < α 1 i 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: 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]: 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: Note that the symbol r ,ij = ∂ 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 ni = 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: The Lamè parameters A i = A 1 , A 2 are used for the computation of the infinitesimal arch lengths ds i = ds 1 , ds 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: 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: 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]:

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 along the thickness direction.To this end, three thickness functions 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 3 is introduced for each τ = 0, . . ., N + 1.One gets the following relation [16]: The ESL formulation of the previous relation allows one to derive a generalized twodimensional 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: 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] is the three-dimensional strain vector of the k-th layer, the following condensed relation can be taken into account: In the previous relation, D is the kinematic differential operator, which is split in those denoted by D ζ and D Ω , which contain the derivatives with respect to the coordinate ζ and α 1 , α 2 , respectively.In an extended form, D ζ can be expressed as: On the other hand, the kinematic operators Ω of size 9 × 3 take the following aspect: The sub-operators D Ω are written in an extended form as: 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]: 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 , whose definition is reported in the following: 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):

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]: In the previous relation, E (k) denotes the three-dimensional stiffness matrix, whose ij 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 consti- tutive 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 threedimensional stress and strain components, respectively, written with respect to O ′ α 1 α 2 α 3 material reference system, whereas E (k) is the corresponding stiffness matrix: In the previous equation, E ij 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 ij .However, when the kinematic field assumption of Equation ( 8) neglects the stretching effect, E (k) ij turn out to be the reduced elastic coefficients Q (k) ij , which are calculated from the three-dimensional ones according to what is exerted in Ref. [16].

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]: 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 The virtual work of external actions, denoted by δL , is computed as the sum of the virtual work of the actions q ia and q ia 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: In the previous relation, U 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 (τ) 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: According to the static equivalence principle, the virtual work δL of Equation ( 35) turns out to be equal to the virtual work δL e of Equation (36): 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 3 [16]: In the previous equation, the quantity 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 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: In the previous relation, Ω denote the equilibrium operators, which can be expressed with a matrix notation as follows: An extended version of the quantities where each term D * α i Ω with i = 1, 2, 3 reads as: 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]: 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]: 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 α More specifically, the boundary conditions reported below are applied at the edges of the shell located at α 1 = α 0 1 or α 1 = α 1 1 : or u or u where , and 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 = α 0 2 or α 2 = α 1 2 : or u or u In the previous relation, the fixed values of the generalized stress resultants and generalized displacement field components are denoted by 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:

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: 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: In the case of a cylindrical surface with k n1 = 0 and R 2 = R, the quantities L 1 , L 2 read as follows: When a rectangular plate is studied, the principal curvatures are null, namely k n1 = k n2 = 0, therefore, the lengths of the parametric lines are calculated from the following expression: The semi-analytical Navier solution accounts for the harmonic distribution of the unknown field variable within the physical domain (s 2 and u (τ) 3 are thus expressed for each τ-th kinematic expansion order with τ = 0, . . ., N + 1 as follows: 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 2nm and U 3nm 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 Ñ and M, respectively.According to the Navier's method, these quantities are assumed as Ñ = M = ∞.It can be seen that the harmonic expansion of Equation ( 53) respects the following boundary conditions along the edges of the physical domain: 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: with Ñ = M = ∞.The wave amplitudes of the external loads, defined for each wave number n, m, are denoted by In Figures 1 and 2 we report the expressions of the wave amplitudes of some load shapes of practical interest., , , , , , s q q q q ± ± ± ± = applied on a doubly-curved shell panel.
In addition, the material orientation angle ( ) 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 ( ) ( ) ( ) ( ) of the of external surface loads q 3 of reference magnitudes q The generalized external actions 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: being 3snm the amplitudes of the generalized external actions 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 are obtained in terms of the wave amplitudes Q 3snm of the surface loads applied at the top and bottom surfaces of the shell: 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: 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 T of the wave amplitude of the displacement field components is derived for each τ = 0, . . ., N + 1 from the following expression: where the vector T collects the amplitudes of the generalized external actions of the τ-th expansion order.Furthermore, the quantity L(τη) nm 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: The complete expression of the coefficients L(τη)α i α j ijnm of the fundamental matrix L(τη) nm , with i, j = 1, 2, 3, is reported below for each τ, η = 0, . . ., N + 1: [11] [11]α 3 α 3 33(00) (61) 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 of Equation ( 29) are calculated for each τ, η = 0, . . ., N + 1 from the following relation: In this case, the fundamental coefficients L(τη)α i α j ijnm , 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 ijnm for the case of a cylindrical surface.
Finally, in the case of a rectangular plate with are derived as follows: The expression of the coefficients L(τη)α i α j ijnm of the fundamental matrix L(τη) nm can be found for the case of a rectangular plate in Appendix B.

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: 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 ∈ R, 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]: 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 ς (n) ij 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 ς (n) ij 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: 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 ς (0) ij = δ ij with i, j = 1, . . ., I Q should be introduced, being δ ij 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]: As can be seen, the definition interval of f is discretized with a grid of I Q points.The GIQ coefficients w ik , w jk 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 ς (1) ij with i, j = 1, . . .I Q , which are defined as follows, setting ε = 1 × 10 −10 : 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) : 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 w1I Q k of Equation ( 67) are transformed into those w by means of the following GIQ mapping technique: In this way, the integral of f = f (x) restricted to the interval [a, b] are evaluated as follows [16]:

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]: where the parameters a (k) , b (k) and c (k) are written in an extended form as follows: At this point, a two-dimensional grid is extracted from the physical domain s 0 1 , s 1 1 × s 0 2 , s 1 2 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 Finally, the elements ζ m are arranged in the vector 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 The through-the-thickness displacement field profile can be evaluated for each point s 1i , s 2j 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 : In the same way, from Equation ( 14) the profiles are derived of the three-dimensional strain components of the vector ε T for each point s 1i , s 2j of the reference surface of the shell: The quantity ε , s 2j 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 ij for i, j = 1, . . ., 6, the distribution of the membrane stresses σ tational grid the three-dimensional strain components of Equation ( 76), remembering the hypotheses made for the derivation of the Navier closed-form solution: 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: (78) where the notations ς 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 τ (k) 13 and τ (k) 23 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: The parameters a (k) and b (k) , dependent on the membrane stresses σ 12 , are written in a discrete form as follows: 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 (+) The boundary conditions introduced in the previous equation are written below in discrete form: τ (1) 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: 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 τ   , , s q q q q + + + + = acting at the top surface of the shell.The external pressure at the bottom surface, denoted by ( ) ( ) ( ) ( ) , , s q q q q − − − − = , has been previously modeled as the boundary condition of the three-dimensional elasticity equation for the derivation of , , At this point, the derivative of ( ) , s s can be evaluated with the GDQ method as follows: acting at the top surface of the shell.The external pressure at the bottom surface, denoted by q At this point, the derivative of τ (k) 13 and τ (k) 23 shear stresses with respect to s 1 , s 2 can be evaluated with the GDQ method as follows: 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 σ (k) 3 : where the parameter c (k) (ijm) accounts for the following expression: Two possible sets of boundary conditions are considered for the derivation of the numerical solution of Equation ( 87), setting q (+) 3 and q the value of the external actions oriented along the normal direction applied at the top and bottom surfaces of the shell, respectively: In discrete form, the last two relations become: (90) 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: The second boundary condition of Equation ( 90) is applied by means of a linear correction [16] of the profile of the normal stress σ (k) 3 derived from Equation (87): 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 , 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 throughthe-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 , 0.25 L L .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 :  In the same way, the mechanical properties of the graphite-epoxy soft10 ρ (k) = 145 kg/m 3 can be written as: 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 Ñ 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-thethickness 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 % : where w Ñ 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 Ñ = 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 Figures 5-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 Ñ = 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 Figures 8-10.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 ( ) For each load configuration, the semi-analytical solution has been calculated with . 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 3 U , reported in Figure 8, whereas the EDZ4 theory perfectly matches the threedimensional predictions of all quantities in each load case.With particular reference to both in-plane and out-of-plane strain and stress profiles in Figures 9 and 10, it should be noted that for each case, the loading conditions are perfectly respected.Furthermore, the          For each load configuration, the semi-analytical solution has been calculated with Ñ = 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 Figures 9 and 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 α 0 2 , α 1 2 = (−π/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 Figures 11-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 Ñ = 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  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 Ñ = M = 500.Thickness plots of displacement field, strain, and stress components are reported in Figures 14-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.For each load case, the thickness plots are evaluated ( ) 0.25 , 0.25 L L and collec , and a 3D finite element reference solution has been derived commercial software.In addition, a GDQ-based numerical solution has been eva with the EDZ4 theory, which matches the 3D FEM predictions.The semi-ana solution is evaluated with the EDZ4 displacement field assumption, set each load case.The results are in line with the reference solution in terms of displac field components, as can be seen in Figure 17.It should be noted that the se lamination scheme presents some zigzag effects, especially for case C5.The adopt Murakami's zigzag function in the present model allows one to consider the piec inclination of the profile of the displacement field components.Similar consideratio made for the three-dimensional strain components of Figure 18, where the 3D reference solution is perfectly predicted in a reduced amount of time by the higher semi-analytical model, and a good level of accuracy is also reached in the central r of the structure.In order to reduce the computational effort of the simulation, the r of more complicated load cases like C3 are obtained from the algebraic sum of C1 a simulations.In the same way, load case C5 is derived from the sum of C3 and C4 consequence, the results of the simulations referred to in C3 and C5 can be effic obtained as an algebraic sum of the profiles of all kinematic and mechanical qua recovered previously in the corresponding two-dimensional semi-analytical solutio shown in Figure 19 in the case of the three-dimensional stress components.The rec procedure is not applied in C3 and C5 because the equilibrium equations in the thic For each load case, the thickness plots are evaluated (0.25 L 1 , 0.25 L 2 ) and collected in Figures 17-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 Ñ = 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 threedimensional 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.
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 α 0 1 , α 1 1 = (5π/12, 7π/12) and α 0 2 , α 1 2 = (−π/9, π/9).Three different lamination schemes are considered made of two external layers of graphiteepoxy ( 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 Ñ = 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 293,847 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 Figures 20-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 ( Ñ = 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 Ñ = 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 Figures 23-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 throughthe-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 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 1 U and 2 U .Furthermore, when the central layer stiffness is reduced, the vertical deflection 3 U 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.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   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 threedimensional model are well predicted when higher-order thickness functions are considered in the two-dimensional semi-analytical solution.

Conclusions
In the present manuscript, an efficient two-dimensional semi-analytical mode been presented for the evaluation of the static response of curved laminated pa subjected to arbitrary loads.The fundamental governing equations have been wr 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.

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 uniformlydistributed patch and hydrostatic loads are applied to the panel, the convergence of results is seen for Ñ = M = 150, while further terms are required in the case of concentrated and line loads, namely Ñ = 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 semianalytical 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.
2 denotes the thickness function associated to U (+) i and U (−) i , whereas H (j) 1 , H (j) 2 with j = 1, 2 are the scaling parameters calculated at the top (ζ = h/2) and the bottom (ζ = −h/2) surfaces of the shell.

Figure 1 .
Figure 1.Coefficients for the harmonic expansion of the wave amplitudes ( )

Figure 1 .
Figure 1.Coefficients for the harmonic expansion of the wave amplitudes Q

Figure 2 .
Figure 2.Further coefficients for the harmonic expansion of the wave amplitudes( )

Figure 2 .
Figure 2. Further coefficients for the harmonic expansion of the wave amplitudes Q

Figure 3 .
Figure 3. Linear correction σ ij = τ 13 , τ 23 , σ 3 of the profile of the out-of-plane stress components by means of the external load q

3 ,
has been previously modeled as the boundary condition of the three-dimensional elasticity equation for the derivation of σij = τ 13 , τ 23 , σ 3 . τ

Figure 4 .
Figure 4. Representation of 3D FEM models of structures of various curvatures, and 2D GDQ-based models for the present semi-analytical formulation.The first example consists of a simply-supported rectangular plate of dimensions the vertical deflection of the central point of the three-dimensional solid derived from Equation (75), while FEM 0.00150236 m w =is the corresponding value derived from the finite element simulation.

Figure 4 .
Figure 4. Representation of 3D FEM models of structures of various curvatures, and 2D GDQ-based models for the present semi-analytical formulation.

Figure 5 .
Figure 5. Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

. 1 α
Thickness plots have been provided for the point located at next example takes into account the same rectangular panel subjected to hydrostatic loads.More specifically, the first load of magnitude principal direction is applied at the top surface, whereas the second load of magnitude the bottom surface, is directed along 2 α

Figure 8 .U
Figure 8. Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

Figure 9 .
Figure 9. Reconstruction of the profiles of the components of the three-dimensional strain vector ( ) 1 2 , , α α ζ ε q − = − ⋅ , respectively.Thickness plots have been provided for the point located at ( )

Figure 11 .U
Figure 11.Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

Figure 12 .
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 ).

55 Figure 13 .σ
Figure 13.Reconstruction of the profiles of the components of the three-dimensional stress vector ( ) 1 2 , , α α ζ σ GDQ numerical solution has been calculated with

55 Figure 14 .U
Figure 14.Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

Figure 14 .) 3 = 55 Figure 15 ..
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

Figure 15 .) 3 = 3 = − 4 × 3 = − 2 × 98 )Figure 16 .
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 ).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 (C1) → H1 Case 02 (C2) → H2 Case 03 (C3) → H1 + H2 Case 04 (C4) → U Case 05 (C5) → H1 + H2 + U

Figure 16 .) 3 =
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 17 .U
Figure 17.Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

Figure 18 .
Figure 18.Reconstruction of the profiles of the components of the three-dimensional strain vector ( ) 1 2 , , α α ζ ε different lamination schemes are considered made of two external layers of graphite-epoxy(94), whereas the central core consists of graphite-epoxy (94), graphiteepoxy 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 solutions are provided, developed with 3D finite element solution and a two-dimensional GDQ-based formulation, accounting for the EDZ4 displacement field theory.

Figure 20 .U
Figure 20.Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

.
Thickness plots have been provided for the point located at ( )

Figure 23 .U
Figure 23.Reconstruction of the profiles of the components of the three-dimensional displacement field vector ( ) 1 2 , , α α ζ U

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 Ñ is expressed in 10 −4 m.