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.
3. Kinematic Relations
In the present section, a generalized ESL approach is presented for the expansion up to the
-th order of the three-dimensional displacement field vector
along the thickness direction. To this end, three thickness functions
are introduced for each
-th kinematic expansion order with
. Remembering that the thickness coordinate is defined so that
for
, a generalized formulation can be derived, and the vector
of the generalized displacement field components
is introduced for each
. One gets the following relation [
16]:
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
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
, together with Murakami’s zigzag function associated to
:
When the thickness functions of Equation (9) are used, the notation
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
, the letter “Z” is used. Finally,
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
is the three-dimensional strain vector of the
-th layer, the following condensed relation can be taken into account:
In the previous relation,
is the kinematic differential operator, which is split in those denoted by
and
, which contain the derivatives with respect to the coordinate
and
, respectively. In an extended form,
can be expressed as:
On the other hand, the kinematic operators
of size
take the following aspect:
The sub-operators
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
has been expressed in terms of the generalized strain vector, for each
, denoted by
, whose definition is reported in the following:
On the other hand, the ESL kinematic operator
is defined for each
in terms of the generalized thickness functions of Equation (8):
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
-th layer of the solid, a three-dimensional linear elastic constitutive relation is considered, defined as follows [
16]:
In the previous relation,
denotes the three-dimensional stiffness matrix, whose generic component
with
relates the
-th element of the stress vector, denoted by
, to the
-th element of the strain vector
. 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
, whose axes are oriented along the material symmetry directions. For this reason, the equation reported in the following should be considered, where
and
are the vectors of the three-dimensional stress and strain components, respectively, written with respect to
material reference system, whereas
is the corresponding stiffness matrix:
In the previous equation,
with
are the elements of the matrix
. More specifically, they can be taken as equal to the three-dimensional stiffness components of the material of the -th layer, namely
. However, when the kinematic field assumption of Equation (8) neglects the stretching effect,
turn out to be the reduced elastic coefficients
, 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
is rotated by means of the transformation matrix
, as defined in Ref. [
16]. One gets:
In the present study, the material axis
is taken along the thickness direction of the shell; therefore, the equivalence
is considered. As a consequence, the rotation matrix
depends only on the angle
between
and
reference axes such that
. 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:
At this point, the previous relation can be used for the evaluation of the virtual variation
of the elastic strain energy of the shell:
The computation of the integrals of Equation (21) in the closed interval
allows one to introduce the generalized stress resultants, which are collected for each
and
in the generalized stress resultant vector
[
16]:
Finally, the higher-order ESL elastic constitutive relationship can be written in terms of
and
by substituting Equation (20) in the final expression of the virtual variation of the elastic strain energy
, as reported in Equation (22). One gets:
In the previous equation,
denotes the generalized higher-order constitutive operator, which is defined for each
and
as follows [
16]:
For the sake of clarity, the sub-matrices
,
,
,
are reported below in an extended form:
The generalized elastic coefficients of Equations (25)–(28) can be computed with the following condensed expression, setting the definitions
and
:
The quantities
with
occurring in the previous expression denote the three-dimensional elastic stiffness coefficients
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:
Accordingly, when a linear distribution of the displacement field components is considered in Equation (8), the value
is assumed, whereas in the case of higher-order theories, a unitary value
is assigned to this quantity. The value
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
and the generalized displacement field vector
. Substituting the definition of the vector
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:
The matrix
of size
, defined for each
and
is reported in the following in an extended form:
The complete expression of the coefficients
with
and
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
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
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
, is computed as the sum of the virtual work of the actions
and
with
acting at the top
and the bottom
of the shell, respectively. Referring to a doubly-curved three-dimensional solid, one gets:
In the previous relation,
and
with
are the virtual variations of the displacement field components at the top
and the bottom
surfaces of the three-dimensional solid, respectively. The application of the static equivalence principle introduces a set of generalized loads for each
kinematic expansion order, which are collected in the vector
. They are associated to the virtual variation of the vector
of the generalized displacement field components. The following relation is thus obtained:
According to the static equivalence principle, the virtual work
of Equation (35) turns out to be equal to the virtual work
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
[
16]:
In the previous equation, the quantity with denotes the thickness function associated to and , whereas with are the scaling parameters calculated at the top and the bottom 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
is considered, whereas the magnitude of the external load at issue is denoted by
In the case of a constant distribution of the external load, the distribution 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
is reported below:
where each term
with
reads as:
Introducing in Equation (40) the definition of
in terms of
, 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
[
16]:
The fundamental matrix
referred to an arbitrary
occurring in Equation (44) turns out to be of size
, 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
allows one to also derive the boundary conditions of the problem, which are applied at the edges of the rectangular physical domain
of extremes
with
. More specifically, the boundary conditions reported below are applied at the edges of the shell located at
or
:
where
and
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
or
:
In the previous relation, the fixed values of the generalized stress resultants and generalized displacement field components are denoted by and , 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:
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
and of the principal radii of curvature
. As a result, the relations reported in the following are considered:
As a consequence, the lengths
of the curvilinear parametric lines can be calculated in terms of the radii
as follows:
In the case of a cylindrical surface with
and
, the quantities
read as follows:
When a rectangular plate is studied, the principal curvatures are null, namely
, 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
. The displacement field components
and
are thus expressed for each
-th kinematic expansion order with
as follows:
In the previous equation, the quantities
and
denote the wave number of the solution along
and
, respectively, whereas the quantities
and
are the wave amplitudes associated to each wave number
. The total number of waves along
and
is denoted by
and
, respectively. According to the Navier’s method, these quantities are assumed as
. 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
and
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
. The wave amplitudes of the external loads, defined for each wave number
, are denoted by
and
. 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
and
introduced in Equation (38) which are on the shell reference surface, are expanded for each
with trigonometric series:
being
and
the amplitudes of the generalized external actions
associated to the
-th kinematic expansion order, with
. Introducing Equations (55) and (56) in the static equivalence principle of Equation (37), the following definitions of the generalized amplitudes
of the actions
are obtained in terms of the wave amplitudes
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
occurring in Equation (19) is selected so that
or
. 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 wave amplitude of the displacement field components is derived for each
from the following expression:
where the vector
collects the amplitudes of the generalized external actions of the
-th expansion order. Furthermore, the quantity
is the fundamental matrix related to the wave numbers
, defined for each
. In a more expanded form, the previous relation becomes:
The complete expression of the coefficients
of the fundamental matrix
with
, is reported below for each
:
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
, the generalized coefficients
of Equation (29) are calculated for each
from the following relation:
In this case, the fundamental coefficients
, 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
for the case of a cylindrical surface.
Finally, in the case of a rectangular plate with
, the generalized coefficients
are derived as follows:
The expression of the coefficients
of the fundamental matrix
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
, the location of the generic element
of the grid at issue is derived as follows:
where
is the total number of discrete points. As stated previously, the GDQ method allows one to evaluate the derivative of a generic
-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
defined in the closed interval
with
, its
-th order derivative in a point
with
is thus calculated with the expression reported in the following [
16]:
In the previous relation, the quantity
with
denotes the values assumed by the function in the discrete grid, whereas
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
. The GDQ weighting coefficients
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 and denote the first order derivatives of the Lagrange polynomials at the points and , respectively. On the other hand, the definition with should be introduced, being 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
with
, of a smooth function
with
can be evaluated as follows [
16]:
As can be seen, the definition interval of
is discretized with a grid of
points. The GIQ coefficients
with
are collected in the matrix
of size
. The coefficients at issue are derived from the GDQ shifted coefficients of the first order derivative, denoted by
with
, which are defined as follows, setting
:
The shifted coefficients of Equation (68) are collected in the matrix
, whose size is
. It can be shown that the matrix of the GIQ coefficients is the inverse of the matrix
:
When the numerical integration is restricted to an arbitrary interval
, the domain
becomes a parent interval. For this reason, the coefficients
of Equation (67) are transformed into those
by means of the following GIQ mapping technique:
In this way, the integral of
restricted to the interval
are evaluated as follows [
16]:
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
along the physical domain and
, the three-dimensional equilibrium equations assume the following aspect [
16]:
where the parameters
and
are written in an extended form as follows:
At this point, a two-dimensional grid is extracted from the physical domain made of 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
points is defined in each interval
of length
related to an arbitrary
-th layer of the stacking sequence, setting
. The generic element of this grid is denoted by
for
, with
. The quantity
is defined from an arbitrary distribution within the dimensionless interval
as follows:
Finally, the elements are arranged in the vector of size , defined in each -th lamina. At this point, a new vector of size is introduced, which contains all the discrete points belonging to the interval that are located in the thickness direction. To this end, the index is introduced, defined as . As a consequence, the vector is arranged in the global vector of size .
The through-the-thickness displacement field profile can be evaluated for each point
of the reference surface, remembering the kinematic assumption of Equation (8). The relation reported in the following is thus considered for each
:
In the same way, from Equation (14) the profiles are derived of the three-dimensional strain components of the vector
for each point
of the reference surface of the shell:
The quantity
with
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
for
, the distribution of the membrane stresses
and
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:
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
parametric lines. These derivations are performed numerically by means of the GDQ method of Equation (65). One gets:
where the notations
and
mean that the GDQ rule is applied along
and
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
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
and
are evaluated from the first two three-dimensional equilibrium Equation (72) of a doubly-curved shell, which are written in each
-th layer of the shell in a discrete form as follows:
The parameters
and
, dependent on the membrane stresses
and
, 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
and
. In the same way, at the top surface of the shell, the shear stress profile must guarantee equilibrium with the external loads
and
:
The boundary conditions introduced in the previous equation are written below in discrete form:
At this point, a numerical solution is found for
, 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
-th layer with
are assessed starting from the results obtained for
. Note that the discrete points associated to the indexes
and
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
and
, derived numerically, are corrected by adding a linear term dependent on the thickness coordinate
[
16], as shown in
Figure 3:
At this point, the derivative of
and
shear stresses with respect to
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
:
where the parameter
accounts for the following expression:
Two possible sets of boundary conditions are considered for the derivation of the numerical solution of Equation (87), setting
and
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:
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
, the following boundary condition is considered for the numerical assessment of the equilibrium Equation (87) in the remaining laminae, namely for
:
The second boundary condition of Equation (90) is applied by means of a linear correction [
16] of the profile of the normal stress
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 updated profile of the out-of-plane strain components
and
. To this end, the following relation is considered at each point
of the three-dimensional solid:
Starting from the previous equation, the expression of the quantities and 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 and 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
, whose engineering constants, obtained from Ref. [
101], are reported in the following:
In the previous relation, the quantities
denote the Poisson’s coefficients of the orthotropic materials, whereas
and
are the elastic moduli and the shear moduli, respectively. They are related to the three-dimensional stiffness constants
with
of Equation (18) according to the procedure extensively detailed in Ref. [
16]. The stiffness constants
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
, the following aspect:
In the same way, the mechanical properties of the graphite-epoxy soft10
can be written as:
In each simulation, the lamination scheme consists of three layers of thicknesses and with material orientation angles equal to and . Finally, the recovery procedure has been applied setting discrete points in the thickness direction for each -th layer. Note that the computational cost of the present semi-analytical formulation is here intended as the number of terms, denoted by and , 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
and
made of three layers of graphite-epoxy with properties as in Equation (94) of thicknesses
and
subjected to two different patch loads, one at the top surface and the other at the bottom surface with magnitudes
and
, 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
. 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
:
where
denotes the vertical deflection of the central point of the three-dimensional solid derived from Equation (75), while
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
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
discrete points, following the CGL distribution. The present semi-analytical theory has been used with the EDZ4 displacement field assumption, setting
. As can be seen from
Figure 5, very accurate results are provided in terms of the displacement field components
. 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
, especially in the central layer. Furthermore, very accurate results are obtained for in-plane quantities
. 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
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
and directed along
principal direction is applied at the top surface, whereas the second load of magnitude
, applied at the bottom surface, is directed along
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
have been collected in
Figure 8,
Figure 9 and
Figure 10.
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
, 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 is considered with and . Two externally concentrated loads of magnitude and are applied at the top and bottom surfaces, respectively. The position parameters are and , 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
, 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
, 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
distributed along
principal direction, located at
. The thickness plots are evaluated at
. The GDQ numerical solution has been calculated with
, accounting for FSDT, TSDT, and EDZ4 theories, whereas the semi-analytical Navier solution is derived with the EDZ4 kinematic assumption setting
. 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
, 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
and
are considered, which are denoted by H1 and H2, respectively. In addition, a uniform pressure (U) of magnitude
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:
For each load case, the thickness plots are evaluated
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
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 , whose physical domain is defined so that and . 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 and with . 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
. 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
within the physical domain, where the quantities
and
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
and
. Furthermore, when the central layer stiffness is reduced, the vertical deflection
increases. Similar considerations can be made for the strain components of
Figure 21, where
and
assume in the central lamina a non-linear profile in the case of C2 and C3, whereas the value of
and
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 (
), the linear static response of the same structure is derived for the case of uniform transverse loads
and
applied at the top and bottom surfaces, respectively. In this case, the results obtained with the present semi-analytical theory have been derived setting
, taking into account the EDZ4 higher-order theory. The thickness plots are provided for the point
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.